pounce_presolve/

lib.rs

//! Algorithmic NLP preprocessing exposed as a composable TNLP wrapper.
//!
//! Tracks pounce issue #20.
//!
//! * **Phase 0** — scaffolding, options table, no-op identity path.
//! * **Phase 1** — Andersen-style bound tightening against linear rows.
//! * **Phase 2** — redundant linear-constraint removal: rows whose
//!   activity interval is implied by the (possibly Phase-1-tightened)
//!   variable box are dropped from the problem the solver sees, then
//!   reinstated with `λ=0` when forwarding `eval_h` / `finalize_solution`
//!   to the inner TNLP.
//!
//!   **Dual attribution (issue M24).** A row that was *itself* the
//!   reason a bound got tightened in Phase 1 — e.g. `x ≥ 2` tightening
//!   `x_l = 2` — then has an activity interval flush against its own
//!   bound, so Phase 2 sees it as redundant and drops it. If that bound is
//!   active at the solution the interior-point method reports the dual on
//!   the *variable bound* multiplier (`z_l`/`z_u`) against a bound that did
//!   not exist in the original problem. For a **singleton** row Phase 2
//!   dropped, `finalize_solution` hands that multiplier back to the row
//!   (`λ = (z_u − z_l)/a`, the bound's side zeroed), and
//!   `get_starting_point` does the inverse, so the duals agree with a
//!   no-presolve solve and a `.sol` round-trips; see `SingletonBound`. A
//!   declared bound tied with the row keeps the multiplier, row `λ = 0`
//!   (a degenerate pair: `presolve=no`'s interior point splits it roughly
//!   1:1, and both are valid KKT points; the rule is "the declared bound is
//!   the incumbent", pinned by `tests/presolve_derived_bound_duals.rs`).
//!
//!   **Bounds derived from multi-variable rows** (Phase 1 / FBBT, e.g.
//!   `x0 = 2·x1, x0 ≥ 4` ⇒ `x1 ≥ 2`, or `x0 = 0.5` then `x0 + x1 ≤ 2` ⇒
//!   `x1 ≤ 1.5`) are handled by a general postsolve recovery rather than
//!   per-derivation provenance: `finalize_solution` zeroes every multiplier
//!   on a *derived* side (presolve's bound strictly tighter than the
//!   declared one) that the original model does not have active at `x*`,
//!   and re-solves for the original problem's multipliers on the active
//!   rows and declared bounds the leftover stationarity residual reaches —
//!   a sign-constrained least squares, dense over that connected component,
//!   by a proximal-point iteration started from the reduced solve's own
//!   multipliers (`dual_recovery`). Where the original multipliers are
//!   unique this is exactly the `presolve=no` answer; where they are not it
//!   is a valid KKT multiplier set of the original problem. It runs only
//!   when a derived side carries a multiplier, and a component past the
//!   dense cap (600 multipliers / 5000 variables) or a solve that does not
//!   close stationarity leaves the reduced attribution in place with a
//!   warning. Inbound, `get_starting_point` moves a Phase-2-dropped row's
//!   supplied multiplier onto the reduced bounds of the variables it
//!   touches (`a_ij·λ_i`), so the reduced start is stationary where the
//!   supplied one was. A `.sol` round-trips through
//!   `warm_start_init_point=yes max_iter=0` under any warm-start push: the
//!   recovery above, the singleton hand-back, Phase 0's frames and
//!   `linear_eq_elim`'s sweep all re-derive multipliers from stationarity
//!   at the returned `x`, which the default push moved 1e-3 off its bounds,
//!   so they would be `O(push)` off the supplied duals. Instead each layer
//!   records the warm start it served and, when the solver hands back the
//!   served multipliers unchanged, reports the ones it was handed — its
//!   inbound map's exact inverse on that map's image (`warm_echo`). A cold
//!   solve records none, so its postsolve is the recovery alone.
//!   Primal point and objective are unaffected throughout.
//! * **Phase 3** — structural LICQ check on the surviving equality
//!   rows. Verdict is published via [`PresolveTnlp::licq_verdict`].
//! * **Phase 4** — bound-multiplier warm-start hints for variables
//!   whose bounds were tightened by Phase 1. Hints are emitted on
//!   `init_z` and exposed via [`PresolveTnlp::z_warm_starts`].
//! * **Phase 5** — sensitivity-aware passthrough: projects
//!   user-supplied constraint metadata and scaling through the row
//!   reduction on the way in, and expands outer→inner on the way out
//!   in `finalize_metadata`.
//! * **Phase 6** — linear-equality variable elimination
//!   (`presolve_linear_eq_reduction`, issue #487). The only phase that
//!   removes *columns*: singleton and two-variable linear equality rows
//!   determine variables, which are substituted away and recovered at
//!   `finalize_solution`. It lives in its own wrapper
//!   ([`linear_eq_elim::LinearEqElimTnlp`]) stacked **outside**
//!   [`PresolveTnlp`], so Phases 0–5 keep working in the variable space
//!   they were written against, and it is the outer layer that
//!   re-presents the reduced one. Off by default.

#![cfg_attr(test, allow(clippy::unwrap_used, clippy::expect_used))]

use std::cell::RefCell;
use std::rc::Rc;

use pounce_common::exception::SolverException;
use pounce_common::options_list::OptionsList;
use pounce_common::reg_options::RegisteredOptions;
use pounce_common::tolerance::is_negligible;
use pounce_common::types::{Index, Number, lower_bound_present, upper_bound_present};
use pounce_nlp::expression_provider::ExpressionProvider;
use pounce_nlp::tnlp::{
    BoundsInfo, IndexStyle, InfeasibilityProof, IpoptCq, IpoptData, IterStats, Linearity, MetaData,
    NlpInfo, ScalingRequest, Solution, SparsityRequest, StartingPoint, TNLP,
};

pub mod auxiliary;
pub mod block_solve;
pub mod bound_tighten;
pub mod btf;
pub mod components;
pub mod coupling;
pub mod diagnostics;
mod dual_recovery;
pub mod dulmage_mendelsohn;
pub mod fbbt;
pub mod incidence;
pub mod inequality_projection;
pub mod licq;
pub mod linear_eq_elim;
pub mod linear_eq_plan;
mod live_lift;
pub mod matching;
pub mod options;
pub mod reduction_frame;
pub mod redundant;
pub mod trivial_elim;
pub mod warm;
mod warm_echo;

pub use block_solve::{
    BlockEquations, BlockSolveError, BlockSolveOptions, BlockSolveOutcome, BlockSolver,
    DampedNewtonSolver,
};
pub use bound_tighten::{INF_BOUND, LinearRow, TightenReport, tighten_bounds};
pub use btf::{BlockTriangularBlock, BlockTriangularForm};
pub use components::{SquareComponent, SquareComponents};
pub use coupling::{AuxiliaryCouplingClass, classify_block, objective_gradient_support};
pub use diagnostics::{AuxiliaryPreprocessingDiagnostics, AuxiliaryRejectionReason};
pub use dulmage_mendelsohn::{DMPart, DulmageMendelsohnPartition};
pub use incidence::{EqualityIncidence, InequalityIncidence, ProbeView};
pub use licq::{EqRow, LicqVerdict, licq_check};
pub use linear_eq_elim::{FullSolution, LinearEqElimTnlp, recover_dropped_multipliers};
pub use linear_eq_plan::{
    ElimStep, EliminationPlan, LinearEqElimReport, PlanConfig, PlanInput, VarRecovery, build_plan,
};
pub use options::{AuxiliaryCouplingPolicy, LicqAction, PresolveOptions, register_options};
pub use reduction_frame::{ReductionFrame, ReductionStack};
pub use redundant::find_redundant_rows;
pub use warm::{
    PresolveFingerprint, PresolveMap, ProjectedWarm, WarmPoint, WarmProjectionReport, clamp_seed,
    compute_fingerprint, project_warm_point, project_warm_point_full,
};

/// Errors that can arise while building a presolved TNLP.
#[derive(Debug)]
pub enum PresolveError {
    OptionsError(SolverException),
}

impl std::fmt::Display for PresolveError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::OptionsError(e) => write!(f, "presolve options error: {e}"),
        }
    }
}

impl std::error::Error for PresolveError {}

impl From<SolverException> for PresolveError {
    fn from(e: SolverException) -> Self {
        Self::OptionsError(e)
    }
}

/// Stack Phase 6 on top of an already-built presolve wrapper, when the
/// option asks for it.
///
/// Phase 6 goes *outside*, not inside: it is the one pass that changes the
/// variable count, and Phases 0–5 are written against a fixed column
/// space. Putting it on top also means it sees the row set Phase 2 already
/// reduced, and its `finalize_solution` hands a full-length `x` back down
/// before `PresolveTnlp` lifts the row duals — so each layer only has to
/// undo its own reduction.
fn stack_linear_eq_elim(
    wrapped: Rc<RefCell<dyn TNLP>>,
    opts: PresolveOptions,
) -> Rc<RefCell<dyn TNLP>> {
    if !opts.linear_eq_reduction {
        return wrapped;
    }
    Rc::new(RefCell::new(LinearEqElimTnlp::new(wrapped, opts)))
}

/// Top-level entry: returns a TNLP wrapping `inner` with whatever
/// presolve passes the option table has enabled. When the master
/// switch is off, returns `inner` unchanged.
pub fn wrap_with_presolve(
    inner: Rc<RefCell<dyn TNLP>>,
    opts: PresolveOptions,
) -> Result<Rc<RefCell<dyn TNLP>>, PresolveError> {
    if !opts.enabled || inner.borrow().is_presolve_wrapper() {
        return Ok(inner);
    }
    let wrapped: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(PresolveTnlp::new(inner, opts)));
    Ok(stack_linear_eq_elim(wrapped, opts))
}

/// Same as [`wrap_with_presolve`] but also installs an
/// [`ExpressionProvider`] so passes like FBBT (issue #62) can see
/// constraint expression trees. Callers who have the concrete inner
/// TNLP type (`pounce-cli` with `NlTnlp`) should prefer this; the
/// plain `wrap_with_presolve` leaves `presolve_fbbt` as a silent
/// no-op.
pub fn wrap_with_presolve_provider(
    inner: Rc<RefCell<dyn TNLP>>,
    expr_provider: Rc<RefCell<dyn ExpressionProvider>>,
    opts: PresolveOptions,
) -> Result<Rc<RefCell<dyn TNLP>>, PresolveError> {
    if !opts.enabled || inner.borrow().is_presolve_wrapper() {
        return Ok(inner);
    }
    let wrapped: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(
        PresolveTnlp::with_expression_provider(inner, expr_provider, opts),
    ));
    Ok(stack_linear_eq_elim(wrapped, opts))
}

/// Convenience: read the `presolve_*` keys out of an `OptionsList`
/// and call [`wrap_with_presolve`].
pub fn wrap_from_options(
    inner: Rc<RefCell<dyn TNLP>>,
    options: &OptionsList,
) -> Result<Rc<RefCell<dyn TNLP>>, PresolveError> {
    let opts = PresolveOptions::from_options_list(options)?;
    wrap_with_presolve(inner, opts)
}

/// Which accepting test the witness-refutation gate uses to decide that a
/// sampled point satisfies a row.
///
/// The two forms exist because the question the witness answers is not always
/// the same question. Both are *accepting* tests, so both must fail closed —
/// when in doubt, accept the point, withdraw the verdict, keep the proof
/// unclaimed.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum WitnessRule {
    /// **Default, and the only rule for any path where the solve can run.**
    ///
    /// Accepts a row when the residual is negligible at the row's *live*
    /// magnitude through the clamped form `tol * max(scale, 1)` — i.e. exactly
    /// what the solver's own acceptance test would wave through. This is the
    /// #380 rule: never certify what the solver itself would accept as
    /// feasible, or the same model reports "proved infeasible" with presolve on
    /// and `Solve_Succeeded` with it off.
    #[default]
    SolverAcceptance,
    /// **Only for paths where the solve provably cannot run**, so no
    /// `Solve_Succeeded` counterfactual exists to contradict (today: the
    /// too-few-degrees-of-freedom gate, gh#391).
    ///
    /// Accepts a row when the residual is negligible relative to the row's
    /// **declared** magnitude, `max(|g_l|, |g_u|)` over its *finite* bounds,
    /// with no absolute clamp. Declared, not live: the live value moves with
    /// the sampled point, while the bounds are the magnitude the modeller wrote
    /// the row in — the same "declared, not live" pattern
    /// `fbbt_infeasibility_survives_margin` uses for its per-row margin.
    ///
    /// Why the clamp has to go here: `tol * max(scale, 1)` reinstates an
    /// absolute floor once the row's magnitude drops below 1, so multiplying
    /// every row of an infeasible model by `1e-12` — which changes the feasible
    /// set not at all — makes every point of the box "satisfy" every row and
    /// withdraws a scale-free bound-propagation proof. That is the whole of
    /// gh#391: `s*x == 0.2*s` with `s*x == 0.8*s` crosses by `0.6` at every `s`,
    /// yet the verdict flipped at `s <= ~3e-8`.
    ///
    /// Why it is still safe: the clamp's job is to not demand more precision
    /// than the solver promised, because a solver converges to *absolute*
    /// residuals. On this path the solver never produces a point at all, so
    /// there is no converged residual to be compatible with — the alternative
    /// to the proof is a structural error, not a solution.
    ///
    /// The `b = 0` hazard is handled explicitly and fails closed: a homogeneous
    /// row (`g_l = g_u = 0`, or a row with no finite bound at all) has no
    /// declared magnitude, which would make the relative test unsatisfiable by
    /// construction — `viol == scale` for any nonzero float noise — and a
    /// genuinely feasible point would fail to refute. Such a row keeps
    /// [`WitnessRule::SolverAcceptance`]'s clamped form, i.e. the absolute
    /// floor.
    DeclaredRowRelative,
}

/// Cached, reduced view of the problem after presolve passes have
/// run. Exposed for inspection from integration tests.
/// How far a model must be from satisfiable before an emptiness detection is
/// promoted from a *detection* to a *proof*.
///
/// The rule: **never certify what the solver itself would accept as feasible.**
/// POUNCE returns `Solve_Succeeded` for a violation up to its feasibility
/// tolerance `tol`, so certifying anything smaller makes the same model report
/// "proved infeasible" with presolve on and "solved" with it off — a
/// contradiction no caller can reconcile.
///
/// Presolve's internal emptiness tests are far tighter than this: bound
/// propagation uses a `1e-12` crossing margin and FBBT's tests are exact
/// (`is_empty`, `new_lo > new_hi` — no margin at all). Those are right for
/// deciding whether to keep propagating, and wrong for the word "proved".
///
/// Two models forced this. `x >= 0.1 + 0.2` with `x <= 0.3` is infeasible by
/// `5.5e-17` — pure binary-float noise — yet was reported proved infeasible.
/// And `x >= 1 + 1e-11` with `x <= 1` is solved to `Solve_Succeeded` on the
/// default path while presolve called it proved infeasible.
/// Try to *refute* a candidate infeasibility by finding a point that satisfies
/// every constraint. Returns `true` when a witness is found — the verdict must
/// then be withdrawn.
///
/// This is the check interval arithmetic cannot do for itself. FBBT and bound
/// propagation reason about *over-approximated ranges*; at extreme coefficient
/// scale those ranges lose the precision the conclusion rests on, and a
/// feasible model can be reported empty. Evaluating the real constraints at a
/// concrete point is exact by comparison — if the point satisfies them, the
/// region demonstrably is not empty, whatever the intervals concluded.
///
/// Measured on a 400-instance feasible-by-construction sweep, Pyomo's
/// `contrib.fbbt` — same tolerance idea, no refutation step — reports 57 false
/// positives where POUNCE reports 3. This gate is where the difference comes
/// from, and is the part worth keeping ahead of it.
///
/// Deliberately one-directional: it can only ever *withdraw* a verdict, never
/// create one, so a failed evaluation or an unsampled witness costs nothing but
/// a missed certification. It is also a single gate covering both certification
/// paths — the previous two safeguards each landed on one path and not its
/// twin, which is how a bound-propagation hole survived.
#[allow(clippy::too_many_arguments)]
fn witness_refutes_infeasibility(
    inner: &Rc<RefCell<dyn TNLP>>,
    n: usize,
    m_in: usize,
    x_l: &[Number],
    x_u: &[Number],
    g_l: &[Number],
    g_u: &[Number],
    tol: Number,
    rule: WitnessRule,
) -> bool {
    if m_in == 0 || n == 0 {
        return false;
    }
    let clamp = |v: Number, fallback: Number| -> Number {
        if v.is_finite() && v.abs() < crate::bound_tighten::INF_BOUND {
            v
        } else {
            fallback
        }
    };
    let lo: Vec<Number> = (0..n).map(|j| clamp(x_l[j], -1.0)).collect();
    let hi: Vec<Number> = (0..n).map(|j| clamp(x_u[j], 1.0)).collect();
    let mid: Vec<Number> = (0..n).map(|j| 0.5 * (lo[j] + hi[j])).collect();

    // The model's own starting point is the first candidate. It is a witness
    // the modeller supplied, and at extreme coefficient scale it is not
    // reproducible from the bounds: `0.5*(lo+hi)` and the modeller's
    // `lo + 0.5*(hi-lo)` differ by an ulp, and against a 1e30 coefficient that
    // ulp moves the row value by ~1e5 — enough to turn a satisfied row into a
    // violated one. Sampling geometry alone therefore missed witnesses that
    // `pounce check-x0` reports as exactly feasible.
    let mut x0 = vec![0.0; n];
    let mut z_l = vec![0.0; n];
    let mut z_u = vec![0.0; n];
    let mut lam = vec![0.0; m_in];
    let have_x0 = inner.borrow_mut().get_starting_point(StartingPoint {
        init_x: true,
        x: &mut x0,
        init_z: false,
        z_l: &mut z_l,
        z_u: &mut z_u,
        init_lambda: false,
        lambda: &mut lam,
    });

    let mut g = vec![0.0; m_in];
    let mut candidates: Vec<&Vec<Number>> = Vec::with_capacity(4);
    if have_x0 && x0.iter().all(|v| v.is_finite()) {
        candidates.push(&x0);
    }
    candidates.extend([&mid, &lo, &hi]);
    for x in candidates {
        if !inner.borrow_mut().eval_g(x, true, &mut g) {
            continue;
        }
        let feasible = (0..m_in).all(|i| {
            let v = g[i];
            if !v.is_finite() {
                return false;
            }
            // An absent bound is stored as the sentinel (`±INF_BOUND`), not as
            // an infinity, so every use of a bound has to ask whether it is
            // real — and the question is **directional**. A lower bound is
            // absent when it is at or below `-INF_BOUND`; an upper bound when
            // it is at or above `+INF_BOUND`. This is the convention
            // [`crate::bound_tighten`] already uses (`x_l[j] <= -INF_BOUND`,
            // `x_u[j] >= INF_BOUND`).
            //
            // gh #396 is what the symmetric magnitude test (`|b| < INF_BOUND`)
            // cost, and it went wrong in both directions at once on the same
            // row. Property-test seed 223 carries
            // `-1e30·x[0] <= -5.0000000000000007e20` over `x[0] ∈ [0, 1e-9]`,
            // with a starting point landing *exactly* on that upper bound:
            //
            //   * `g_l = -1e19` is the sentinel for an absent lower bound, but
            //     `viol` had no presence test at all, so it scored
            //     `g_l - v = 4.9e20` — a violation of a bound the row does not
            //     have, against a point whose true violation is zero.
            //   * `g_u = -5e20` is a genuine, finite upper bound, but
            //     `|−5e20| < 1e19` is false, so the magnitude test called it
            //     absent — discarding a real constraint.
            //
            // Either error alone sinks the witness. With no candidate able to
            // pass, presolve certified a *feasible* model infeasible and
            // reported `solve_result_num` 201, which asserts a proof.
            let lo_present = |b: Number| b.is_finite() && b > -crate::bound_tighten::INF_BOUND;
            let up_present = |b: Number| b.is_finite() && b < crate::bound_tighten::INF_BOUND;
            // Only *real* bounds inform the declared magnitude. The sentinel
            // would otherwise set the slack around 1e11 and make every row look
            // satisfied — which withdrew even correct verdicts.
            let declared = (if lo_present(g_l[i]) {
                g_l[i].abs()
            } else {
                0.0
            })
            .max(if up_present(g_u[i]) {
                g_u[i].abs()
            } else {
                0.0
            });
            let scale = v.abs().max(declared);
            let lo_viol = if lo_present(g_l[i]) { g_l[i] - v } else { 0.0 };
            let hi_viol = if up_present(g_u[i]) { v - g_u[i] } else { 0.0 };
            let viol = lo_viol.max(hi_viol).max(0.0);
            match rule {
                // Accepting direction, so the shared clamped form: never
                // stricter than the absolute `tol`.
                WitnessRule::SolverAcceptance => is_negligible(viol, scale, tol),
                // Pure relative against the *declared* magnitude — no clamp, so
                // the test moves with the row under scaling. A row with no
                // declared magnitude (homogeneous, or unbounded on both sides)
                // has nothing to be relative to and keeps the absolute floor,
                // which is the fail-closed direction here: it accepts the
                // point, and accepting withdraws the verdict.
                //
                // The directional presence tests above make "no declared
                // magnitude" mean what it says. Under the previous symmetric
                // test a row could have a perfectly real bound of `-5e20` and
                // still report `declared == 0`, dropping into the absolute
                // floor by accident rather than by the reasoning above.
                WitnessRule::DeclaredRowRelative if declared > 0.0 => {
                    viol.is_finite() && viol <= tol * declared
                }
                WitnessRule::DeclaredRowRelative => is_negligible(viol, scale, tol),
            }
        });
        if feasible {
            return true;
        }
    }
    false
}

/// Whether the model's numbers are representable enough to base an
/// infeasibility verdict on.
///
/// **Only the FBBT path consults this, deliberately.** Wiring the
/// bound-propagation path to it was tried and measured: every genuinely
/// infeasible model dropped from a certified verdict to the numerical one,
/// because POUNCE stores an unbounded constraint bound as the finite sentinel
/// `INF_BOUND` rather than `f64::INFINITY` — so `is_infinite()` is false, the
/// `< INF_BOUND` test fails on equality, and any model with a one-sided
/// constraint (nearly all of them) is rejected. Making the comparison inclusive
/// admits the sentinel, but then it also admits `1e300`: under POUNCE's own
/// convention `|v| >= INF_BOUND` *means* unbounded, so a bounds-magnitude test
/// cannot separate "no bound" from "absurd finite value". A guard that cannot
/// make that distinction should not be spread to a second path.
///
/// The overflow it does catch on the FBBT path comes from constraint
/// *coefficients*, not bounds; a check on the Jacobian magnitudes would be the
/// mechanism that generalises. See the commit history for the measured
/// alternatives (1e150 -> 99/400, 1e20 -> 37/400, sentinel-inclusive -> 37/400,
/// against a 3/400 baseline). Rejects NaN and finite magnitudes at or beyond the
/// unbounded sentinel, where interval arithmetic overflows; genuine infinities
/// stay admissible because they mean "no bound", not "a number".
fn certificate_data_is_admissible(
    x_l: &[Number],
    x_u: &[Number],
    g_l: &[Number],
    g_u: &[Number],
) -> bool {
    let sane = |v: &Number| {
        let v = *v;
        !v.is_nan() && (v.is_infinite() || v.abs() < crate::bound_tighten::INF_BOUND)
    };
    x_l.iter().all(sane) && x_u.iter().all(sane) && g_l.iter().all(sane) && g_u.iter().all(sane)
}

fn certify_margin(tol: Number) -> Number {
    // Guard against a degenerate or absurdly tight `tol`; below ~1e-12 the
    // margin would drop into the float-noise band the whole check exists to
    // exclude.
    tol.max(1e-12)
}

/// Whether some variable's crossed bounds are crossed by more than the solver
/// itself would accept as feasible — the condition for promoting a Phase-1
/// detection to a certified proof.
///
/// The scale is the crossed pair's own magnitude, through `is_negligible`'s
/// clamped accepting form — deliberately NOT the pure-relative proving form
/// (`is_significant`). The #380 rule binds here: a crossing so small that the
/// solver's acceptance test would wave a point through must never be certified,
/// or the same model reports "proved infeasible" with presolve on and
/// `Solve_Succeeded` with it off. The clamp hazard that rules the proving form
/// out for *row residuals* does not arise for bound crossings: row scaling
/// divides out of a propagated bound (`s*x >= 2*s` implies `x >= 2` at every
/// `s`), so crossings do not shrink in the down-scaled direction. What the
/// relative part fixes is the *up-scaled* model — a crossing of `6e-5` on
/// bounds near `3e11` is float noise (relative `2e-16`), and the absolute
/// margin used to certify it.
fn crossing_is_certifiable(x_l: &[Number], x_u: &[Number], tol: Number) -> bool {
    x_l.iter().zip(x_u).any(|(&l, &u)| {
        // Both sides present before a crossing counts (gh #402) — otherwise an
        // absent lower bound at the `-1e19` sentinel against a real upper bound
        // of `-5e20` reads as a `5e20` crossing, far too large for
        // `is_negligible` to dismiss, and certifies a feasible model infeasible.
        lower_bound_present(l)
            && upper_bound_present(u)
            && l > u
            && !is_negligible(l - u, l.abs().max(u.abs()), tol)
    })
}

/// Re-run FBBT with every constraint's bounds widened by that row's own
/// acceptance slack, `tol * max(|bound|, 1)` — the same clamped form
/// `is_negligible` uses. Returns `true` only if the problem is *still*
/// infeasible with that slack — i.e. the infeasibility is robust at the scale
/// the row is written in, not a hair's-breadth modelling artifact.
///
/// Per-row rather than one absolute margin: an absolute `1e-8` widening is
/// invisible on a row whose bounds sit near `5e13`, so an "infeasibility" of
/// `3e3` there (eleven relative digits — a violation `pounce verify` accepts)
/// would survive the probe and be certified, contradicting the verifier on the
/// same model.
///
/// Sound in the conservative direction: a false `false` merely withholds the
/// proof and leaves the previous behavior (let the IPM decide), while a false
/// `true` would be a wrong "proved infeasible" — so the test is written to fail
/// closed.
/// The acceptance slack for one row: `tol * max(|bound|, 1)` over the row's
/// *present* bounds.
///
/// Only a present bound carries magnitude information, and presence is
/// **directional** (gh #402). The symmetric `|b| < INF_BOUND` test this used to
/// run failed **open**: a row whose real bound is `-5e20` contributed `0.0`, so
/// the margin collapsed to `tol * 1.0` on a row written at `5e20` scale — the
/// infeasibility then survived an absurdly small widening and certified, which
/// is precisely the direction `fbbt_infeasibility_survives_margin` exists to
/// rule out. (The two `relaxed` maps in that function already used the
/// directional form; its halves disagreed.)
fn row_margin_for(g_l: Number, g_u: Number, tol: Number) -> Number {
    let lo_mag = if lower_bound_present(g_l) {
        g_l.abs()
    } else {
        0.0
    };
    let hi_mag = if upper_bound_present(g_u) {
        g_u.abs()
    } else {
        0.0
    };
    tol * lo_mag.max(hi_mag).max(1.0)
}

#[allow(clippy::too_many_arguments)]
fn fbbt_infeasibility_survives_margin(
    provider: &dyn ExpressionProvider,
    n: usize,
    m_in: usize,
    x_l: &[Number],
    x_u: &[Number],
    g_l: &[Number],
    g_u: &[Number],
    row_kept: &[bool],
    cfg: &crate::fbbt::FbbtConfig,
    tol: Number,
) -> bool {
    let row_margin = |i: usize| -> Number { row_margin_for(g_l[i], g_u[i], tol) };
    let g_l_relaxed: Vec<Number> = g_l
        .iter()
        .enumerate()
        .map(|(i, &v)| {
            if v <= -crate::bound_tighten::INF_BOUND {
                v
            } else {
                v - row_margin(i)
            }
        })
        .collect();
    let g_u_relaxed: Vec<Number> = g_u
        .iter()
        .enumerate()
        .map(|(i, &v)| {
            if v >= crate::bound_tighten::INF_BOUND {
                v
            } else {
                v + row_margin(i)
            }
        })
        .collect();
    // Overflow guard. Interval arithmetic on absurd-but-finite magnitudes
    // overflows to +/-inf, and `Interval::is_empty` treats a NaN endpoint as an
    // empty range — so `1e300 * 1e300` reads as "this constraint cannot be
    // satisfied" and would be certified as a *proof*. It is not: the model
    // `x in [-1e300, 1e300], 1e300*x >= 1e300` is plainly feasible (`x = 1`),
    // yet was reported proved infeasible while `presolve_fbbt=no` solved it.
    // Widening the bounds cannot expose this — an overflow is unaffected by a
    // 1e-9 relaxation — so the check has to be on the inputs.
    //
    // A genuine infinity is fine and common (`x >= 5` has `g_u = +inf`): those
    // are the INF sentinel and mean "unbounded". What cannot be trusted is a
    // *finite* magnitude at or beyond the sentinel, which is where the
    // arithmetic stops being representable.
    if !certificate_data_is_admissible(x_l, x_u, g_l, g_u) {
        return false;
    }
    let mut probe_x_l = x_l.to_vec();
    let mut probe_x_u = x_u.to_vec();
    let report = crate::fbbt::run_fbbt(
        provider,
        n,
        m_in,
        &mut probe_x_l,
        &mut probe_x_u,
        &g_l_relaxed,
        &g_u_relaxed,
        Some(row_kept),
        cfg,
    );
    report.infeasibility_witness.is_some()
}

pub struct CachedBounds {
    pub x_l: Vec<Number>,
    pub x_u: Vec<Number>,
    /// Reduced (post-row-drop) constraint lower bounds.
    pub g_l: Vec<Number>,
    /// Reduced constraint upper bounds.
    pub g_u: Vec<Number>,
}

/// TNLP wrapper that re-presents the inner problem after presolve.
pub struct PresolveTnlp {
    inner: Rc<RefCell<dyn TNLP>>,
    /// Optional structural-expression handle on the inner TNLP for
    /// passes (FBBT, issue #62) that need DAG-level access. Callers
    /// who know the concrete inner type (e.g. `pounce-cli` with
    /// `NlTnlp`) install this via
    /// [`Self::with_expression_provider`]. Callers without (e.g.
    /// callback-based bridges) leave it `None` and the expression-
    /// hungry passes silently become no-ops.
    expr_provider: Option<Rc<RefCell<dyn ExpressionProvider>>>,
    opts: PresolveOptions,

    /// Project served primal seeds into the presolved box. Enabled only by
    /// `TnlpPresolveSession`. Direct wrapper users retain non-projecting
    /// behavior.
    project_seed: bool,

    /// Which accepting test the final witness-refutation gate uses. Left at
    /// [`WitnessRule::SolverAcceptance`] for every wrapper that is actually
    /// solved through; raised to [`WitnessRule::DeclaredRowRelative`] only by
    /// callers that have already established the solve cannot run (gh#391).
    /// Deliberately *not* a `PresolveOptions` field: it is a property of the
    /// call site, not a user-tunable knob, and no user should be able to switch
    /// a live solve onto the stricter rule.
    witness_rule: WitnessRule,

    /// `None` until init has run; afterwards `Some(state)`.
    state: Option<PresolveState>,

    /// Full-space `(x, lambda)` forwarded to the inner TNLP at the last
    /// `finalize_solution` — i.e. after row-drop multiplier recovery, in
    /// the original `.nl` row order (length `info_inner.m`). `None` until
    /// a solve finalizes. The CLI prefers this over its reduced-space
    /// `on_converged` / counting capture so the `.sol` / JSON dual block
    /// regains the original constraint count: the kept-row-space lambda
    /// the solver produces otherwise mis-aligns against the `.nl`'s `m`.
    finalized_full_solution: Option<(Vec<Number>, Vec<Number>)>,

    /// `(variable, upper)` for every bound multiplier the last
    /// `finalize_solution` handed back to the dropped singleton row that
    /// implied the bound (M24). A caller whose bound-multiplier capture sits
    /// outside this wrapper zeroes these entries so the `.sol` agrees with
    /// the `lambda` in [`Self::finalized_full_solution`].
    bound_multipliers_handed_to_rows: Vec<(usize, bool)>,

    /// `(variable, z_l, z_u)` for every variable whose bound multipliers the
    /// last `finalize_solution` changed from what the solver handed it — the
    /// singleton hand-back above, Phase 0's fixed-variable zeroing, and the
    /// derived-bound recovery (`dual_recovery`). A caller that captured the
    /// solver's own `z` outside this wrapper overwrites these entries.
    bound_multiplier_overrides: Vec<(usize, Number, Number)>,

    /// `true` when the last `finalize_solution` changed any multiplier
    /// beyond lifting it to the original row space (see
    /// [`Self::postsolve_rewrote_duals`]).
    duals_rewritten: bool,

    /// What the derived-bound recovery did at the last `finalize_solution`.
    derived_bound_recovery: Option<dual_recovery::Outcome>,
    /// The inner evaluations the iteration made, for the live-iterate lift.
    live_record: live_lift::EvalRecord,

    /// The duals the last `get_starting_point` was handed and served, when
    /// it served any: `finalize_solution` returns the handed ones when the
    /// solver gives the served ones back unchanged (`warm_echo`).
    served_warm_start: Option<warm_echo::ServedWarmStart>,

    /// What the most recent starting-point calls actually changed while
    /// projecting the inner seed into this wrapper's row space and box.
    /// Reset by session callers before each solve.
    starting_point_projection_report: WarmProjectionReport,
}

struct PresolveState {
    info_inner: NlpInfo,
    info_outer: NlpInfo,
    bounds: CachedBounds,

    /// Maps outer (reduced) row index → inner row index. Length
    /// equals `info_outer.m`.
    rows_kept: Vec<usize>,

    /// For each outer nnz, the position in the inner nnz array.
    jac_kept_idx: Vec<usize>,
    /// Cached outer (reduced + renumbered) Jacobian sparsity, served
    /// on `eval_jac_g(Structure)`.
    jac_irow_outer: Vec<Index>,
    jac_jcol_outer: Vec<Index>,

    /// Phase 1 report.
    tighten_report: TightenReport,
    /// A proof that the feasible region is empty, when presolve derived one on
    /// an *un-clamped* box. `None` unless one of the four certifiable sites in
    /// `ensure_init` fired — in particular an infeasibility detected while a
    /// Phase-0 auxiliary elimination is in force is deliberately NOT recorded
    /// here, because a dropped row can fabricate one (the #53 / F6 hazard).
    /// Such a detection is re-checked on the rolled-back box first, and only
    /// the surviving result is certified.
    certified_infeasible: Option<InfeasibilityProof>,
    /// FBBT report (`None` when `presolve_fbbt` was off or the inner
    /// TNLP did not expose an `ExpressionProvider`).
    fbbt_report: Option<crate::fbbt::FbbtReport>,
    /// Number of rows dropped by Phase 2.
    n_dropped_rows: Index,
    /// Phase 3 verdict (`None` if the LICQ check was disabled).
    licq_verdict: Option<LicqVerdict>,
    /// Phase 4: warm-start values for `z_l` per variable. Entry is
    /// 0.0 where presolve did not tighten the lower bound, else
    /// `bound_mult_init_val`. Same length as `bounds.x_l`.
    z_l_warm: Vec<Number>,
    /// Phase 4: warm-start values for `z_u` per variable.
    z_u_warm: Vec<Number>,

    /// Scratch buffers reused across eval_* calls.
    scratch_g: Vec<Number>,
    scratch_jac: Vec<Number>,
    scratch_lambda: Vec<Number>,

    /// Phase 0 (issue #53) diagnostics. Always present; defaulted to
    /// zeros when the master switch is off.
    aux_diagnostics: AuxiliaryPreprocessingDiagnostics,
    /// Phase 0 postsolve stack. Empty until PR 7 wires real frames.
    #[allow(dead_code)]
    reduction_stack: ReductionStack,
    /// Bounds Phase 1 created from a singleton linear row that Phase 2 then
    /// dropped (M24): the multiplier the solver reports on such a bound
    /// belongs to the row. See [`SingletonBound`].
    singleton_bounds: Vec<SingletonBound>,
    /// The model's declared box and row bounds, before any presolve pass:
    /// postsolve judges a reported multiplier against these (M24, derived
    /// bounds — see `dual_recovery`).
    decl_x_l: Vec<Number>,
    decl_x_u: Vec<Number>,
    decl_g_l: Vec<Number>,
    decl_g_u: Vec<Number>,
    /// Linear rows Phase 2 dropped as redundant that are not handled as a
    /// [`SingletonBound`], with their coefficients: on a warm start their
    /// supplied multiplier's contribution moves onto the reduced bounds of
    /// the variables they touch (the inverse of the `dual_recovery`
    /// postsolve).
    phase2_dropped_rows: Vec<(usize, Vec<(usize, Number)>)>,
}

/// A variable bound that exists only because a singleton linear row
/// `lo ≤ a·x_j ≤ hi` implied it, and that row was dropped as redundant
/// (issue M24).
///
/// The reduced problem has the bound and not the row, so the solver
/// reports the dual on `z_l[j]` / `z_u[j]`, against a bound the model
/// never declared. Postsolve hands it back to the row: the bound
/// contributes `z_u − z_l` to stationarity at `j` and the row contributes
/// `a·λ` (Ipopt's convention, `∇f + Jᵀλ − z_l + z_u = 0`), so
/// `λ = (z_u − z_l)/a` with that side's `z` zeroed closes stationarity
/// exactly as before. Recorded only where the implied bound is strictly
/// tighter than the declared one; on a tie the multiplier stays on the
/// declared bound (a valid split of a degenerate pair, the same posture
/// `linear_eq_elim` takes for transferred bounds).
#[derive(Debug, Clone, Copy, PartialEq)]
struct SingletonBound {
    /// Inner row index.
    row: usize,
    /// Variable index.
    var: usize,
    /// The row's coefficient on `var`; non-zero.
    coeff: Number,
    /// `true` for the variable's upper bound, `false` for its lower.
    upper: bool,
}

/// Find the bounds Phase 1 took from singleton rows that Phase 2 dropped.
///
/// `kept_before` is the row mask Phase 2 started from (so rows Phase 0
/// consumed, whose multipliers its frames recover, are excluded) and
/// `kept_after` the mask it produced. `x_l` / `x_u` are the presolved box
/// and `decl_l` / `decl_u` the declared one.
fn find_singleton_bounds(
    linear_row_map: &[Option<LinearRow>],
    kept_before: &[bool],
    kept_after: &[bool],
    x_l: &[Number],
    x_u: &[Number],
    decl_l: &[Number],
    decl_u: &[Number],
) -> Vec<SingletonBound> {
    let same = |a: Number, b: Number| (a - b).abs() <= 1e-9 * a.abs().max(b.abs()).max(1.0);
    let mut out = Vec::new();
    // A side already claimed by an earlier row stays with that row.
    let mut claimed: std::collections::BTreeSet<(usize, bool)> = Default::default();
    for (row, lr) in linear_row_map.iter().enumerate() {
        let Some(lr) = lr else { continue };
        if !kept_before.get(row).copied().unwrap_or(false)
            || kept_after.get(row).copied().unwrap_or(true)
        {
            continue;
        }
        let mut nz = lr.entries.iter().filter(|(_, a)| *a != 0.0);
        let (Some(&(j, a)), None) = (nz.next(), nz.next()) else {
            continue;
        };
        let j = j as usize;
        if j >= x_l.len() || !a.is_finite() {
            continue;
        }
        // The variable bounds `lo ≤ a·x ≤ hi` implies, by side.
        let (imp_l, imp_u) = if a > 0.0 {
            (
                lower_bound_present(lr.lo).then(|| lr.lo / a),
                upper_bound_present(lr.hi).then(|| lr.hi / a),
            )
        } else {
            (
                upper_bound_present(lr.hi).then(|| lr.hi / a),
                lower_bound_present(lr.lo).then(|| lr.lo / a),
            )
        };
        if let Some(b) = imp_l
            && lower_bound_present(x_l[j])
            && same(x_l[j], b)
            && (!lower_bound_present(decl_l[j]) || x_l[j] > decl_l[j])
            && claimed.insert((j, false))
        {
            out.push(SingletonBound {
                row,
                var: j,
                coeff: a,
                upper: false,
            });
        }
        if let Some(b) = imp_u
            && upper_bound_present(x_u[j])
            && same(x_u[j], b)
            && (!upper_bound_present(decl_u[j]) || x_u[j] < decl_u[j])
            && claimed.insert((j, true))
        {
            out.push(SingletonBound {
                row,
                var: j,
                coeff: a,
                upper: true,
            });
        }
    }
    out
}

impl PresolveTnlp {
    /// Build a presolve wrapper directly. Prefer
    /// [`wrap_with_presolve`] in production code; this constructor
    /// is exposed so integration tests can keep a typed handle for
    /// accessors like [`Self::licq_verdict`].
    pub fn new(inner: Rc<RefCell<dyn TNLP>>, opts: PresolveOptions) -> Self {
        Self {
            inner,
            expr_provider: None,
            opts,
            project_seed: false,
            witness_rule: WitnessRule::default(),
            state: None,
            finalized_full_solution: None,
            bound_multipliers_handed_to_rows: Vec::new(),
            bound_multiplier_overrides: Vec::new(),
            duals_rewritten: false,
            derived_bound_recovery: None,
            served_warm_start: None,
            live_record: live_lift::EvalRecord::default(),
            starting_point_projection_report: WarmProjectionReport::default(),
        }
    }

    /// Switch the final witness-refutation gate to
    /// [`WitnessRule::DeclaredRowRelative`].
    ///
    /// **Only legitimate when the caller has already established that the solve
    /// cannot run**, so the "never certify what the solver would accept"
    /// counterfactual has nothing to compare against. Today that is the
    /// too-few-degrees-of-freedom gate alone (gh#391). Calling this on a
    /// wrapper that will be solved through reintroduces the #380 defect: a
    /// small-magnitude model would report "proved infeasible" with presolve on
    /// and `Solve_Succeeded` with it off.
    pub fn probing_without_a_solve(mut self) -> Self {
        self.witness_rule = WitnessRule::DeclaredRowRelative;
        self
    }

    /// Enable projection of served warm-start primals into the presolved box.
    pub fn set_project_seed(&mut self, project_seed: bool) {
        self.project_seed = project_seed;
    }

    /// Build a presolve wrapper with an `ExpressionProvider` handle on
    /// the same inner TNLP. The two handles should reference the
    /// *same* object (typical pattern: clone an `Rc<RefCell<NlTnlp>>`
    /// twice, once as `dyn TNLP` and once as `dyn ExpressionProvider`).
    /// Required for `presolve_fbbt=yes` to fire — without a provider,
    /// FBBT silently becomes a no-op.
    pub fn with_expression_provider(
        inner: Rc<RefCell<dyn TNLP>>,
        expr_provider: Rc<RefCell<dyn ExpressionProvider>>,
        opts: PresolveOptions,
    ) -> Self {
        Self {
            inner,
            expr_provider: Some(expr_provider),
            opts,
            project_seed: false,
            witness_rule: WitnessRule::default(),
            state: None,
            finalized_full_solution: None,
            bound_multipliers_handed_to_rows: Vec::new(),
            bound_multiplier_overrides: Vec::new(),
            duals_rewritten: false,
            derived_bound_recovery: None,
            served_warm_start: None,
            live_record: live_lift::EvalRecord::default(),
            starting_point_projection_report: WarmProjectionReport::default(),
        }
    }

    /// FBBT report (`None` until init runs, or when FBBT was disabled
    /// or the inner TNLP did not expose an `ExpressionProvider`).
    pub fn fbbt_report(&self) -> Option<crate::fbbt::FbbtReport> {
        self.state.as_ref().and_then(|s| s.fbbt_report.clone())
    }

    /// Phase 1 report (zeroed until init has run).
    pub fn tighten_report(&self) -> TightenReport {
        self.state
            .as_ref()
            .map(|s| s.tighten_report.clone())
            .unwrap_or_default()
    }

    /// A proof that the feasible region is empty, if presolve derived one on
    /// an un-clamped box. `None` means "not proved", not "feasible".
    ///
    /// Distinct from `tighten_report().infeasible` /
    /// `fbbt_report().infeasibility_witness`, which are raw detections: those
    /// also fire for a contradiction manufactured by a Phase-0 auxiliary
    /// elimination, which is why they were never safe to act on. This returns
    /// only detections that survive on the original box.
    pub fn certified_infeasible(&self) -> Option<InfeasibilityProof> {
        self.state.as_ref().and_then(|s| s.certified_infeasible)
    }

    /// Number of constraint rows dropped by Phase 2 (0 if presolve
    /// has not yet run or no rows are redundant).
    pub fn n_dropped_rows(&self) -> Index {
        self.state.as_ref().map(|s| s.n_dropped_rows).unwrap_or(0)
    }

    /// The full-space `(x, lambda)` captured at the last
    /// `finalize_solution`, lifted back to the original `.nl` row space
    /// (length `info_inner.m`) with dropped-row multipliers recovered.
    /// `None` until a solve finalizes. The CLI consumes this so its
    /// `.sol` / JSON dual block carries the original constraint count
    /// rather than the reduced kept-row count. See the field docs.
    pub fn finalized_full_solution(&self) -> Option<(Vec<Number>, Vec<Number>)> {
        self.finalized_full_solution.clone()
    }

    /// The bound multipliers the last `finalize_solution` moved onto a
    /// dropped singleton row's `lambda`, as `(variable, upper)` (M24). The
    /// inner TNLP already received them zeroed; this is for a caller that
    /// captured the solver's own `z` outside the wrapper.
    pub fn bound_multipliers_handed_to_rows(&self) -> &[(usize, bool)] {
        &self.bound_multipliers_handed_to_rows
    }

    /// `(variable, z_l, z_u)` for every variable whose bound multipliers the
    /// last `finalize_solution` rewrote before forwarding them (M24: the
    /// singleton hand-back, Phase 0's fixed-variable zeroing, and the
    /// derived-bound recovery). The values are what the inner TNLP
    /// received, in its frame (Ipopt's internal convention, both `≥ 0`); a
    /// caller whose bound-multiplier capture sits outside this wrapper
    /// writes them over its own.
    pub fn bound_multiplier_overrides(&self) -> &[(usize, Number, Number)] {
        &self.bound_multiplier_overrides
    }

    /// `true` when the last `finalize_solution` rewrote any multiplier
    /// beyond lifting it into the original row space — so a caller that
    /// holds its own reduced-space capture must take
    /// [`Self::finalized_full_solution`]'s `lambda` even when no row was
    /// dropped.
    pub fn postsolve_rewrote_duals(&self) -> bool {
        self.duals_rewritten
    }

    /// How many multipliers the last `finalize_solution` found on derived
    /// bounds the original model does not have active, and whether it
    /// recovered the original problem's multipliers in their place
    /// (`Some((n_spurious, recovered))`; `None` when there were none). See
    /// the module docs, "Dual attribution".
    pub fn derived_bound_dual_recovery(&self) -> Option<(usize, bool)> {
        self.derived_bound_recovery
            .as_ref()
            .filter(|o| o.n_spurious > 0 || o.rolled_back)
            .map(|o| (o.n_spurious, !o.rolled_back))
    }

    /// Which sides of each variable's box presolve derived: its bound there
    /// is strictly tighter than the declared one, or the declared one is
    /// absent. `None` before init.
    fn derived_sides(&self) -> Option<(Vec<bool>, Vec<bool>)> {
        let s = self.state.as_ref()?;
        let tighter = |p: Number, d: Number, lower: bool| -> bool {
            let present = if lower {
                lower_bound_present(p)
            } else {
                upper_bound_present(p)
            };
            if !present {
                return false;
            }
            let declared = if lower {
                lower_bound_present(d)
            } else {
                upper_bound_present(d)
            };
            if !declared {
                return true;
            }
            let gap = if lower { p - d } else { d - p };
            gap > 1e-9 * d.abs().max(1.0)
        };
        let n = s.bounds.x_l.len();
        let dl = (0..n)
            .map(|j| tighter(s.bounds.x_l[j], s.decl_x_l[j], true))
            .collect();
        let du = (0..n)
            .map(|j| tighter(s.bounds.x_u[j], s.decl_x_u[j], false))
            .collect();
        Some((dl, du))
    }

    /// Evaluate the original problem's `∇f` and Jacobian triplets
    /// (0-based) at `x`.
    fn original_derivatives(
        &mut self,
        x: &[Number],
    ) -> Option<(Vec<Number>, Vec<(usize, usize, Number)>)> {
        let (n, nnz, one_based) = {
            let s = self.state.as_ref()?;
            (
                s.info_inner.n as usize,
                s.info_inner.nnz_jac_g as usize,
                matches!(s.info_inner.index_style, IndexStyle::Fortran),
            )
        };
        let mut grad_f = vec![0.0; n];
        if !self.inner.borrow_mut().eval_grad_f(x, true, &mut grad_f) {
            return None;
        }
        let mut irow = vec![0 as Index; nnz];
        let mut jcol = vec![0 as Index; nnz];
        let mut vals = vec![0.0; nnz];
        if nnz > 0 {
            let mut inner = self.inner.borrow_mut();
            if !inner.eval_jac_g(
                None,
                false,
                SparsityRequest::Structure {
                    irow: &mut irow,
                    jcol: &mut jcol,
                },
            ) || !inner.eval_jac_g(
                Some(x),
                false,
                SparsityRequest::Values { values: &mut vals },
            ) {
                return None;
            }
        }
        let off: Index = if one_based { 1 } else { 0 };
        let jac = (0..nnz)
            .filter_map(|k| {
                let i = irow[k] - off;
                let j = jcol[k] - off;
                (i >= 0 && j >= 0).then_some((i as usize, j as usize, vals[k]))
            })
            .collect();
        Some((grad_f, jac))
    }

    /// The postsolve half of the derived-bound rule (M24; see
    /// `dual_recovery`). Evaluates the original problem only when some
    /// derived side carries a multiplier, so a solve that never leans on a
    /// derived bound is untouched.
    fn recover_derived_bound_duals(
        &mut self,
        x: &[Number],
        g: &[Number],
        lambda: &mut [Number],
        z_l: &mut [Number],
        z_u: &mut [Number],
    ) -> Option<dual_recovery::Outcome> {
        let (derived_l, derived_u) = self.derived_sides()?;
        let n = x.len();
        if z_l.len() != n || z_u.len() != n || derived_l.len() != n {
            return None;
        }
        let touched = (0..n).any(|j| {
            (derived_l[j] && z_l[j].is_finite() && z_l[j] != 0.0)
                || (derived_u[j] && z_u[j].is_finite() && z_u[j] != 0.0)
        });
        if !touched {
            return None;
        }
        let (grad_f, jac) = self.original_derivatives(x)?;
        let s = self.state.as_ref()?;
        let view = dual_recovery::KktView {
            x,
            g,
            grad_f: &grad_f,
            jac: &jac,
            x_l: &s.decl_x_l,
            x_u: &s.decl_x_u,
            g_l: &s.decl_g_l,
            g_u: &s.decl_g_u,
            derived_l: &derived_l,
            derived_u: &derived_u,
        };
        let out = dual_recovery::reattribute(&view, lambda, z_l, z_u);
        if out.rolled_back {
            tracing::warn!(
                target: "pounce_presolve",
                "presolve postsolve: multipliers on derived bounds could not be \
                 re-attributed to the original problem (component past the dense cap \
                 or residual not reduced); the reported duals keep the reduced \
                 problem's attribution"
            );
        } else if out.n_spurious > 0 {
            tracing::debug!(
                target: "pounce_presolve",
                n_spurious = out.n_spurious,
                n_unknowns = out.n_unknowns,
                residual_before = out.residual_before,
                residual_after = out.residual_after,
                "presolve postsolve: re-attributed derived-bound multipliers"
            );
        }
        Some(out)
    }

    /// The original problem's duals when the solver handed back the warm
    /// start this wrapper served without moving its duals: the ones the
    /// caller supplied, `(lambda, z_l, z_u)` (see `warm_echo`). Every
    /// re-derivation postsolve would otherwise do reads `∇f` at the returned
    /// `x`, which the interior-point start pushed off the bounds the supplied
    /// duals belong to. An entry the solver had to fill (a bound multiplier
    /// supplied as 0 or not at all) is reported as filled when the model
    /// declares that bound — what `presolve=no` reports — and as 0 on a
    /// bound only presolve has. `None` on any other return.
    fn echoed_warm_start(
        &self,
        sol: &Solution<'_>,
    ) -> Option<(Vec<Number>, Vec<Number>, Vec<Number>)> {
        let w = self.served_warm_start.as_ref()?;
        if !w.echoed_by(sol.z_l, sol.z_u, sol.lambda) {
            return None;
        }
        let s = self.state.as_ref()?;
        let n = w.z_l_in.len();
        if w.z_u_in.len() != n || s.decl_x_l.len() != n || s.bounds.x_l.len() != n {
            return None;
        }
        let side = |z_in: &[Number], z_out: &[Number], got: &[Number], lower: bool| {
            (0..n)
                .map(|j| {
                    let mut z = warm_echo::as_reported_z(z_in[j]);
                    let declared = if lower {
                        lower_bound_present(s.decl_x_l[j])
                    } else {
                        upper_bound_present(s.decl_x_u[j])
                    };
                    // A variable presolve fixed is not one the solver saw:
                    // whatever it reports there is not a filled seed.
                    let pinned =
                        s.bounds.x_l[j] == s.bounds.x_u[j] && s.decl_x_l[j] != s.decl_x_u[j];
                    if !warm_echo::is_seed(z_out[j]) && declared && !pinned && got[j].is_finite() {
                        z += got[j];
                    }
                    z
                })
                .collect::<Vec<_>>()
        };
        Some((
            w.lambda_in
                .iter()
                .map(|&l| warm_echo::as_reported_lambda(l))
                .collect(),
            side(&w.z_l_in, &w.z_l_out, sol.z_l, true),
            side(&w.z_u_in, &w.z_u_out, sol.z_u, false),
        ))
    }

    /// Record what `finalize_solution` is about to hand the inner TNLP
    /// (overrides against the solver's `z`, whether any dual was rewritten,
    /// the full-space `(x, λ)`), then hand it down.
    #[allow(clippy::too_many_arguments)]
    fn forward_finalized(
        &mut self,
        sol: Solution<'_>,
        g_full: &[Number],
        lambda_full: &[Number],
        lambda_lifted: &[Number],
        z_l_full: &[Number],
        z_u_full: &[Number],
        ip_data: &IpoptData,
        ip_cq: &IpoptCq,
    ) {
        self.bound_multiplier_overrides = (0..z_l_full.len().min(sol.z_l.len()))
            .filter(|&j| {
                z_l_full[j].to_bits() != sol.z_l[j].to_bits()
                    || z_u_full.get(j).map(|z| z.to_bits()) != sol.z_u.get(j).map(|z| z.to_bits())
            })
            .map(|j| (j, z_l_full[j], z_u_full.get(j).copied().unwrap_or(0.0)))
            .collect();
        self.duals_rewritten = !self.bound_multiplier_overrides.is_empty()
            || lambda_full
                .iter()
                .zip(lambda_lifted)
                .any(|(a, b)| a.to_bits() != b.to_bits());
        // Stash the full-space (x, lambda) the inner TNLP is about to
        // receive so the CLI can prefer it over its reduced-space
        // capture when sizing the `.sol` / JSON dual block.
        self.finalized_full_solution = Some((sol.x.to_vec(), lambda_full.to_vec()));
        self.inner.borrow_mut().finalize_solution(
            Solution {
                status: sol.status,
                x: sol.x,
                z_l: z_l_full,
                z_u: z_u_full,
                g: g_full,
                lambda: lambda_full,
                obj_value: sol.obj_value,
            },
            ip_data,
            ip_cq,
        );
    }

    /// Clear the record of transformations applied by
    /// [`TNLP::get_starting_point`]. A persistent caller resets this before
    /// each solve, then reads [`Self::starting_point_projection_report`].
    pub fn reset_starting_point_projection_report(&mut self) {
        self.starting_point_projection_report = WarmProjectionReport::default();
    }

    /// Projection effects observed while this wrapper actually served the
    /// current solve's starting point.
    pub fn starting_point_projection_report(&self) -> WarmProjectionReport {
        self.starting_point_projection_report
    }

    /// Cached reduced bounds, if presolve has run.
    pub fn cached_bounds(&self) -> Option<&CachedBounds> {
        self.state.as_ref().map(|s| &s.bounds)
    }

    /// Phase 3 verdict — `Some` only if the LICQ check was enabled
    /// and presolve has run.
    pub fn licq_verdict(&self) -> Option<&LicqVerdict> {
        self.state.as_ref().and_then(|s| s.licq_verdict.as_ref())
    }

    /// Phase 4 warm-start hints `(z_l, z_u)`. Each entry is 0.0 if
    /// no hint is set for that variable, else the configured
    /// `bound_mult_init_val`. `None` until init has run.
    pub fn z_warm_starts(&self) -> Option<(&[Number], &[Number])> {
        self.state
            .as_ref()
            .map(|s| (&s.z_l_warm[..], &s.z_u_warm[..]))
    }

    /// Phase 0 (issue #53) summary. Returns a zero-valued struct
    /// until init has run; afterwards, populated by
    /// [`auxiliary::run_auxiliary_phase0`]. PR 1 always returns
    /// zeros even with the master switch on.
    pub fn auxiliary_diagnostics(&self) -> AuxiliaryPreprocessingDiagnostics {
        self.state
            .as_ref()
            .map(|s| s.aux_diagnostics.clone())
            .unwrap_or_default()
    }

    /// Snapshot the transformation for warm-start mapping. Runs init first.
    /// `None` when init fails. See [`warm`] for the contract.
    pub fn transformation(&mut self) -> Option<PresolveMap> {
        let s = self.ensure_init()?;
        let mut fixed: std::collections::BTreeMap<usize, Number> =
            std::collections::BTreeMap::new();
        for frame in s.reduction_stack.iter_bottom_up() {
            for (k, &i) in frame.fixed_vars.iter().enumerate() {
                if let Some(&v) = frame.fixed_values.get(k) {
                    fixed.insert(i, v);
                }
            }
        }
        let (fixed_vars, fixed_values): (Vec<usize>, Vec<Number>) = fixed.into_iter().unzip();
        Some(PresolveMap {
            n_inner: s.info_inner.n.max(0) as usize,
            m_inner: s.info_inner.m.max(0) as usize,
            m_outer: s.info_outer.m.max(0) as usize,
            rows_kept: s.rows_kept.clone(),
            x_l: s.bounds.x_l.clone(),
            x_u: s.bounds.x_u.clone(),
            fixed_vars,
            fixed_values,
        })
    }

    /// Fingerprint what this transformation was computed from. A session
    /// rebuilds the wrapper when this moves. `None` means "always rebuild".
    pub fn fingerprint(&mut self) -> Option<PresolveFingerprint> {
        crate::warm::compute_fingerprint(&self.inner, &self.opts)
    }

    /// Drop the cached transformation, the next query recomputes it. Needed
    /// after fingerprint-invisible changes (FBBT tapes).
    pub fn invalidate(&mut self) {
        self.state = None;
        self.served_warm_start = None;
        self.reset_starting_point_projection_report();
    }

    /// Lazy initialization: pull inner dims, bounds, linearity tags,
    /// Jacobian, and starting point; run Phase 1 + Phase 2 passes;
    /// cache everything needed to translate later eval_* calls.
    fn ensure_init(&mut self) -> Option<&PresolveState> {
        if self.state.is_some() {
            return self.state.as_ref();
        }

        let info_inner = self.inner.borrow_mut().get_nlp_info()?;
        let n = info_inner.n as usize;
        let m_in = info_inner.m as usize;
        let nnz_in = info_inner.nnz_jac_g as usize;

        // Inner bounds.
        let mut x_l = vec![0.0; n];
        let mut x_u = vec![0.0; n];
        let mut g_l_inner = vec![0.0; m_in];
        let mut g_u_inner = vec![0.0; m_in];
        {
            let mut inner = self.inner.borrow_mut();
            if !inner.get_bounds_info(BoundsInfo {
                x_l: &mut x_l,
                x_u: &mut x_u,
                g_l: &mut g_l_inner,
                g_u: &mut g_u_inner,
            }) {
                return None;
            }
        }

        // Jacobian sparsity.
        let mut jac_irow_inner = vec![0 as Index; nnz_in];
        let mut jac_jcol_inner = vec![0 as Index; nnz_in];
        if nnz_in > 0 {
            let mut inner = self.inner.borrow_mut();
            if !inner.eval_jac_g(
                None,
                false,
                SparsityRequest::Structure {
                    irow: &mut jac_irow_inner,
                    jcol: &mut jac_jcol_inner,
                },
            ) {
                return None;
            }
        }

        // Linearity tags (presolve is dormant without them).
        let mut linearity = vec![Linearity::NonLinear; m_in];
        let have_linearity = if m_in > 0 {
            self.inner
                .borrow_mut()
                .get_constraints_linearity(&mut linearity)
        } else {
            true
        };

        // Per-variable linearity (H11): lets Phase-0 objective-coupling
        // classification distinguish a genuinely objective-free variable from
        // one that is merely zero-gradient at the single probe point. The
        // objective-scoped query is preferred — it is exactly the set the
        // guard needs. The global tags are only a fallback: they are a
        // conservative superset (a variable nonlinear only in a *constraint*
        // is tagged NonLinear too), which keeps the guard sound but blocks
        // legitimate eliminations of objective-free blocks (the gas-network
        // case). Most TNLPs decline both (default `false`), in which case
        // Phase 0 falls back to the probe gradient alone.
        let mut var_linearity = vec![Linearity::NonLinear; n];
        let have_var_linearity = {
            let mut inner = self.inner.borrow_mut();
            inner.get_objective_variables_linearity(&mut var_linearity)
                || inner.get_variables_linearity(&mut var_linearity)
        };

        // Probe point for Jacobian values (linear rows have constant
        // Jacobians; this `x` is only needed because some inner
        // TNLPs assert on receipt).
        let mut x_probe = vec![0.0; n];
        let mut z_l_probe = vec![0.0; n];
        let mut z_u_probe = vec![0.0; n];
        let mut lambda_probe = vec![0.0; m_in];
        let started = self.inner.borrow_mut().get_starting_point(StartingPoint {
            init_x: true,
            x: &mut x_probe,
            init_z: false,
            z_l: &mut z_l_probe,
            z_u: &mut z_u_probe,
            init_lambda: false,
            lambda: &mut lambda_probe,
        });
        if !started {
            return None;
        }

        // Jacobian values at the probe.
        let mut jac_values_inner = vec![0.0; nnz_in];
        if nnz_in > 0 {
            let ok = self.inner.borrow_mut().eval_jac_g(
                Some(&x_probe),
                true,
                SparsityRequest::Values {
                    values: &mut jac_values_inner,
                },
            );
            if !ok {
                return None;
            }
        }

        // Build LinearRow list from the inner Jacobian + linearity.
        let one_based = matches!(info_inner.index_style, IndexStyle::Fortran);
        let mut by_row: Vec<Vec<(Index, Number)>> = vec![Vec::new(); m_in];
        for k in 0..nnz_in {
            let i = if one_based {
                (jac_irow_inner[k] - 1) as usize
            } else {
                jac_irow_inner[k] as usize
            };
            let j = if one_based {
                jac_jcol_inner[k] - 1
            } else {
                jac_jcol_inner[k]
            };
            if i < m_in && (j as usize) < n {
                by_row[i].push((j, jac_values_inner[k]));
            }
        }
        let linear_row_map: Vec<Option<LinearRow>> = (0..m_in)
            .map(|i| {
                if have_linearity && matches!(linearity[i], Linearity::Linear) {
                    Some(LinearRow {
                        entries: by_row[i].clone(),
                        lo: g_l_inner[i],
                        hi: g_u_inner[i],
                    })
                } else {
                    None
                }
            })
            .collect();
        // NOTE: `linear_rows` is materialised AFTER Phase 0 so it
        // can filter out rows Phase 0 dropped — propagating bounds
        // through aux-eliminated rows lets tighten_bounds derive
        // contradictions (see issue #53 PR review).

        // Snapshot inner bounds before Phase 1 mutates them; needed
        // for Phase 4 warm-start hints AND for rolling back Phase 0
        // if its clamps later prove infeasible against the kept
        // linear rows.
        let inner_x_l = x_l.clone();
        let inner_x_u = x_u.clone();

        // Phase 0 (issue #53): auxiliary-equality preprocessing.
        // PR 8 wires the real pipeline (incidence → matching → DM →
        // BTF → classify → linear block solve → frame). Variables it
        // fixes are clamped via x_l/x_u; rows it drops are recorded
        // in `row_kept_inner` so the existing remapping logic below
        // picks them up.
        let mut row_kept_inner: Vec<bool> = vec![true; m_in];
        let mut reduction_stack = ReductionStack::default();
        let aux_diagnostics = if self.opts.auxiliary && m_in > 0 {
            // Probe extra quantities Phase 0 needs.
            let mut g_at_probe = vec![0.0; m_in];
            let g_ok = self
                .inner
                .borrow_mut()
                .eval_g(&x_probe, true, &mut g_at_probe);
            if !g_ok {
                return None;
            }
            let mut grad_f_probe = vec![0.0; n];
            let grad_ok = self
                .inner
                .borrow_mut()
                .eval_grad_f(&x_probe, false, &mut grad_f_probe);
            if !grad_ok {
                return None;
            }
            // Plug everything into the orchestrator.
            let linearity_for_phase0: Vec<Linearity> = if have_linearity {
                linearity.clone()
            } else {
                vec![Linearity::NonLinear; m_in]
            };
            let probe_view = auxiliary::Phase0Probe {
                n_vars: n,
                n_rows: m_in,
                jac_irow: &jac_irow_inner,
                jac_jcol: &jac_jcol_inner,
                jac_values: &jac_values_inner,
                g_l: &g_l_inner,
                g_u: &g_u_inner,
                g_at_probe: &g_at_probe,
                linearity: &linearity_for_phase0,
                one_based,
                eq_tol: 1e-12,
                x_probe: &x_probe,
                grad_f: &grad_f_probe,
                var_linearity: if have_var_linearity {
                    Some(&var_linearity)
                } else {
                    None
                },
                x_l: &x_l,
                x_u: &x_u,
            };
            // Adapter: wrap `self.inner` so the orchestrator can
            // call eval_g / eval_jac_g for nonlinear blocks.
            struct TnlpCallbackAdapter {
                inner: Rc<RefCell<dyn TNLP>>,
            }
            impl auxiliary::Phase0TnlpCallback for TnlpCallbackAdapter {
                fn eval_g_full(&mut self, x: &[Number], g: &mut [Number]) -> bool {
                    self.inner.borrow_mut().eval_g(x, true, g)
                }
                fn eval_jac_g_values(&mut self, x: &[Number], values: &mut [Number]) -> bool {
                    self.inner.borrow_mut().eval_jac_g(
                        Some(x),
                        true,
                        SparsityRequest::Values { values },
                    )
                }
            }
            let mut adapter = TnlpCallbackAdapter {
                inner: Rc::clone(&self.inner),
            };
            let mut large_solver = block_solve::RelaxedNewtonSolver;
            let plan = auxiliary::run_auxiliary_phase0(
                &self.opts,
                &probe_view,
                Some(&mut adapter),
                Some(&mut large_solver),
            );
            if let Some(frame) = plan.frame {
                // Clamp fixed variables.
                for (k, &i) in frame.fixed_vars.iter().enumerate() {
                    x_l[i] = frame.fixed_values[k];
                    x_u[i] = frame.fixed_values[k];
                }
                // Drop dropped rows.
                for &r in &frame.dropped_rows {
                    row_kept_inner[r] = false;
                }
                reduction_stack.push(frame);
            }
            // When `presolve_auxiliary_diagnostics=yes`, emit the
            // summary through tracing (pounce#71).
            if self.opts.auxiliary_diagnostics {
                tracing::info!(target: "pounce::presolve", "{}", plan.diagnostics);
            }
            plan.diagnostics
        } else {
            AuxiliaryPreprocessingDiagnostics::default()
        };

        // Build `linear_rows` excluding rows Phase 0 dropped. This
        // is the headline fix from the #53 PR review: propagating
        // bounds through an aux-dropped row lets tighten_bounds
        // derive `x_l[j] > x_u[j]` for an aux-clamped variable and
        // then hand corrupted bounds to the IPM.
        let mut linear_rows: Vec<LinearRow> = linear_row_map
            .iter()
            .enumerate()
            .filter_map(|(i, r)| if row_kept_inner[i] { r.clone() } else { None })
            .collect();

        // Phase 1: bound tightening using linear rows.
        let mut tighten_report = TightenReport::default();
        // Set only where a contradiction is derived on an un-clamped box; see
        // the field doc on `PresolveState::certified_infeasible`.
        let mut certified_infeasible: Option<InfeasibilityProof> = None;
        if self.opts.bound_tightening && !linear_rows.is_empty() {
            tighten_report = tighten_bounds(
                &linear_rows,
                &mut x_l,
                &mut x_u,
                self.opts.max_passes,
                1e-12,
            );
        }

        // Defence in depth: if Phase 1 still flags infeasibility AND
        // Phase 0 made changes, those changes are presumed to blame
        // (aux solved a block to a point inconsistent with bounds
        // from kept rows). Roll back Phase 0 — restore bounds, undo
        // row drops, clear the reduction stack — and re-run Phase 1
        // on the un-filtered linear rows. Without this guard,
        // `report.infeasible` was previously never inspected and
        // corrupted bounds reached the IPM (#53 PR review).
        if tighten_report.infeasible && !reduction_stack.is_empty() {
            tracing::warn!(
                target: "pounce::presolve",
                "auxiliary-equality elimination produced bounds inconsistent \
                 with kept linear rows; rolling back the elimination for this solve."
            );
            x_l.copy_from_slice(&inner_x_l);
            x_u.copy_from_slice(&inner_x_u);
            for kept in row_kept_inner.iter_mut() {
                *kept = true;
            }
            reduction_stack = ReductionStack::default();
            // Re-run tighten on the unfiltered linear rows now that
            // the aux clamps are gone. Rebuild `linear_rows` to the
            // full set too, so Phase 2's redundancy mask stays aligned
            // with the now all-kept `row_kept_inner` (C1).
            let full_linear_rows: Vec<LinearRow> =
                linear_row_map.iter().filter_map(|r| r.clone()).collect();
            tighten_report = TightenReport::default();
            if self.opts.bound_tightening && !full_linear_rows.is_empty() {
                tighten_report = tighten_bounds(
                    &full_linear_rows,
                    &mut x_l,
                    &mut x_u,
                    self.opts.max_passes,
                    1e-12,
                );
            }
            linear_rows = full_linear_rows;
        }

        // M25: a *genuine* Phase-1 infeasibility (empty feasible box) must not
        // reach the IPM as crossed bounds `x_l > x_u`. `tighten_pass` returns
        // the moment it detects `x_l[j] > x_u[j]`, leaving those crossed
        // bounds in place; the rollback above only fires when Phase 0 made
        // changes (`!reduction_stack.is_empty()`), so an infeasibility found
        // with an empty reduction stack — or one that survives the rollback
        // re-tighten — would otherwise hand a degenerate box to the solver,
        // which reports an invalid-problem error instead of a clean
        // infeasibility verdict. Presolve has no channel to certify
        // infeasibility, so — mirroring the aux rollback above and the FBBT
        // handling below — restore the original inner box (always a valid box)
        // and let the IPM run on it and certify infeasibility itself. The
        // `tighten_report.infeasible` flag is preserved and surfaced via
        // `tighten_report()` for diagnostics.
        if tighten_report.infeasible {
            // Measure the crossing *before* the restore below overwrites it.
            let crossing = (0..n).map(|j| x_l[j] - x_u[j]).fold(0.0_f64, f64::max);
            // Reaching here implies the tightening ran on an un-clamped box:
            // the aux rollback above either did not fire (empty reduction
            // stack ⇒ Phase 0 changed nothing) or fired and recomputed
            // `tighten_report` on the restored box. Either way this is a
            // genuine contradiction in the user's model, not a presolve
            // artifact, so it is certifiable. The crossed bounds are still
            // discarded — a degenerate box must not reach the solver — but the
            // *verdict* now survives.
            // Phase 1 flags a crossing above its own `1e-12` test, which is
            // far below what the solver treats as feasible. Certify only a
            // crossing the solver's own acceptance would call a real violation
            // at that variable's scale (see `crossing_is_certifiable`), so a
            // model POUNCE would solve to `Solve_Succeeded` can never be
            // called proved infeasible.
            let robust = crossing_is_certifiable(&x_l, &x_u, certify_margin(self.opts.certify_tol));
            x_l.copy_from_slice(&inner_x_l);
            x_u.copy_from_slice(&inner_x_u);
            if robust {
                certified_infeasible = Some(InfeasibilityProof::BoundPropagation);
            }
            tracing::warn!(
                target: "pounce::presolve",
                crossing,
                certified = robust,
                "Phase 1 bound tightening found the feasible region empty; its \
                 crossed bounds are being discarded."
            );
        }

        // Phase 1 declares infeasibility only when the crossing exceeds its
        // `1e-12` margin, so a *sub-margin* crossing survives the block above:
        // `infeasible` is false, nothing is restored, and `x_l[j] > x_u[j]`
        // by a hair is handed to the solver, which rejects it as
        // `Invalid_Problem_Definition`. The model that surfaced this is
        // `x >= 0.1 + 0.2` with `x <= 0.3` — in binary floating point the two
        // differ by `5.5e-17`, so the box arrives crossed by that much and a
        // model POUNCE otherwise solves cleanly (`presolve=no` gives
        // `Solve_Succeeded`, the LP route gives "Optimal Solution Found")
        // failed with an invalid-problem error the moment presolve was on.
        //
        // Collapse those to a point. The crossing is below the margin at which
        // the tightening itself is willing to call the region empty, so the
        // honest reading is a single feasible value, not an empty set — and it
        // keeps every route's answer consistent.
        // The collapsed point must stay inside the *declared* box. The crossing
        // exists because tightening moved a bound past its partner, so the
        // midpoint of the crossed pair lies strictly outside the original range
        // — for `x >= 0.1 + 0.2`, `x <= 0.3` it lands on `0.30000000000000004`,
        // above the user's `x <= 0.3`. Returning that would silently violate a
        // declared bound, which is worse than the error it replaces and defeats
        // `honor_original_bounds`. Clamping to the original box keeps the
        // reported point admissible.
        //
        // Only *negligible* crossings collapse — `is_negligible` at the
        // crossed pair's own magnitude, the exact complement of
        // `crossing_is_certifiable`, so every crossing is either collapsed
        // here or certifiable there, never both. A crossing the solver itself
        // would call a real violation is not float noise; it is a
        // user-declared empty box (`x in [5, 3]`) on a variable no linear row
        // ever propagated into, which Phase 1 therefore never flagged.
        // Collapsing that would turn `Invalid_Problem_Definition` into
        // `Solve_Succeeded` at a point the model excludes — a wrong answer
        // replacing a correct error. Those stay crossed so the solver rejects
        // them exactly as it does with presolve off.
        let collapse_tol = certify_margin(self.opts.certify_tol);
        for j in 0..n {
            if x_l[j] > x_u[j]
                && is_negligible(
                    x_l[j] - x_u[j],
                    x_l[j].abs().max(x_u[j].abs()),
                    collapse_tol,
                )
            {
                let mid = (0.5 * (x_l[j] + x_u[j]))
                    .max(inner_x_l[j])
                    .min(inner_x_u[j]);
                x_l[j] = mid;
                x_u[j] = mid;
            }
        }

        // Phase 1b — FBBT (issue #62). Runs interval arithmetic over
        // each nonlinear constraint's expression DAG to tighten
        // variable bounds further. No-op when (a) `presolve_fbbt` is
        // off, (b) the inner TNLP did not supply an
        // `ExpressionProvider`, or (c) the problem has zero
        // constraints. Honors `fbbt_tol`, `fbbt_max_iter`, and
        // `fbbt_max_constraints`.
        let mut fbbt_report: Option<crate::fbbt::FbbtReport> = None;
        if self.opts.fbbt && m_in > 0 {
            if let Some(provider) = self.expr_provider.as_ref() {
                let cfg = crate::fbbt::FbbtConfig {
                    tol: self.opts.fbbt_tol,
                    max_iter: self.opts.fbbt_max_iter.max(1) as usize,
                    max_constraints: self.opts.fbbt_max_constraints.max(0) as usize,
                };
                // H12: snapshot the bounds before FBBT. The `FbbtReport`
                // contract states that on detected infeasibility the
                // variable bounds are "undefined and must not be trusted"
                // — FBBT may have partially tightened several variables
                // before a later constraint proved the box empty. Pass the
                // Phase-0 `row_kept_inner` mask so propagation skips any row
                // an auxiliary elimination dropped (over the aux-clamped box
                // a dropped row can fabricate a spurious infeasibility — the
                // same #53 hazard Phase 1 guards against with filtered rows).
                let fbbt_x_l_pre = x_l.clone();
                let fbbt_x_u_pre = x_u.clone();
                let provider_borrow = provider.borrow();
                let mut report = crate::fbbt::run_fbbt(
                    &*provider_borrow,
                    n,
                    m_in,
                    &mut x_l,
                    &mut x_u,
                    &g_l_inner,
                    &g_u_inner,
                    Some(&row_kept_inner),
                    &cfg,
                );
                drop(provider_borrow);
                if report.infeasibility_witness.is_some() && !reduction_stack.is_empty() {
                    // F6: an FBBT infeasibility found while a Phase-0
                    // auxiliary elimination is in force may be an *artifact*
                    // of that elimination — an aux clamp can make a kept
                    // *nonlinear* row infeasible over the clamped box even
                    // though the original problem is feasible. Restoring only
                    // the pre-FBBT box (the `else` branch below) leaves the
                    // aux clamps in place, so the IPM would then cleanly
                    // certify infeasibility of a problem presolve itself broke
                    // — a wrong "infeasible" verdict on a feasible original.
                    // Mirror the Phase-1 aux rollback (the `tighten_report`
                    // branch above): undo Phase 0 entirely — restore the inner
                    // box, re-keep every dropped row, clear the reduction
                    // stack, rebuild the full linear-row set — then re-run
                    // Phase 1 and FBBT on the un-clamped box. Only an
                    // infeasibility that *survives* on the un-clamped box is
                    // genuine; presolve still cannot certify it, so it then
                    // falls through to the same "discard FBBT bounds, let the
                    // IPM certify" handling.
                    tracing::warn!(
                        target: "pounce::presolve",
                        witness = report.infeasibility_witness,
                        "FBBT reported a constraint infeasibility while auxiliary \
                         elimination was active; rolling back the elimination and \
                         re-running FBBT on the un-clamped box to avoid certifying \
                         a presolve-induced infeasibility on a feasible original."
                    );
                    x_l.copy_from_slice(&inner_x_l);
                    x_u.copy_from_slice(&inner_x_u);
                    for kept in row_kept_inner.iter_mut() {
                        *kept = true;
                    }
                    reduction_stack = ReductionStack::default();
                    // Rebuild `linear_rows` to the full set so Phase 2's
                    // redundancy mask stays aligned with the now all-kept
                    // `row_kept_inner` (C1), mirroring the Phase-1 rollback.
                    let full_linear_rows: Vec<LinearRow> =
                        linear_row_map.iter().filter_map(|r| r.clone()).collect();
                    tighten_report = TightenReport::default();
                    if self.opts.bound_tightening && !full_linear_rows.is_empty() {
                        tighten_report = tighten_bounds(
                            &full_linear_rows,
                            &mut x_l,
                            &mut x_u,
                            self.opts.max_passes,
                            1e-12,
                        );
                    }
                    linear_rows = full_linear_rows;
                    // M25: a genuine Phase-1 infeasibility on the un-clamped
                    // box must not reach the IPM as crossed bounds.
                    if tighten_report.infeasible {
                        // Measure the crossing before the restore wipes it.
                        let crossing = (0..n).map(|j| x_l[j] - x_u[j]).fold(0.0_f64, f64::max);
                        // Derived on the rolled-back, un-clamped box, so the
                        // aux elimination cannot be to blame — certifiable when
                        // the crossing is real at the variable's own scale.
                        let robust = crossing_is_certifiable(
                            &x_l,
                            &x_u,
                            certify_margin(self.opts.certify_tol),
                        );
                        x_l.copy_from_slice(&inner_x_l);
                        x_u.copy_from_slice(&inner_x_u);
                        if robust {
                            certified_infeasible = Some(InfeasibilityProof::BoundPropagation);
                        }
                        tracing::warn!(
                            target: "pounce::presolve",
                            crossing,
                            certified = robust,
                            "Phase 1 bound tightening on the rolled-back (un-clamped) \
                             box found the feasible region empty; its crossed bounds \
                             are being discarded."
                        );
                    }
                    // Re-run FBBT on the un-clamped, all-kept box.
                    let rerun_x_l_pre = x_l.clone();
                    let rerun_x_u_pre = x_u.clone();
                    let provider_borrow = provider.borrow();
                    report = crate::fbbt::run_fbbt(
                        &*provider_borrow,
                        n,
                        m_in,
                        &mut x_l,
                        &mut x_u,
                        &g_l_inner,
                        &g_u_inner,
                        Some(&row_kept_inner),
                        &cfg,
                    );
                    drop(provider_borrow);
                    if let Some(witness) = report.infeasibility_witness {
                        // Survives without the aux clamp: a genuine
                        // infeasibility of the original problem, so it is
                        // certifiable *if* it is more than a hair's breadth —
                        // see `certify_margin`. FBBT's own bounds are
                        // still undefined per the `FbbtReport` contract and are
                        // discarded either way.
                        let provider_borrow = provider.borrow();
                        let robust = fbbt_infeasibility_survives_margin(
                            &*provider_borrow,
                            n,
                            m_in,
                            &inner_x_l,
                            &inner_x_u,
                            &g_l_inner,
                            &g_u_inner,
                            &row_kept_inner,
                            &cfg,
                            certify_margin(self.opts.certify_tol),
                        );
                        drop(provider_borrow);
                        if robust {
                            certified_infeasible =
                                Some(InfeasibilityProof::IntervalArithmetic { witness });
                        }
                        tracing::warn!(
                            target: "pounce::presolve",
                            witness,
                            certified = robust,
                            "FBBT still reports infeasibility on the un-clamped box; \
                             treating it as genuine."
                        );
                        x_l.copy_from_slice(&rerun_x_l_pre);
                        x_u.copy_from_slice(&rerun_x_u_pre);
                    }
                } else if let Some(witness) = report.infeasibility_witness {
                    // No aux elimination active (empty reduction stack): the
                    // witnessed infeasibility cannot be a Phase-0 artifact, so
                    // it is certifiable *if* it is more than a hair's breadth —
                    // see `certify_margin`. FBBT's tightened bounds are
                    // undefined per the `FbbtReport` contract and are discarded
                    // either way.
                    let provider_borrow = provider.borrow();
                    let robust = fbbt_infeasibility_survives_margin(
                        &*provider_borrow,
                        n,
                        m_in,
                        &inner_x_l,
                        &inner_x_u,
                        &g_l_inner,
                        &g_u_inner,
                        &row_kept_inner,
                        &cfg,
                        certify_margin(self.opts.certify_tol),
                    );
                    drop(provider_borrow);
                    if robust {
                        certified_infeasible =
                            Some(InfeasibilityProof::IntervalArithmetic { witness });
                    }
                    tracing::warn!(
                        target: "pounce::presolve",
                        witness,
                        certified = robust,
                        "FBBT reported a constraint's feasible range empty; its \
                         tightened bounds are undefined and are being discarded."
                    );
                    x_l.copy_from_slice(&fbbt_x_l_pre);
                    x_u.copy_from_slice(&fbbt_x_u_pre);
                }
                fbbt_report = Some(report);
            }
        }

        // Final admissibility gate for the infeasibility verdict, covering
        // BOTH certification paths in one place.
        //
        // Interval reasoning over-approximates; at extreme coefficient scale it
        // can report a feasible region empty. Before the verdict escapes, try
        // to refute it by evaluating the real constraints at concrete points in
        // the declared box. A witness means the region is demonstrably not
        // empty, whatever the intervals concluded.
        //
        // Placed here, once, rather than beside each detection: the two
        // previous safeguards were each added to one path and not its twin,
        // and a hole survived both times.
        //
        // Which accepting test decides "satisfies" is the call site's choice,
        // not an option — see `WitnessRule`. Every wrapper that is solved
        // through uses the solver's own acceptance; only a caller that has
        // established the solve cannot run may raise it to the declared,
        // scale-relative form.
        if certified_infeasible.is_some()
            && witness_refutes_infeasibility(
                &self.inner,
                n,
                m_in,
                &inner_x_l,
                &inner_x_u,
                &g_l_inner,
                &g_u_inner,
                self.opts.certify_tol,
                self.witness_rule,
            )
        {
            tracing::warn!(
                target: "pounce::presolve",
                "presolve flagged the feasible region empty, but a point in the \
                 declared box satisfies every constraint — withdrawing the \
                 verdict and letting the solver decide."
            );
            certified_infeasible = None;
        }

        // Phase 4: any variable whose lower (upper) bound moved
        // strictly inward is a candidate for a bound-multiplier warm
        // start. Zero entries leave that bound's multiplier on the
        // global default (`bound_mult_init_val` from upstream).
        let warm_tol: Number = 1e-12;
        let (z_l_warm, z_u_warm) = if self.opts.warm_z_bounds {
            let v0 = self.opts.bound_mult_init_val;
            let mut zl = vec![0.0; n];
            let mut zu = vec![0.0; n];
            for i in 0..n {
                if x_l[i] > inner_x_l[i] + warm_tol {
                    zl[i] = v0;
                }
                if x_u[i] < inner_x_u[i] - warm_tol {
                    zu[i] = v0;
                }
            }
            (zl, zu)
        } else {
            (vec![0.0; n], vec![0.0; n])
        };

        // Phase 2: detect redundant linear rows in the (possibly
        // tightened) box. Non-linear rows are never dropped. The
        // `row_kept_inner` mask was initialised above by Phase 0.
        let mut n_dropped_rows: Index = 0;
        let kept_before_phase2 = row_kept_inner.clone();
        if self.opts.redundant_constraint_removal {
            let redundant_mask = find_redundant_rows(&linear_rows, &x_l, &x_u, 1e-9);
            // `redundant_mask` aligns with the *kept* linear rows
            // (`linear_rows`); map it back onto inner rows, skipping
            // rows Phase 0 already dropped (C1).
            n_dropped_rows =
                apply_redundant_verdicts(&linear_row_map, &redundant_mask, &mut row_kept_inner);
        }

        let singleton_bounds = find_singleton_bounds(
            &linear_row_map,
            &kept_before_phase2,
            &row_kept_inner,
            &x_l,
            &x_u,
            &inner_x_l,
            &inner_x_u,
        );
        let phase2_dropped_rows: Vec<(usize, Vec<(usize, Number)>)> = linear_row_map
            .iter()
            .enumerate()
            .filter(|(i, _)| {
                kept_before_phase2.get(*i).copied().unwrap_or(false)
                    && !row_kept_inner.get(*i).copied().unwrap_or(true)
                    && !singleton_bounds.iter().any(|sb| sb.row == *i)
            })
            .filter_map(|(i, lr)| {
                lr.as_ref().map(|lr| {
                    (
                        i,
                        lr.entries
                            .iter()
                            .filter(|(_, a)| *a != 0.0 && a.is_finite())
                            .map(|&(j, a)| (j as usize, a))
                            .collect(),
                    )
                })
            })
            .collect();

        // Phase 3: structural LICQ check on the kept equality rows.
        let licq_verdict = if self.opts.licq_check {
            let eq_tol: Number = 1e-12;
            let mut eq_rows: Vec<EqRow> = Vec::new();
            for (i, &kept) in row_kept_inner.iter().enumerate() {
                if !kept {
                    continue;
                }
                if (g_u_inner[i] - g_l_inner[i]).abs() > eq_tol {
                    continue;
                }
                use std::collections::BTreeSet;
                let mut cols: BTreeSet<Index> = BTreeSet::new();
                for &(j, v) in &by_row[i] {
                    if v != 0.0 {
                        cols.insert(j);
                    }
                }
                eq_rows.push(EqRow {
                    cols: cols.into_iter().collect(),
                });
            }
            Some(licq_check(&eq_rows, info_inner.n))
        } else {
            None
        };

        // Build outer row mapping.
        let mut rows_kept: Vec<usize> = Vec::with_capacity(m_in);
        let mut row_inner_to_outer = vec![usize::MAX; m_in];
        for (i, &kept) in row_kept_inner.iter().enumerate() {
            if kept {
                row_inner_to_outer[i] = rows_kept.len();
                rows_kept.push(i);
            }
        }
        let m_out = rows_kept.len();

        // Build outer Jacobian sparsity: keep entries whose row is
        // kept, renumber rows.
        let mut jac_kept_idx = Vec::new();
        let mut jac_irow_outer = Vec::new();
        let mut jac_jcol_outer = Vec::new();
        for k in 0..nnz_in {
            let i_inner = if one_based {
                (jac_irow_inner[k] - 1) as usize
            } else {
                jac_irow_inner[k] as usize
            };
            if i_inner >= m_in {
                continue;
            }
            if !row_kept_inner[i_inner] {
                continue;
            }
            let outer = row_inner_to_outer[i_inner];
            let outer_row_index = if one_based {
                (outer as Index) + 1
            } else {
                outer as Index
            };
            jac_irow_outer.push(outer_row_index);
            jac_jcol_outer.push(jac_jcol_inner[k]);
            jac_kept_idx.push(k);
        }
        let nnz_out = jac_kept_idx.len();

        // Reduced g_l, g_u in outer ordering.
        let g_l: Vec<Number> = rows_kept.iter().map(|&i| g_l_inner[i]).collect();
        let g_u: Vec<Number> = rows_kept.iter().map(|&i| g_u_inner[i]).collect();

        let info_outer = NlpInfo {
            n: info_inner.n,
            m: m_out as Index,
            nnz_jac_g: nnz_out as Index,
            // Linear rows contribute zero to the Hessian, so dropping
            // them does not change `nnz_h_lag`. We carry the inner
            // sparsity through unchanged.
            nnz_h_lag: info_inner.nnz_h_lag,
            index_style: info_inner.index_style,
        };

        self.state = Some(PresolveState {
            info_inner,
            info_outer,
            bounds: CachedBounds { x_l, x_u, g_l, g_u },
            rows_kept,
            jac_kept_idx,
            jac_irow_outer,
            jac_jcol_outer,
            tighten_report,
            certified_infeasible,
            fbbt_report,
            n_dropped_rows,
            licq_verdict,
            z_l_warm,
            z_u_warm,
            scratch_g: vec![0.0; m_in],
            scratch_jac: vec![0.0; nnz_in],
            scratch_lambda: vec![0.0; m_in],
            aux_diagnostics,
            reduction_stack,
            singleton_bounds,
            decl_x_l: inner_x_l,
            decl_x_u: inner_x_u,
            decl_g_l: g_l_inner,
            decl_g_u: g_u_inner,
            phase2_dropped_rows,
        });
        self.state.as_ref()
    }
}

/// Map `find_redundant_rows`' verdict mask back onto `row_kept_inner`,
/// returning the number of rows newly dropped.
///
/// The mask is aligned to the **kept** linear rows in inner-row order
/// — i.e. it has exactly one entry per `(i)` where `linear_row_map[i]`
/// is `Some` **and** `row_kept_inner[i]` is still `true` (that is the
/// filter `linear_rows` was built with). The mapping iterator must
/// therefore advance only on rows that are both linear and still kept.
///
/// Regression guard for C1: the original loop advanced the mask on
/// every `Some` row regardless of `row_kept_inner`, so once Phase 0 had
/// dropped any linear row every later verdict landed on its
/// predecessor's row — silently deleting a binding constraint (and
/// keeping a redundant one).
fn apply_redundant_verdicts(
    linear_row_map: &[Option<LinearRow>],
    redundant_mask: &[bool],
    row_kept_inner: &mut [bool],
) -> Index {
    let mut mask = redundant_mask.iter();
    let mut n_dropped: Index = 0;
    for (i, lr) in linear_row_map.iter().enumerate() {
        if lr.is_some() && row_kept_inner[i] {
            if *mask.next().unwrap_or(&false) {
                row_kept_inner[i] = false;
                n_dropped += 1;
            }
        }
    }
    n_dropped
}

// Inside this impl every `.expect("inited")` is invariant-protected
// by the preceding `ensure_init` (which is the only way state ever
// becomes `Some`).
#[allow(clippy::expect_used)]
impl TNLP for PresolveTnlp {
    fn is_presolve_wrapper(&self) -> bool {
        true
    }

    fn presolve_infeasibility_proof(&self) -> Option<InfeasibilityProof> {
        self.certified_infeasible()
    }

    fn get_nlp_info(&mut self) -> Option<NlpInfo> {
        let s = self.ensure_init()?;
        Some(s.info_outer)
    }

    fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
        let Some(s) = self.ensure_init() else {
            return false;
        };
        b.x_l.copy_from_slice(&s.bounds.x_l);
        b.x_u.copy_from_slice(&s.bounds.x_u);
        b.g_l.copy_from_slice(&s.bounds.g_l);
        b.g_u.copy_from_slice(&s.bounds.g_u);
        true
    }

    fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
        // n is unchanged by presolve; lambda warm-start is masked to
        // kept rows only (caller already sized `sp.lambda` to m_out).
        self.served_warm_start = None;
        let Some(_) = self.ensure_init() else {
            return false;
        };
        // For now, ask the inner TNLP for its starting point in full
        // and project lambda down. Most users don't warm-start
        // duals, so this hits the no-op path.
        let m_in = self.state.as_ref().expect("inited").info_inner.m as usize;
        // Start from the caller's buffers, not zeros: `OrigIpoptNlp` fills
        // them with NaN to mean "unseeded" (gh#622), and an inner TNLP that
        // supplies no bound multipliers must leave that marker standing. A
        // zero fill turned a primal-only warm start into "z = 0 supplied".
        let mut z_l_full = sp.z_l.to_vec();
        let mut z_u_full = sp.z_u.to_vec();
        let mut lambda_full = vec![0.0; m_in];
        let ok = self.inner.borrow_mut().get_starting_point(StartingPoint {
            init_x: sp.init_x,
            x: sp.x,
            init_z: sp.init_z,
            z_l: &mut z_l_full,
            z_u: &mut z_u_full,
            init_lambda: sp.init_lambda,
            lambda: &mut lambda_full,
        });
        if !ok {
            return false;
        }
        // What the caller handed, before any seeding below (`warm_echo`).
        let handed = (sp.init_z || sp.init_lambda)
            .then(|| (z_l_full.clone(), z_u_full.clone(), lambda_full.clone()));
        let s = self.state.as_ref().expect("inited");
        // M24, inbound: a dropped singleton row's supplied multiplier is the
        // bound multiplier of the bound presolve made from it — the inverse
        // of the hand-back in `finalize_solution`, so a `.sol` round-trips.
        if sp.init_z && sp.init_lambda {
            for sb in &s.singleton_bounds {
                let Some(&lam) = lambda_full.get(sb.row) else {
                    continue;
                };
                // Contribution to stationarity at the variable: `+a·λ`,
                // which the bound carries as `+z_u` or `−z_l`.
                let c = sb.coeff * lam;
                if !c.is_finite() {
                    continue;
                }
                let (side, z) = if sb.upper {
                    (&mut z_u_full, c)
                } else {
                    (&mut z_l_full, -c)
                };
                if z > 0.0
                    && let Some(dst) = side.get_mut(sb.var)
                {
                    *dst = if dst.is_finite() { *dst + z } else { z };
                }
            }
            // M24, derived bounds, inbound: a multi-variable row Phase 2
            // dropped is not in the reduced problem, so its supplied
            // multiplier's contribution to stationarity, `a_ij·λ_i` at each
            // variable it touches, can only be carried there by the reduced
            // bound on that variable — the derived bound the row helped
            // imply. Seed it so the reduced start is as stationary as the
            // supplied one; `finalize_solution`'s recovery moves it back.
            let n = z_l_full.len();
            let mut net = vec![0.0; n];
            for (row, entries) in &s.phase2_dropped_rows {
                let Some(&lam) = lambda_full.get(*row) else {
                    continue;
                };
                if !lam.is_finite() || lam == 0.0 {
                    continue;
                }
                for &(j, a) in entries {
                    if let Some(v) = net.get_mut(j) {
                        *v += a * lam;
                    }
                }
            }
            for (j, &c) in net.iter().enumerate() {
                if !c.is_finite() || c == 0.0 {
                    continue;
                }
                // `+c` in stationarity is carried as `+z_u` or `−z_l`.
                let (side, z, present) = if c > 0.0 {
                    (&mut z_u_full, c, upper_bound_present(s.bounds.x_u[j]))
                } else {
                    (&mut z_l_full, -c, lower_bound_present(s.bounds.x_l[j]))
                };
                if present && let Some(dst) = side.get_mut(j) {
                    *dst = if dst.is_finite() { *dst + z } else { z };
                }
            }
        }
        sp.z_l.copy_from_slice(&z_l_full);
        sp.z_u.copy_from_slice(&z_u_full);
        let mut x_fixed_overridden_count = None;
        let mut x_clamped_count = None;
        // Project the served primal into the reduced box the solver
        // sees, mirroring `warm::project_warm_point`.
        if self.project_seed && sp.init_x {
            let mut n_fixed = 0;
            for frame in s.reduction_stack.iter_bottom_up() {
                for (k, &i) in frame.fixed_vars.iter().enumerate() {
                    if let (Some(dst), Some(&v)) = (sp.x.get_mut(i), frame.fixed_values.get(k)) {
                        if *dst != v {
                            n_fixed += 1;
                        }
                        *dst = v;
                    }
                }
            }
            let mut n_clamped = 0;
            for (i, v) in sp.x.iter_mut().enumerate() {
                if let (Some(&lo), Some(&hi)) = (s.bounds.x_l.get(i), s.bounds.x_u.get(i))
                    && let Some(c) = crate::warm::clamp_seed(*v, lo, hi)
                {
                    *v = c;
                    n_clamped += 1;
                }
            }
            x_fixed_overridden_count = Some(n_fixed);
            x_clamped_count = Some(n_clamped);
        }
        // Phase 4: overlay presolve hints onto any zero/unset
        // entries. User-supplied warm-start values always win.
        if sp.init_z && self.opts.warm_z_bounds {
            for (i, &hint) in s.z_l_warm.iter().enumerate() {
                if hint > 0.0 && sp.z_l[i] <= 0.0 {
                    sp.z_l[i] = hint;
                }
            }
            for (i, &hint) in s.z_u_warm.iter().enumerate() {
                if hint > 0.0 && sp.z_u[i] <= 0.0 {
                    sp.z_u[i] = hint;
                }
            }
        }
        for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
            sp.lambda[outer] = lambda_full[i_inner];
        }
        let n_dropped_rows = m_in.saturating_sub(s.rows_kept.len());
        let dropped_dual_l1 = sp.init_lambda.then(|| {
            let mut kept = vec![false; m_in];
            for &inner in &s.rows_kept {
                if inner < kept.len() {
                    kept[inner] = true;
                }
            }
            lambda_full
                .iter()
                .enumerate()
                .filter(|(i, _)| !kept[*i])
                .map(|(_, value)| value.abs())
                .sum()
        });
        self.starting_point_projection_report.n_dropped_rows = n_dropped_rows;
        if let Some(value) = dropped_dual_l1 {
            self.starting_point_projection_report.dropped_dual_l1 = value;
        }
        if let Some(value) = x_fixed_overridden_count {
            self.starting_point_projection_report
                .x_fixed_overridden_count = value;
        }
        if let Some(value) = x_clamped_count {
            self.starting_point_projection_report.x_clamped_count = value;
        }
        // Record what was served, so `finalize_solution` can recognise it
        // coming back unchanged. An entry the solver ignores (`init_z` /
        // `init_lambda` off) is recorded as NaN, i.e. "not a seed".
        let served = handed.map(|(z_l_in, z_u_in, lambda_in)| warm_echo::ServedWarmStart {
            z_l_in,
            z_u_in,
            lambda_in,
            z_l_out: warm_echo::served_bound_multipliers(
                sp.init_z,
                sp.z_l,
                &s.bounds.x_l,
                &s.bounds.x_u,
            ),
            z_u_out: warm_echo::served_bound_multipliers(
                sp.init_z,
                sp.z_u,
                &s.bounds.x_l,
                &s.bounds.x_u,
            ),
            lambda_out: warm_echo::served_row_multipliers(sp.init_lambda, sp.lambda),
        });
        self.served_warm_start = served;
        true
    }

    fn eval_f(&mut self, x: &[Number], new_x: bool) -> Option<Number> {
        self.inner.borrow_mut().eval_f(x, new_x)
    }

    fn eval_grad_f(&mut self, x: &[Number], new_x: bool, grad_f: &mut [Number]) -> bool {
        let ok = self.inner.borrow_mut().eval_grad_f(x, new_x, grad_f);
        if ok && self.live_record.enabled(&self.inner) {
            self.live_record.grad_f(x, grad_f);
        }
        ok
    }

    fn eval_g(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        let s = self.state.as_mut().expect("inited");
        if !self.inner.borrow_mut().eval_g(x, new_x, &mut s.scratch_g) {
            return false;
        }
        if self.live_record.enabled(&self.inner) {
            self.live_record.g(x, &s.scratch_g);
        }
        for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
            g[outer] = s.scratch_g[i_inner];
        }
        true
    }

    fn eval_jac_g(&mut self, x: Option<&[Number]>, new_x: bool, mode: SparsityRequest<'_>) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        match mode {
            SparsityRequest::Structure { irow, jcol } => {
                let s = self.state.as_ref().expect("inited");
                irow.copy_from_slice(&s.jac_irow_outer);
                jcol.copy_from_slice(&s.jac_jcol_outer);
                true
            }
            SparsityRequest::Values { values } => {
                let s = self.state.as_mut().expect("inited");
                if !self.inner.borrow_mut().eval_jac_g(
                    x,
                    new_x,
                    SparsityRequest::Values {
                        values: &mut s.scratch_jac,
                    },
                ) {
                    return false;
                }
                if let Some(x) = x
                    && self.live_record.enabled(&self.inner)
                {
                    self.live_record.jac(x, &s.scratch_jac);
                }
                for (outer_k, &inner_k) in s.jac_kept_idx.iter().enumerate() {
                    values[outer_k] = s.scratch_jac[inner_k];
                }
                true
            }
        }
    }

    fn eval_h(
        &mut self,
        x: Option<&[Number]>,
        new_x: bool,
        obj_factor: Number,
        lambda: Option<&[Number]>,
        new_lambda: bool,
        mode: SparsityRequest<'_>,
    ) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        // Hessian sparsity is untouched: linear rows (the only ones
        // we drop) contribute zero. Forward `lambda` after expanding
        // outer → inner with zeros at dropped rows.
        let lambda_full_opt = if let Some(lam) = lambda {
            let s = self.state.as_mut().expect("inited");
            for v in s.scratch_lambda.iter_mut() {
                *v = 0.0;
            }
            for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
                s.scratch_lambda[i_inner] = lam[outer];
            }
            Some(&s.scratch_lambda[..])
        } else {
            None
        };
        // Re-borrow inner after dropping the state borrow.
        let lam_ref: Option<&[Number]> = lambda_full_opt;
        // SAFETY: `lam_ref` borrows from `self.state`'s scratch; the
        // call to `inner.borrow_mut()` does not touch `self.state`.
        self.inner
            .borrow_mut()
            .eval_h(x, new_x, obj_factor, lam_ref, new_lambda, mode)
    }

    fn finalize_solution(&mut self, sol: Solution<'_>, ip_data: &IpoptData, ip_cq: &IpoptCq) {
        let Some(_) = self.ensure_init() else {
            // Init failed earlier — best effort: just forward as-is.
            self.inner
                .borrow_mut()
                .finalize_solution(sol, ip_data, ip_cq);
            return;
        };
        // Reconstruct inner-sized g and lambda.
        let (g_full, mut lambda_full, n_inner, m_inner, nnz_inner, one_based) = {
            let s = self.state.as_mut().expect("inited");
            // Recompute g at sol.x — the solver gave us reduced g.
            let ok_g = self
                .inner
                .borrow_mut()
                .eval_g(sol.x, true, &mut s.scratch_g);
            if !ok_g {
                // L45: the final constraint re-eval failed, so `scratch_g` is
                // unreliable — it holds whatever the failing call left behind
                // (partial garbage, or a stale value from an earlier iterate).
                // Don't forward that. Rebuild g from the solver's own
                // (trustworthy) reduced `sol.g`, mapped back to the kept inner
                // rows exactly as the multiplier mapping below does; rows
                // dropped by presolve are left at 0 (no reliable value exists
                // for them once the re-eval has failed).
                for v in s.scratch_g.iter_mut() {
                    *v = 0.0;
                }
                for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
                    if i_inner < s.scratch_g.len() && outer < sol.g.len() {
                        s.scratch_g[i_inner] = sol.g[outer];
                    }
                }
            }
            for v in s.scratch_lambda.iter_mut() {
                *v = 0.0;
            }
            for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
                s.scratch_lambda[i_inner] = sol.lambda[outer];
            }
            (
                s.scratch_g.clone(),
                s.scratch_lambda.clone(),
                s.info_inner.n as usize,
                s.info_inner.m as usize,
                s.info_inner.nnz_jac_g as usize,
                matches!(s.info_inner.index_style, IndexStyle::Fortran),
            )
        };

        // A warm start the solver handed back without moving its duals: the
        // original problem's duals are the ones the caller supplied, not a
        // re-derivation from stationarity at the pushed `x` (`warm_echo`).
        if let Some((lambda_e, z_l_e, z_u_e)) = self.echoed_warm_start(&sol) {
            tracing::debug!(
                target: "pounce_presolve",
                "presolve postsolve: the solver returned the served warm start unchanged; \
                 reporting the supplied duals"
            );
            self.bound_multipliers_handed_to_rows.clear();
            self.derived_bound_recovery = None;
            self.forward_finalized(
                sol,
                &g_full,
                &lambda_e,
                &lambda_full,
                &z_l_e,
                &z_u_e,
                ip_data,
                ip_cq,
            );
            return;
        }

        // Phase-0 (issue #53) multiplier recovery for dropped rows.
        // Walk the reduction stack top-down; for each frame, compute
        // the full-space ∇f and J at sol.x, then solve the k×k LU
        // system to recover λ for the frame's dropped rows. Splice
        // the result into `lambda_full` at the dropped indices.
        let frames: Vec<reduction_frame::ReductionFrame> = {
            let s = self.state.as_ref().expect("inited");
            s.reduction_stack.iter_top_down().cloned().collect()
        };
        // Bound multipliers forwarded to the inner finalize. At each frame's
        // Phase-0 aux-fixed variables these are zeroed once the variable's
        // dropped-row multipliers are recovered below:
        // `recover_dropped_multipliers` folds the entire stationarity residual
        // into λ under the documented assumption `z_l = z_u = 0` there, so
        // forwarding the IPM's clamp multipliers unchanged would double-count
        // that contribution and break `∇f − Jᵀλ − z_l + z_u = 0`.
        let mut z_l_full = sol.z_l.to_vec();
        let mut z_u_full = sol.z_u.to_vec();
        if !frames.is_empty() && m_inner > 0 {
            let mut grad_f = vec![0.0; n_inner];
            let ok_grad = self
                .inner
                .borrow_mut()
                .eval_grad_f(sol.x, true, &mut grad_f);
            // Sparsity + values for the full inner Jacobian.
            let mut jac_irow_inner = vec![0 as Index; nnz_inner];
            let mut jac_jcol_inner = vec![0 as Index; nnz_inner];
            let ok_struct = if nnz_inner > 0 {
                self.inner.borrow_mut().eval_jac_g(
                    None,
                    false,
                    SparsityRequest::Structure {
                        irow: &mut jac_irow_inner,
                        jcol: &mut jac_jcol_inner,
                    },
                )
            } else {
                true
            };
            let mut jac_values = vec![0.0; nnz_inner];
            let ok_vals = if nnz_inner > 0 {
                self.inner.borrow_mut().eval_jac_g(
                    Some(sol.x),
                    false,
                    SparsityRequest::Values {
                        values: &mut jac_values,
                    },
                )
            } else {
                true
            };
            if ok_grad && ok_struct && ok_vals {
                // `recover_dropped_multipliers` reads the Jacobian only at the
                // frames' fixed-var columns, so materialize just those columns
                // (the union across all frames) instead of the full
                // `m_inner × n_inner` dense block — O(m·k) rather than O(m·n),
                // where k = total distinct fixed vars is tiny next to n. The
                // old full densification cost ~80 GB at 100k×100k (issue M26).
                let mut orig_to_compact = vec![usize::MAX; n_inner];
                let mut n_cols = 0usize;
                for frame in &frames {
                    for &c in &frame.fixed_vars {
                        if c < n_inner && orig_to_compact[c] == usize::MAX {
                            orig_to_compact[c] = n_cols;
                            n_cols += 1;
                        }
                    }
                }
                let mut jac_cols = vec![0.0; m_inner * n_cols];
                for k in 0..nnz_inner {
                    let i = if one_based {
                        (jac_irow_inner[k] as isize - 1) as usize
                    } else {
                        jac_irow_inner[k] as usize
                    };
                    let j = if one_based {
                        (jac_jcol_inner[k] as isize - 1) as usize
                    } else {
                        jac_jcol_inner[k] as usize
                    };
                    if i < m_inner && j < n_inner {
                        let cc = orig_to_compact[j];
                        if cc != usize::MAX {
                            jac_cols[i * n_cols + cc] = jac_values[k];
                        }
                    }
                }
                for frame in &frames {
                    if let Ok(lam_dropped) = frame.recover_dropped_multipliers_cols(
                        &grad_f,
                        &jac_cols,
                        n_cols,
                        &orig_to_compact,
                        &lambda_full,
                    ) {
                        for (idx, &r) in frame.dropped_rows.iter().enumerate() {
                            lambda_full[r] = lam_dropped[idx];
                        }
                        // The recovered λ now carries this frame's fixed-var
                        // stationarity; drop the IPM's clamp multipliers there.
                        for &i in &frame.fixed_vars {
                            if i < z_l_full.len() {
                                z_l_full[i] = 0.0;
                            }
                            if i < z_u_full.len() {
                                z_u_full[i] = 0.0;
                            }
                        }
                    }
                }
            }
        }
        let lambda_lifted = lambda_full.clone();
        // M24: a bound that exists only because a dropped singleton row
        // implied it hands its multiplier back to that row.
        self.bound_multipliers_handed_to_rows.clear();
        {
            let s = self.state.as_ref().expect("inited");
            for sb in &s.singleton_bounds {
                let side = if sb.upper {
                    &mut z_u_full
                } else {
                    &mut z_l_full
                };
                let (Some(z), Some(lam)) = (side.get_mut(sb.var), lambda_full.get_mut(sb.row))
                else {
                    continue;
                };
                if !z.is_finite() || *z == 0.0 {
                    continue;
                }
                let signed = if sb.upper { *z } else { -*z };
                *lam += signed / sb.coeff;
                *z = 0.0;
                self.bound_multipliers_handed_to_rows
                    .push((sb.var, sb.upper));
            }
        }
        // M24, derived bounds: a multiplier the solver put on a bound
        // presolve derived from rows, where the original model has no
        // active bound, is not a multiplier of the original problem.
        // Recover the original problem's in its place (`dual_recovery`).
        self.derived_bound_recovery = self.recover_derived_bound_duals(
            sol.x,
            &g_full,
            &mut lambda_full,
            &mut z_l_full,
            &mut z_u_full,
        );
        self.forward_finalized(
            sol,
            &g_full,
            &lambda_full,
            &lambda_lifted,
            &z_l_full,
            &z_u_full,
            ip_data,
            ip_cq,
        );
    }

    fn get_var_con_metadata(&mut self, var: &mut MetaData, con: &mut MetaData) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        // Variable count is unchanged by presolve, so var metadata
        // flows through. Constraint metadata is per-inner-row; if we
        // dropped rows, subset the per-row vectors to kept rows.
        let mut inner_var = MetaData::default();
        let mut inner_con = MetaData::default();
        if !self
            .inner
            .borrow_mut()
            .get_var_con_metadata(&mut inner_var, &mut inner_con)
        {
            return false;
        }
        *var = inner_var;
        let s = self.state.as_ref().expect("inited");
        let m_in = s.info_inner.m as usize;
        *con = project_con_metadata(&inner_con, &s.rows_kept, m_in);
        true
    }

    fn get_scaling_parameters(&mut self, req: ScalingRequest<'_>) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        let s = self.state.as_ref().expect("inited");
        let m_in = s.info_inner.m as usize;
        // Allocate inner-sized g_scaling and forward.
        let mut inner_g = vec![1.0; m_in];
        let mut use_x = false;
        let mut use_g = false;
        let mut obj_scaling = 1.0;
        let inner_x_scaling_len = req.x_scaling.len();
        let mut inner_x = vec![1.0; inner_x_scaling_len];
        let ok = self
            .inner
            .borrow_mut()
            .get_scaling_parameters(ScalingRequest {
                obj_scaling: &mut obj_scaling,
                use_x_scaling: &mut use_x,
                x_scaling: &mut inner_x,
                use_g_scaling: &mut use_g,
                g_scaling: &mut inner_g,
            });
        if !ok {
            return false;
        }
        *req.obj_scaling = obj_scaling;
        *req.use_x_scaling = use_x;
        *req.use_g_scaling = use_g;
        req.x_scaling.copy_from_slice(&inner_x);
        for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
            req.g_scaling[outer] = inner_g[i_inner];
        }
        true
    }

    fn get_variables_linearity(&mut self, types: &mut [Linearity]) -> bool {
        self.inner.borrow_mut().get_variables_linearity(types)
    }

    fn get_objective_variables_linearity(&mut self, types: &mut [Linearity]) -> bool {
        self.inner
            .borrow_mut()
            .get_objective_variables_linearity(types)
    }

    fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
        let Some(_) = self.ensure_init() else {
            return false;
        };
        let m_in = self.state.as_ref().expect("inited").info_inner.m as usize;
        let mut full = vec![Linearity::NonLinear; m_in];
        if !self.inner.borrow_mut().get_constraints_linearity(&mut full) {
            return false;
        }
        let s = self.state.as_ref().expect("inited");
        for (outer, &i_inner) in s.rows_kept.iter().enumerate() {
            types[outer] = full[i_inner];
        }
        true
    }

    fn get_number_of_nonlinear_variables(&mut self) -> Index {
        self.inner.borrow_mut().get_number_of_nonlinear_variables()
    }

    fn get_list_of_nonlinear_variables(&mut self, pos_nonlin_vars: &mut [Index]) -> bool {
        self.inner
            .borrow_mut()
            .get_list_of_nonlinear_variables(pos_nonlin_vars)
    }

    fn wants_intermediate_callback(&self) -> bool {
        self.inner.borrow().wants_intermediate_callback()
    }

    fn wants_current_iterate(&self) -> bool {
        self.inner.borrow().wants_current_iterate()
    }

    fn intermediate_callback(
        &mut self,
        stats: IterStats,
        ip_data: &IpoptData,
        ip_cq: &IpoptCq,
    ) -> bool {
        // The live iterate arrives in this layer's (row-dropped) space;
        // hand the inner TNLP its own, lifted by `finalize_solution`'s
        // rule. The lift runs that function, so the state it records about
        // the *final* solution is put back afterwards.
        let lifted = match self.get_nlp_info() {
            Some(info) => {
                // Every field `finalize_solution` writes, including the
                // derived-bound dual recovery's, so a lift never leaves a
                // live iterate's record where the final solution's belongs.
                let saved = (
                    self.finalized_full_solution.take(),
                    std::mem::take(&mut self.bound_multipliers_handed_to_rows),
                    std::mem::take(&mut self.bound_multiplier_overrides),
                    self.duals_rewritten,
                    self.derived_bound_recovery.take(),
                );
                // Columns are the inner ones; outer row `k` is inner
                // `rows_kept[k]`.
                let cols: Vec<usize> = (0..info.n as usize).collect();
                let rows: Vec<usize> = self
                    .state
                    .as_ref()
                    .map(|s| s.rows_kept.clone())
                    .unwrap_or_else(|| (0..info.m as usize).collect());
                let lifted = live_lift::lift(
                    self,
                    |s| &mut s.inner,
                    (info.n as usize, info.m as usize),
                    // Every column is kept, so no declared bound is needed.
                    live_lift::LiftMaps {
                        cols: &cols,
                        rows: &rows,
                        declared_x_l: &[],
                        declared_x_u: &[],
                        record: self.live_record.clone(),
                    },
                );
                (
                    self.finalized_full_solution,
                    self.bound_multipliers_handed_to_rows,
                    self.bound_multiplier_overrides,
                    self.duals_rewritten,
                    self.derived_bound_recovery,
                ) = saved;
                lifted
            }
            None => None,
        };
        let _lifted = lifted.map(pounce_nlp::live_iterate::UserIterateGuard::install);
        self.inner
            .borrow_mut()
            .intermediate_callback(stats, ip_data, ip_cq)
    }

    fn finalize_metadata(&mut self, var: &MetaData, con: &MetaData) {
        let Some(_) = self.ensure_init() else {
            self.inner.borrow_mut().finalize_metadata(var, con);
            return;
        };
        let s = self.state.as_ref().expect("inited");
        let m_in = s.info_inner.m as usize;
        let con_full = expand_con_metadata(con, &s.rows_kept, m_in);
        self.inner.borrow_mut().finalize_metadata(var, &con_full);
    }
}

/// Subset every per-row vector of `inner` to the rows in `rows_kept`.
///
/// Per-row-ness is inferred from `v.len() == m_in`, because `MetaData`
/// (mirroring upstream Ipopt's untyped `*MetaDataMapType`) carries no
/// per-key arity tag — so the length is the only signal available. By
/// the TNLP contract every vector in the *constraint* `MetaData` bucket
/// is per-constraint (length `m_in`), so this is exact for any
/// conforming inner TNLP (e.g. `nl_reader`'s per-row `con_names`).
///
/// L46 caveat: a vector that is semantically *global* yet happens to
/// have length `m_in` would be wrongly subset, and a genuinely per-row
/// vector whose length differs from `m_in` (a contract violation) is
/// passed through unchanged. Neither can be distinguished from the data
/// alone; fully resolving it needs an arity-aware metadata API. The
/// fast path `m_in == rows_kept.len()` (no rows dropped) is identity,
/// so the hazard only exists when presolve actually drops a row.
fn project_con_metadata(inner: &MetaData, rows_kept: &[usize], m_in: usize) -> MetaData {
    let mut out = MetaData::default();
    for (k, v) in &inner.strings {
        out.strings.insert(
            k.clone(),
            if v.len() == m_in {
                rows_kept.iter().map(|&i| v[i].clone()).collect()
            } else {
                v.clone()
            },
        );
    }
    for (k, v) in &inner.integers {
        out.integers.insert(
            k.clone(),
            if v.len() == m_in {
                rows_kept.iter().map(|&i| v[i]).collect()
            } else {
                v.clone()
            },
        );
    }
    for (k, v) in &inner.numerics {
        out.numerics.insert(
            k.clone(),
            if v.len() == m_in {
                rows_kept.iter().map(|&i| v[i]).collect()
            } else {
                v.clone()
            },
        );
    }
    out
}

/// Expand every per-(outer-row) vector back to `m_in` rows by
/// inserting empty / 0 / 0.0 defaults at dropped rows. The inverse of
/// [`project_con_metadata`]; per-row-ness is inferred from
/// `v.len() == m_out` and carries the same L46 caveat (see that
/// function): a coincidentally-`m_out`-length global vector is wrongly
/// expanded, inherent to the untyped `MetaData` API.
fn expand_con_metadata(outer: &MetaData, rows_kept: &[usize], m_in: usize) -> MetaData {
    let m_out = rows_kept.len();
    let mut full = MetaData::default();
    for (k, v) in &outer.strings {
        let mut buf: Vec<String> = vec![String::new(); m_in];
        if v.len() == m_out {
            for (outer_i, val) in v.iter().enumerate() {
                buf[rows_kept[outer_i]] = val.clone();
            }
            full.strings.insert(k.clone(), buf);
        } else {
            full.strings.insert(k.clone(), v.clone());
        }
    }
    for (k, v) in &outer.integers {
        let mut buf: Vec<Index> = vec![0; m_in];
        if v.len() == m_out {
            for (outer_i, &val) in v.iter().enumerate() {
                buf[rows_kept[outer_i]] = val;
            }
            full.integers.insert(k.clone(), buf);
        } else {
            full.integers.insert(k.clone(), v.clone());
        }
    }
    for (k, v) in &outer.numerics {
        let mut buf: Vec<Number> = vec![0.0; m_in];
        if v.len() == m_out {
            for (outer_i, &val) in v.iter().enumerate() {
                buf[rows_kept[outer_i]] = val;
            }
            full.numerics.insert(k.clone(), buf);
        } else {
            full.numerics.insert(k.clone(), v.clone());
        }
    }
    full
}

/// Re-export for callers that already imported
/// `pounce_presolve::register_options` directly.
pub fn register(reg: &RegisteredOptions) -> Result<(), SolverException> {
    register_options(reg)
}

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

    /// **gh #402.** A crossing between an absent bound and a real one past the
    /// *opposite* sentinel is not a crossing, and must never certify.
    ///
    /// The test was a bare `l > u`, so `x_l = -1e19` (absent) against
    /// `x_u = -5e20` (real) read as a `5e20` gap — far too large for
    /// `is_negligible` to dismiss — and certified a feasible model as *proved*
    /// infeasible, the strongest claim POUNCE makes.
    #[test]
    fn a_sentinel_against_a_real_bound_never_certifies_a_crossing() {
        let tol = 1e-8;
        assert!(
            !crossing_is_certifiable(&[-1e19], &[-5e20], tol),
            "`x <= -5e20` with no lower bound is feasible, not a crossed box"
        );
        assert!(
            !crossing_is_certifiable(&[5e20], &[1e19], tol),
            "`x >= 5e20` with no upper bound is feasible too"
        );
        // A genuine crossing between two present bounds still certifies.
        assert!(crossing_is_certifiable(&[5.0], &[3.0], tol));
        // And a sub-tolerance one still does not (the #380 rule).
        assert!(!crossing_is_certifiable(&[1.0 + 1e-13], &[1.0], tol));
    }

    /// **gh #402.** The acceptance slack must track the row's own scale, even
    /// when the row's real bound lies past the opposite sentinel.
    ///
    /// The old symmetric test contributed `0.0` for a bound of `-5e20`, so the
    /// margin collapsed to `tol * 1.0` — a widening of `1e-8` on a row written
    /// at `5e20`. Any infeasibility survives that, which is the fail-open
    /// direction this margin exists to prevent.
    #[test]
    fn the_row_margin_tracks_a_bound_past_the_opposite_sentinel() {
        let tol = 1e-8;
        // `g <= -5e20`, no lower bound. Magnitude is 5e20, not 1.
        let m = row_margin_for(-1e19, -5e20, tol);
        assert!(
            (m - tol * 5e20).abs() <= tol * 5e20 * 1e-12,
            "margin should be tol*5e20 = {}, got {m}",
            tol * 5e20
        );
        // Absent on both sides -> the floor of 1.
        assert_eq!(row_margin_for(-1e19, 1e19, tol), tol);
        // An ordinary two-sided row uses the larger magnitude.
        assert_eq!(row_margin_for(-2.0, 7.0, tol), tol * 7.0);
    }

    struct Probe;
    impl TNLP for Probe {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 0,
                nnz_jac_g: 0,
                nnz_h_lag: 1,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, _b: BoundsInfo<'_>) -> bool {
            true
        }
        fn get_starting_point(&mut self, _sp: StartingPoint<'_>) -> bool {
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, _g: &mut [Number]) -> bool {
            true
        }
        fn eval_g(&mut self, _x: &[Number], _new_x: bool, _g: &mut [Number]) -> bool {
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            _mode: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    struct NoopProvider;
    impl ExpressionProvider for NoopProvider {}

    #[test]
    fn c1_redundancy_mask_realigned_after_phase0_drop() {
        // Regression for C1. `find_redundant_rows` returns a verdict
        // mask aligned to the *kept* linear rows (inner-row order).
        // When Phase 0 has dropped an earlier linear row, the mapping
        // back to inner rows must skip already-dropped rows; advancing
        // the mask iterator on them shifts every later verdict onto its
        // predecessor's row.
        let lr = |lo: Number| {
            Some(LinearRow {
                entries: vec![(0, 1.0)],
                lo,
                hi: lo,
            })
        };
        // Three linear rows; Phase 0 dropped inner row 0.
        let linear_row_map = vec![lr(0.0), lr(1.0), lr(2.0)];
        let row_kept = vec![false, true, true];
        // `find_redundant_rows` ran over the *kept* rows [row1, row2]
        // and flagged the SECOND kept row (inner row 2) redundant; the
        // first kept row (inner row 1) is binding.
        let mask = vec![false, true];

        // Correct mapping: only inner row 2 is dropped.
        let mut kept_new = row_kept.clone();
        let n = apply_redundant_verdicts(&linear_row_map, &mask, &mut kept_new);
        assert_eq!(
            kept_new,
            vec![false, true, false],
            "verdict must land on inner row 2, not its predecessor"
        );
        assert_eq!(n, 1);

        // Document the pre-fix misalignment this test guards against:
        // the old loop advanced the mask on every `Some` row regardless
        // of `row_kept_inner`, dropping inner row 1 (binding!) and
        // keeping inner row 2 (redundant) — a silent wrong answer.
        let mut buggy = row_kept.clone();
        let mut it = mask.iter();
        for (i, l) in linear_row_map.iter().enumerate() {
            if l.is_some() && *it.next().unwrap_or(&false) {
                buggy[i] = false;
            }
        }
        assert_eq!(buggy, vec![false, false, true]);
        assert_ne!(buggy, kept_new);
    }

    #[test]
    fn disabled_returns_inner_unchanged() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Probe));
        let opts = PresolveOptions {
            enabled: false,
            ..PresolveOptions::defaults()
        };
        let wrapped = wrap_with_presolve(Rc::clone(&inner), opts).unwrap();
        assert!(Rc::ptr_eq(&inner, &wrapped));
    }

    #[test]
    fn enabled_wraps_and_forwards() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Probe));
        let opts = PresolveOptions {
            enabled: true,
            ..PresolveOptions::defaults()
        };
        let wrapped = wrap_with_presolve(Rc::clone(&inner), opts).unwrap();
        assert!(!Rc::ptr_eq(&inner, &wrapped));
        let info = wrapped.borrow_mut().get_nlp_info().unwrap();
        assert_eq!(info.n, 1);
        assert_eq!(info.m, 0);
    }

    #[test]
    fn provider_wrapper_does_not_nest_an_existing_presolve_wrapper() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Probe));
        let opts = PresolveOptions {
            enabled: true,
            ..PresolveOptions::defaults()
        };
        let manual = wrap_with_presolve(inner, opts.clone()).unwrap();
        let provider: Rc<RefCell<dyn ExpressionProvider>> = Rc::new(RefCell::new(NoopProvider));
        let wrapped = wrap_with_presolve_provider(manual.clone(), provider, opts).unwrap();

        assert!(Rc::ptr_eq(&manual, &wrapped));
    }

    #[test]
    fn register_options_roundtrip() {
        let reg = RegisteredOptions::default();
        register_options(&reg).unwrap();
        let opt = reg.get_option("presolve").expect("presolve registered");
        assert_eq!(opt.name, "presolve");
    }

    #[test]
    fn auxiliary_phase0_noop_when_disabled() {
        // Master switch off → Phase 0 returns zero diagnostics and
        // does not perturb the inner problem's dimensions.
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Probe));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: false,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.n, 1);
        assert_eq!(info.m, 0);
        let diag = wrapped.auxiliary_diagnostics();
        assert_eq!(diag.blocks_eliminated, 0);
        assert_eq!(diag.vars_eliminated, 0);
        assert_eq!(diag.rows_eliminated, 0);
    }

    #[test]
    fn auxiliary_phase0_noop_when_enabled_no_algos_yet() {
        // For a 1-var, 0-row TNLP there are no candidate blocks; the
        // orchestrator skips its work regardless of master-switch state.
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Probe));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Aggressive,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.n, 1);
        assert_eq!(info.m, 0);
        let diag = wrapped.auxiliary_diagnostics();
        assert_eq!(diag.blocks_eliminated, 0);
        assert!(diag.rejection_reasons.is_empty());
    }

    /// 2-variable, 2-equality TNLP that PR 8's orchestrator should
    /// reduce: `x + y = 3, x - y = 1` → unique solution `(2, 1)`.
    /// Zero objective gradient → PureEquality, eligible under Safe.
    struct TwoVarSquareEq;
    impl TNLP for TwoVarSquareEq {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 2,
                nnz_jac_g: 4,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            for v in b.x_l.iter_mut() {
                *v = -1e19;
            }
            for v in b.x_u.iter_mut() {
                *v = 1e19;
            }
            b.g_l[0] = 3.0;
            b.g_u[0] = 3.0;
            b.g_l[1] = 1.0;
            b.g_u[1] = 1.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 0.0;
                sp.x[1] = 0.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for v in g.iter_mut() {
                *v = 0.0;
            }
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] + x[1];
            g[1] = x[0] - x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 0, 1, 1]);
                    jcol.copy_from_slice(&[0, 1, 0, 1]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0, 1.0, -1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    #[test]
    fn phase0_via_tnlp_eliminates_square_block() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(TwoVarSquareEq));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        wrapped.set_project_seed(true);
        let info = wrapped.get_nlp_info().expect("init ok");
        // Variable count unchanged (clamp, not reduce).
        assert_eq!(info.n, 2);
        // Both equality rows dropped.
        assert_eq!(info.m, 0);

        let diag = wrapped.auxiliary_diagnostics();
        assert_eq!(diag.blocks_eliminated, 1);
        assert_eq!(diag.vars_eliminated, 2);
        assert_eq!(diag.rows_eliminated, 2);

        let bounds = wrapped.cached_bounds().expect("inited");
        assert!((bounds.x_l[0] - 2.0).abs() < 1e-12);
        assert!((bounds.x_u[0] - 2.0).abs() < 1e-12);
        assert!((bounds.x_l[1] - 1.0).abs() < 1e-12);
        assert!((bounds.x_u[1] - 1.0).abs() < 1e-12);
    }

    /// [`TwoVarSquareEq`] with linearity tags layered on: the *global*
    /// tags report both variables `NonLinear` (as a TNLP whose equality
    /// rows are nonlinear would — the gas-network shape), while the
    /// objective-scoped tags, when offered, report both `Linear`
    /// (matching the zero objective gradient).
    struct SquareEqWithTags {
        inner: TwoVarSquareEq,
        offer_obj_tags: bool,
    }
    impl TNLP for SquareEqWithTags {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            self.inner.get_nlp_info()
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            self.inner.get_bounds_info(b)
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            self.inner.get_starting_point(sp)
        }
        fn eval_f(&mut self, x: &[Number], new_x: bool) -> Option<Number> {
            self.inner.eval_f(x, new_x)
        }
        fn eval_grad_f(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            self.inner.eval_grad_f(x, new_x, g)
        }
        fn eval_g(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            self.inner.eval_g(x, new_x, g)
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            self.inner.eval_jac_g(x, new_x, mode)
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            self.inner.get_constraints_linearity(types)
        }
        fn get_variables_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::NonLinear;
            types[1] = Linearity::NonLinear;
            true
        }
        fn get_objective_variables_linearity(&mut self, types: &mut [Linearity]) -> bool {
            if !self.offer_obj_tags {
                return false;
            }
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, sol: Solution<'_>, d: &IpoptData, q: &IpoptCq) {
            self.inner.finalize_solution(sol, d, q)
        }
    }

    /// Regression for the gas-network CI failure: objective-scoped tags
    /// must take precedence over the global ones, so a variable that is
    /// nonlinear only in *constraints* does not poison `obj_support` and
    /// block a legitimate Phase-0 elimination. Pre-fix (global tags only)
    /// both variables were unioned into the objective support and the
    /// Safe policy kept both rows.
    #[test]
    fn phase0_objective_scoped_tags_dont_block_constraint_nonlinear_elimination() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(SquareEqWithTags {
            inner: TwoVarSquareEq,
            offer_obj_tags: true,
        }));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(
            info.m, 0,
            "objective-free block must still be eliminated when only the \
             global tags say NonLinear"
        );
        assert_eq!(wrapped.auxiliary_diagnostics().vars_eliminated, 2);
    }

    /// When the TNLP offers only the global tags, the fallback stays
    /// conservative: constraint nonlinearity is treated as possible
    /// objective coupling and the block is kept (sound, just suboptimal).
    #[test]
    fn phase0_global_tags_fallback_remains_conservative() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(SquareEqWithTags {
            inner: TwoVarSquareEq,
            offer_obj_tags: false,
        }));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(
            info.m, 2,
            "global-only tags fall back to the conservative union"
        );
        assert_eq!(wrapped.auxiliary_diagnostics().vars_eliminated, 0);
    }

    #[test]
    fn phase0_via_tnlp_disabled_is_pass_through() {
        // Same inner TNLP, but presolve_auxiliary=no. Orchestrator
        // is byte-identical to today; both rows remain.
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(TwoVarSquareEq));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: false,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.n, 2);
        assert_eq!(info.m, 2);
        let diag = wrapped.auxiliary_diagnostics();
        assert_eq!(diag.blocks_eliminated, 0);
    }

    /// Smoke test: turning on `presolve_auxiliary_diagnostics` still
    /// produces a correct elimination. We can't easily capture
    /// stderr from inside a #[test], but we can verify the option
    /// path doesn't break the orchestrator.
    #[test]
    fn phase0_via_tnlp_diagnostics_flag_does_not_break_solve() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(TwoVarSquareEq));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_diagnostics: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.m, 0);
        let diag = wrapped.auxiliary_diagnostics();
        assert_eq!(diag.blocks_eliminated, 1);
    }

    /// Regression for the #60 PR review blocker. The original
    /// `TwoVarSquareEq` TNLP, when wrapped with BOTH
    /// `presolve_auxiliary=yes` AND `presolve_bound_tightening=yes`
    /// (the new defaults), used to:
    ///   (1) let aux clamp x_l[0..2] = x_u[0..2] = solved values;
    ///   (2) let `tighten_bounds` re-propagate the (dropped)
    ///       equality rows over the clamped bounds, derive a
    ///       contradiction, and set `tighten_report.infeasible`;
    ///   (3) hand `x_l > x_u` to the IPM (because `infeasible` was
    ///       never inspected), crashing it with
    ///       `Invalid Problem Definition`.
    /// The fix: build `linear_rows` AFTER Phase 0, filtered by
    /// `row_kept_inner`, so dropped rows don't propagate. With this
    /// test we just confirm the wrapper init still succeeds and the
    /// final bounds remain self-consistent.
    #[test]
    fn phase0_via_tnlp_no_infeasible_with_default_bound_tightening() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(TwoVarSquareEq));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            // Both phases on — this is the combination that crashed
            // on `gaslib11_steady.nl` before the fix.
            bound_tightening: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.m, 0); // aux dropped both rows
        let bounds = wrapped.cached_bounds().expect("inited");
        // x_l ≤ x_u must hold for every variable.
        for i in 0..(info.n as usize) {
            assert!(
                bounds.x_l[i] <= bounds.x_u[i] + 1e-12,
                "x_l[{i}] = {} > x_u[{i}] = {}",
                bounds.x_l[i],
                bounds.x_u[i]
            );
        }
        // tighten_report must NOT have flagged infeasibility.
        let rpt = wrapped.tighten_report();
        assert!(!rpt.infeasible, "Phase 1 falsely flagged infeasibility");
    }

    /// Same model as [`TwoVarSquareEq`] but records the bound multipliers
    /// `z_l`/`z_u` that `finalize_solution` forwards to the inner TNLP, so a
    /// test can assert what the reported KKT point carries.
    struct RecordingTwoVar {
        rec: Rc<RefCell<Option<(Vec<Number>, Vec<Number>)>>>,
    }
    impl TNLP for RecordingTwoVar {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 2,
                nnz_jac_g: 4,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            for v in b.x_l.iter_mut() {
                *v = -1e19;
            }
            for v in b.x_u.iter_mut() {
                *v = 1e19;
            }
            b.g_l[0] = 3.0;
            b.g_u[0] = 3.0;
            b.g_l[1] = 1.0;
            b.g_u[1] = 1.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 0.0;
                sp.x[1] = 0.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for v in g.iter_mut() {
                *v = 0.0;
            }
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] + x[1];
            g[1] = x[0] - x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 0, 1, 1]);
                    jcol.copy_from_slice(&[0, 1, 0, 1]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0, 1.0, -1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
            *self.rec.borrow_mut() = Some((sol.z_l.to_vec(), sol.z_u.to_vec()));
        }
    }

    /// Regression for H10. Phase 0 clamps `x_l = x_u = v` at the two
    /// eliminated variables, so the IPM reports large bound multipliers
    /// there. `recover_dropped_multipliers` attributes the *whole*
    /// stationarity residual to the recovered row multipliers λ under the
    /// documented assumption `z_l = z_u = 0` at those variables; if
    /// `finalize_solution` forwards the clamp multipliers unchanged the
    /// contribution is double-counted and the reported point violates
    /// `∇f − Jᵀλ − z_l + z_u = 0`. The fix zeros `z_l`/`z_u` at every
    /// frame's `fixed_vars` once their λ has been recovered.
    #[test]
    fn phase0_finalize_zeroes_bound_multipliers_at_fixed_vars() {
        let rec = Rc::new(RefCell::new(None));
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(RecordingTwoVar {
            rec: Rc::clone(&rec),
        }));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        wrapped.set_project_seed(true);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.n, 2, "variable count unchanged (clamp, not reduce)");
        assert_eq!(info.m, 0, "both equality rows dropped by Phase 0");

        let (mut x0, mut z_l0, mut z_u0) = (vec![0.0; 2], vec![0.0; 2], vec![0.0; 2]);
        let mut lambda0 = Vec::new();
        wrapped.reset_starting_point_projection_report();
        assert!(wrapped.get_starting_point(StartingPoint {
            init_x: true,
            x: &mut x0,
            init_z: false,
            z_l: &mut z_l0,
            z_u: &mut z_u0,
            init_lambda: true,
            lambda: &mut lambda0,
        }));
        assert_eq!(x0, vec![2.0, 1.0], "wrapper serves aux-fixed values");
        let start_report = wrapped.starting_point_projection_report();
        assert_eq!(
            start_report.x_fixed_overridden_count, 2,
            "report = {start_report:?}"
        );
        assert_eq!(start_report.x_clamped_count, 0, "report = {start_report:?}");
        assert_eq!(start_report.n_dropped_rows, 2, "report = {start_report:?}");

        // The reduced solve returns the clamped point (2, 1); the IPM places
        // large bound multipliers on the now-fixed variables, and the reduced
        // problem has no rows (so empty g / lambda).
        let x = [2.0, 1.0];
        let z_l = [7.0, 0.0];
        let z_u = [0.0, 3.0];
        let g: [Number; 0] = [];
        let lambda: [Number; 0] = [];
        let sol = Solution {
            status: pounce_nlp::alg_types::SolverReturn::Success,
            x: &x,
            z_l: &z_l,
            z_u: &z_u,
            g: &g,
            lambda: &lambda,
            obj_value: 0.0,
        };
        wrapped.finalize_solution(sol, &IpoptData::default(), &IpoptCq::default());

        let (got_zl, got_zu) = rec.borrow().clone().expect("inner finalize ran");
        // fixed_vars = {0, 1}: both must be zeroed. Pre-fix the inner sees the
        // forwarded clamp multipliers [7,0] / [0,3] verbatim.
        assert_eq!(
            got_zl,
            vec![0.0, 0.0],
            "z_l must be zeroed at aux-fixed vars (H10)"
        );
        assert_eq!(
            got_zu,
            vec![0.0, 0.0],
            "z_u must be zeroed at aux-fixed vars (H10)"
        );
    }

    /// min x^2 s.t. x + y = 3, x/y in [0, 5]. The row tightens both
    /// upper bounds 5 -> 3; the seed sits at the stale corner [5, 5].
    struct TightenSeed {
        x0: Vec<Number>,
    }

    impl TNLP for TightenSeed {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 1,
                nnz_jac_g: 2,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l.copy_from_slice(&[0.0, 0.0]);
            b.x_u.copy_from_slice(&[5.0, 5.0]);
            b.g_l[0] = 3.0;
            b.g_u[0] = 3.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x.copy_from_slice(&self.x0);
            }
            true
        }
        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
            Some(x[0] * x[0])
        }
        fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 2.0 * x[0];
            g[1] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] + x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 0]);
                    jcol.copy_from_slice(&[0, 1]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types.fill(Linearity::Linear);
            true
        }
        fn finalize_solution(&mut self, _s: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// The report is recorded by the same call that projects the served seed.
    #[test]
    fn served_warm_primal_records_its_projection() {
        let inner: Rc<RefCell<dyn TNLP>> =
            Rc::new(RefCell::new(TightenSeed { x0: vec![5.0, 5.0] }));
        let opts = PresolveOptions {
            enabled: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        wrapped.set_project_seed(true);

        let (mut x_l, mut x_u) = (vec![0.0; 2], vec![0.0; 2]);
        let (mut g_l, mut g_u) = (vec![0.0; 1], vec![0.0; 1]);
        assert!(wrapped.get_bounds_info(BoundsInfo {
            x_l: &mut x_l,
            x_u: &mut x_u,
            g_l: &mut g_l,
            g_u: &mut g_u,
        }));
        assert_eq!(x_u, vec![3.0, 3.0], "row x + y = 3 tightens x_u to 3");

        let (mut x, mut z_l, mut z_u) = (vec![0.0; 2], vec![0.0; 2], vec![0.0; 2]);
        let mut lambda = vec![0.0; 1];
        assert!(wrapped.get_starting_point(StartingPoint {
            init_x: true,
            x: &mut x,
            init_z: false,
            z_l: &mut z_l,
            z_u: &mut z_u,
            init_lambda: false,
            lambda: &mut lambda,
        }));
        assert_eq!(
            x,
            vec![3.0, 3.0],
            "served seed must be clamped into the tightened box"
        );

        let report = wrapped.starting_point_projection_report();
        assert_eq!(report.x_clamped_count, 2, "report = {report:?}");
        assert_eq!(report.x_fixed_overridden_count, 0, "report = {report:?}");
    }

    /// Same model as [`RecordingTwoVar`] but (a) records the `g` vector that
    /// `finalize_solution` forwards to the inner TNLP and (b) can be told to
    /// *fail* `eval_g` (returning `false` after scribbling a sentinel), so a
    /// test can check what happens when the final constraint re-eval fails.
    struct GFailRecordingVar {
        rec_g: Rc<RefCell<Option<Vec<Number>>>>,
        fail_g: Rc<RefCell<bool>>,
    }
    impl TNLP for GFailRecordingVar {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 2,
                nnz_jac_g: 4,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            for v in b.x_l.iter_mut() {
                *v = -1e19;
            }
            for v in b.x_u.iter_mut() {
                *v = 1e19;
            }
            b.g_l[0] = 3.0;
            b.g_u[0] = 3.0;
            b.g_l[1] = 1.0;
            b.g_u[1] = 1.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 0.0;
                sp.x[1] = 0.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for v in g.iter_mut() {
                *v = 0.0;
            }
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            if *self.fail_g.borrow() {
                // Simulate a failing evaluator that leaves garbage behind.
                for v in g.iter_mut() {
                    *v = 999.0;
                }
                return false;
            }
            g[0] = x[0] + x[1];
            g[1] = x[0] - x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 0, 1, 1]);
                    jcol.copy_from_slice(&[0, 1, 0, 1]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0, 1.0, -1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
            *self.rec_g.borrow_mut() = Some(sol.g.to_vec());
        }
    }

    /// L45: when the final constraint re-eval in `finalize_solution` fails,
    /// the forwarded `g` must NOT be the garbage/stale buffer the failing
    /// `eval_g` left behind. Both rows here are Phase-0 dropped (reduced
    /// `m = 0`, empty `sol.g`), so the trustworthy fallback yields all-zeros
    /// — never the `999.0` sentinel.
    #[test]
    fn finalize_does_not_forward_stale_g_when_eval_g_fails() {
        let rec_g = Rc::new(RefCell::new(None));
        let fail_g = Rc::new(RefCell::new(false));
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(GFailRecordingVar {
            rec_g: Rc::clone(&rec_g),
            fail_g: Rc::clone(&fail_g),
        }));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.m, 0, "both equality rows dropped by Phase 0");

        // Now make the final re-eval fail.
        *fail_g.borrow_mut() = true;
        let x = [2.0, 1.0];
        let z_l = [0.0, 0.0];
        let z_u = [0.0, 0.0];
        let g: [Number; 0] = [];
        let lambda: [Number; 0] = [];
        let sol = Solution {
            status: pounce_nlp::alg_types::SolverReturn::Success,
            x: &x,
            z_l: &z_l,
            z_u: &z_u,
            g: &g,
            lambda: &lambda,
            obj_value: 0.0,
        };
        wrapped.finalize_solution(sol, &IpoptData::default(), &IpoptCq::default());

        let got_g = rec_g.borrow().clone().expect("inner finalize ran");
        // Fail-first: pre-fix the ignored `eval_g` return forwards the
        // sentinel-filled `scratch_g` (`[999, 999]`) verbatim.
        assert_eq!(
            got_g,
            vec![0.0, 0.0],
            "failed eval_g must not forward stale/garbage constraint values",
        );
    }

    /// gh #396's shape, reduced to its essentials: one `<=`-only row whose
    /// genuine upper bound is *more negative* than the `-INF_BOUND` sentinel
    /// that stands in for its absent lower bound, and a starting point sitting
    /// exactly on that upper bound.
    ///
    /// `-1e30 · x <= -5e20` over `x ∈ [0, 1e-9]`, with `x0 = 5e-10`, so
    /// `g(x0) == g_u` exactly. True violation: zero.
    struct SentinelLowerBoundRow;
    impl TNLP for SentinelLowerBoundRow {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 1,
                nnz_jac_g: 1,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = 0.0;
            b.x_u[0] = 1e-9;
            b.g_l[0] = -1e19; // absent — the sentinel
            b.g_u[0] = -5e20; // real, and beyond the sentinel's magnitude
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 5e-10;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = -1e30 * x[0];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => {
                    values[0] = -1e30;
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// **gh #396.** The absent-bound sentinel must not manufacture a violation,
    /// and a genuine bound beyond the sentinel's magnitude must not be
    /// discarded.
    ///
    /// Both errors landed on this one row pre-fix. `viol` had no presence test,
    /// so it scored `g_l - v = 4.9e20` against a lower bound the row does not
    /// have; and the magnitude test `|b| < INF_BOUND` called the row's *real*
    /// upper bound `-5e20` absent. Either alone sinks the witness, and with no
    /// candidate able to pass, presolve certified a feasible model infeasible
    /// with `solve_result_num` 201 — a claimed proof.
    #[test]
    fn absent_bound_sentinel_does_not_manufacture_a_violation() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(SentinelLowerBoundRow));
        assert!(
            witness_refutes_infeasibility(
                &inner,
                1,
                1,
                &[0.0],
                &[1e-9],
                &[-1e19],
                &[-5e20],
                1e-8,
                WitnessRule::SolverAcceptance,
            ),
            "x0 = 5e-10 puts the row exactly on its only real bound — that is a \
             witness, and the -1e19 sentinel standing in for the absent lower \
             bound must not turn it into a 4.9e20 violation"
        );
    }

    /// The other direction, and the reason the fix is a *directional* presence
    /// test rather than dropping the terms: a point that genuinely violates a
    /// real bound must still fail, so correct verdicts survive.
    ///
    /// Tightening the ceiling to `-5.1e20` over `x ∈ [4e-10, 5e-10]` puts every
    /// candidate — including `x0` — strictly above it, so there is genuinely no
    /// witness and the verdict must stand.
    ///
    /// This one does *not* reproduce #396: pre-fix it also returned `false`,
    /// via the phantom lower-bound violation rather than the real upper-bound
    /// one. It is here to pin the direction the fix could plausibly have broken
    /// — making `up_present` admit the bound must not make real violations
    /// invisible.
    #[test]
    fn a_real_bound_beyond_the_sentinel_still_blocks_the_witness() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(SentinelLowerBoundRow));
        assert!(
            !witness_refutes_infeasibility(
                &inner,
                1,
                1,
                &[4e-10],
                &[5e-10],
                &[-1e19],
                &[-5.1e20],
                1e-8,
                WitnessRule::SolverAcceptance,
            ),
            "no point in [4e-10, 5e-10] satisfies -1e30*x <= -5.1e20, so there \
             is no witness — a real bound of magnitude 5.1e20 must not be \
             discarded as the infinity sentinel"
        );
    }

    /// `min x` over `x ∈ [decl_l, 10]` with one linear row `a·x` bounded so
    /// that it reads `x ≥ 2` for either sign of `a` — the M24 shape: Phase 1
    /// turns the row into `x_l = 2` and Phase 2 drops it as redundant.
    struct SingletonLowerRow {
        a: Number,
        decl_l: Number,
    }
    impl TNLP for SingletonLowerRow {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 1,
                nnz_jac_g: 1,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = self.decl_l;
            b.x_u[0] = 10.0;
            if self.a > 0.0 {
                b.g_l[0] = 2.0 * self.a;
                b.g_u[0] = 1e19;
            } else {
                b.g_l[0] = -1e19;
                b.g_u[0] = 2.0 * self.a;
            }
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 5.0;
            }
            if sp.init_lambda {
                // The row multiplier a `presolve=no` solve reports:
                // `∇f + a·λ = 0` with `∇f = 1`.
                sp.lambda[0] = -1.0 / self.a;
            }
            true
        }
        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
            Some(x[0])
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 1.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = self.a * x[0];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => {
                    values[0] = self.a;
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// Finalize `SingletonLowerRow` at its optimum `x = 2` with the reduced
    /// solve's `z_l = 1` on the presolve-made bound, returning the full-space
    /// row multiplier and the handed-back list.
    fn finalize_singleton(a: Number, decl_l: Number) -> (Number, Vec<(usize, bool)>) {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(SingletonLowerRow { a, decl_l }));
        let opts = PresolveOptions {
            enabled: true,
            bound_tightening: true,
            redundant_constraint_removal: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(inner, opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(wrapped.n_dropped_rows(), 1, "a={a}: the row is dropped");
        assert_eq!(info.m, 0);
        let b = wrapped.cached_bounds().expect("inited");
        assert!((b.x_l[0] - 2.0).abs() < 1e-12, "a={a}: x_l = {}", b.x_l[0]);
        let (x, z_l, z_u) = ([2.0], [1.0], [0.0]);
        let (g, lambda): ([Number; 0], [Number; 0]) = ([], []);
        wrapped.finalize_solution(
            Solution {
                status: pounce_nlp::alg_types::SolverReturn::Success,
                x: &x,
                z_l: &z_l,
                z_u: &z_u,
                g: &g,
                lambda: &lambda,
                obj_value: 2.0,
            },
            &IpoptData::default(),
            &IpoptCq::default(),
        );
        let (_, lamf) = wrapped.finalized_full_solution().expect("finalized");
        assert_eq!(lamf.len(), 1, "full-space lambda regains the original row");
        (lamf[0], wrapped.bound_multipliers_handed_to_rows().to_vec())
    }

    /// M24 (see the module-level "Dual-attribution caveat"): the multiplier
    /// of a bound presolve made from a singleton row it then dropped goes
    /// back to the row, as a `presolve=no` solve reports it. Ipopt's
    /// convention, `∇f + Jᵀλ − z_l + z_u = 0`: `1 + a·λ − 0 = 0`, so
    /// `λ = −1/a`, for either sign of `a` (with `a < 0` the row's *upper*
    /// side implies the variable's *lower* bound).
    #[test]
    fn dropped_singleton_row_takes_back_its_bound_dual() {
        for a in [1.0, 4.0, -1.0, -0.5] {
            let (lam, handed) = finalize_singleton(a, 0.0);
            assert!(
                (lam - (-1.0 / a)).abs() < 1e-12,
                "a={a}: dropped-row λ = {lam}, expected {}",
                -1.0 / a
            );
            assert_eq!(handed, vec![(0, false)], "a={a}: z_l[0] was handed back");
            // Stationarity with the bound multiplier zeroed.
            let stat = 1.0 + a * lam;
            assert!(
                stat.abs() < 1e-12,
                "a={a}: KKT stationarity residual {stat}"
            );
        }
    }

    /// A declared bound that already sits where the row puts it keeps its
    /// multiplier: the two are a degenerate pair and the tie goes to the
    /// incumbent, as for `linear_eq_elim`'s transferred bounds.
    #[test]
    fn a_declared_bound_tied_with_the_row_keeps_its_dual() {
        let (lam, handed) = finalize_singleton(1.0, 2.0);
        assert_eq!(lam, 0.0, "the row is redundant against the declared bound");
        assert!(handed.is_empty(), "{handed:?}");
    }

    /// The inverse, for a warm start: the dropped row's supplied multiplier
    /// becomes the bound multiplier of the bound presolve made from it, so
    /// a `.sol` written through presolve round-trips.
    #[test]
    fn a_warm_row_multiplier_seeds_the_bound_presolve_made_from_it() {
        for a in [1.0, -0.5] {
            let inner: Rc<RefCell<dyn TNLP>> =
                Rc::new(RefCell::new(SingletonLowerRow { a, decl_l: 0.0 }));
            let opts = PresolveOptions {
                enabled: true,
                bound_tightening: true,
                redundant_constraint_removal: true,
                ..PresolveOptions::defaults()
            };
            let mut wrapped = PresolveTnlp::new(inner, opts);
            let _ = wrapped.get_nlp_info().expect("init ok");
            let (mut x, mut z_l, mut z_u) = ([0.0], [Number::NAN], [Number::NAN]);
            let mut lambda: [Number; 0] = [];
            assert!(wrapped.get_starting_point(StartingPoint {
                init_x: true,
                x: &mut x,
                init_z: true,
                z_l: &mut z_l,
                z_u: &mut z_u,
                init_lambda: true,
                lambda: &mut lambda,
            }));
            assert!((z_l[0] - 1.0).abs() < 1e-12, "a={a}: z_l = {}", z_l[0]);
            assert!(z_u[0].is_nan(), "a={a}: the upper side stays unseeded");
        }
    }

    /// `min (x0+1)² + (x1−3)²  s.t.  C0: x0 = 0.5,  C1: x0 + x1 ≤ 2`, no
    /// declared bounds. Presolve fixes x0 from C0, derives `x1 ≤ 1.5` from
    /// C1, and drops both rows (M24, derived bounds). At `x* = (0.5, 1.5)`
    /// the original multipliers are unique: `λ = (−6, 3)`, no bound
    /// multipliers.
    struct DerivedChain;
    impl TNLP for DerivedChain {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 2,
                nnz_jac_g: 3,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l.fill(-1e19);
            b.x_u.fill(1e19);
            b.g_l.copy_from_slice(&[0.5, -1e19]);
            b.g_u.copy_from_slice(&[0.5, 2.0]);
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x.copy_from_slice(&[0.5, 1.5]);
            }
            if sp.init_lambda {
                sp.lambda.copy_from_slice(&[-6.0, 3.0]);
            }
            true
        }
        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
            Some((x[0] + 1.0).powi(2) + (x[1] - 3.0).powi(2))
        }
        fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 2.0 * (x[0] + 1.0);
            g[1] = 2.0 * (x[1] - 3.0);
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0];
            g[1] = x[0] + x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 1, 1]);
                    jcol.copy_from_slice(&[0, 0, 1]);
                }
                SparsityRequest::Values { values } => values.copy_from_slice(&[1.0, 1.0, 1.0]),
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types.fill(Linearity::Linear);
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    fn derived_chain() -> PresolveTnlp {
        let opts = PresolveOptions {
            enabled: true,
            bound_tightening: true,
            redundant_constraint_removal: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::new(RefCell::new(DerivedChain)), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.m, 0, "both rows are dropped");
        let b = wrapped.cached_bounds().expect("inited");
        assert!(
            (b.x_u[1] - 1.5).abs() < 1e-12,
            "x1 ≤ 1.5 derived: {}",
            b.x_u[1]
        );
        wrapped
    }

    /// Outbound: the reduced solve's `z_U(x1) = 3` on the derived bound and
    /// x0's fixed-variable multiplier go back to C1 and C0.
    #[test]
    fn a_derived_bound_multiplier_is_recovered_onto_its_rows() {
        let mut wrapped = derived_chain();
        // Reduced stationarity: x0 (fixed) 3 − z_l = 0; x1 −3 + z_u = 0.
        let (x, z_l, z_u) = ([0.5, 1.5], [3.0, 0.0], [0.0, 3.0]);
        let (g, lambda): ([Number; 0], [Number; 0]) = ([], []);
        wrapped.finalize_solution(
            Solution {
                status: pounce_nlp::alg_types::SolverReturn::Success,
                x: &x,
                z_l: &z_l,
                z_u: &z_u,
                g: &g,
                lambda: &lambda,
                obj_value: 4.5,
            },
            &IpoptData::default(),
            &IpoptCq::default(),
        );
        let (_, lam) = wrapped.finalized_full_solution().expect("finalized");
        assert!(
            (lam[0] + 6.0).abs() < 1e-9 && (lam[1] - 3.0).abs() < 1e-9,
            "λ = {lam:?}"
        );
        assert_eq!(wrapped.derived_bound_dual_recovery(), Some((1, true)));
        let ov = wrapped.bound_multiplier_overrides().to_vec();
        assert!(
            ov.iter().all(|&(_, l, u)| l.abs() < 1e-9 && u.abs() < 1e-9) && ov.len() == 2,
            "{ov:?}"
        );
        assert!(wrapped.postsolve_rewrote_duals());
    }

    /// Inbound: C1's supplied multiplier seeds the derived bound it implies,
    /// so the reduced start is stationary where the supplied one was.
    #[test]
    fn a_dropped_rows_warm_multiplier_seeds_the_derived_bound() {
        let mut wrapped = derived_chain();
        let (mut x, mut z_l, mut z_u) = ([0.0; 2], [Number::NAN; 2], [Number::NAN; 2]);
        let mut lambda: [Number; 0] = [];
        assert!(wrapped.get_starting_point(StartingPoint {
            init_x: true,
            x: &mut x,
            init_z: true,
            z_l: &mut z_l,
            z_u: &mut z_u,
            init_lambda: true,
            lambda: &mut lambda,
        }));
        assert!((z_u[1] - 3.0).abs() < 1e-12, "z_U(x1) = {}", z_u[1]);
        assert!(z_l[1].is_nan(), "x1's lower side stays unseeded");
        // x0: C0 (singleton, λ0 = −6) seeds z_L = 6, C1 adds z_U = 3; the
        // reduced stationarity there, 3 − 6 + 3, closes.
        assert!(
            (z_l[0] - 6.0).abs() < 1e-12 && (z_u[0] - 3.0).abs() < 1e-12,
            "{z_l:?} {z_u:?}"
        );
    }

    /// M25 scratch: one variable, box `x ∈ [0, 10]`, two contradictory linear
    /// rows `x ≥ 5` and `x ≤ 3`. Phase 1 tightens `x_l = 5`, `x_u = 3` and
    /// flags infeasible. Auxiliary is off, so the reduction stack is empty and
    /// the rollback guard does NOT fire.
    struct TwoContradictoryRows;
    impl TNLP for TwoContradictoryRows {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 2,
                nnz_jac_g: 2,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = 0.0;
            b.x_u[0] = 10.0;
            b.g_l[0] = 5.0; // x ≥ 5
            b.g_u[0] = 1e19;
            b.g_l[1] = -1e19; // x ≤ 3
            b.g_u[1] = 3.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 4.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0];
            g[1] = x[0];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 1]);
                    jcol.copy_from_slice(&[0, 0]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// Regression for M25. A genuine Phase-1 infeasibility found with an empty
    /// reduction stack (auxiliary off, so the aux rollback guard never fires)
    /// must NOT hand crossed bounds `x_l > x_u` to the IPM — that triggers an
    /// invalid-problem failure instead of a clean infeasibility verdict.
    /// Presolve restores the original box and lets the IPM certify
    /// infeasibility itself; the `infeasible` flag is still surfaced for
    /// diagnostics.
    #[test]
    fn phase1_infeasible_restores_valid_box_for_ipm() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(TwoContradictoryRows));
        let opts = PresolveOptions {
            enabled: true,
            bound_tightening: true,
            auxiliary: false,
            ..PresolveOptions::defaults()
        };
        let mut wrapped = PresolveTnlp::new(Rc::clone(&inner), opts);
        let _info = wrapped.get_nlp_info().expect("init ok");

        // The infeasibility is still detected and surfaced for diagnostics.
        assert!(
            wrapped.tighten_report().infeasible,
            "Phase 1 must still flag the empty feasible region"
        );

        // ...and, with no Phase-0 elimination in force, it is a genuine
        // contradiction in the model rather than a presolve artifact, so it is
        // certified. This is what lets the solver report "proved infeasible"
        // instead of re-deriving a weaker numerical verdict.
        assert_eq!(
            wrapped.certified_infeasible(),
            Some(InfeasibilityProof::BoundPropagation),
            "a Phase-1 contradiction on an un-clamped box must be certified"
        );

        // But the box handed to the IPM is the original, valid box — not the
        // crossed `[5, 3]` Phase 1 derived.
        let b = wrapped.cached_bounds().expect("inited");
        assert!(
            b.x_l[0] <= b.x_u[0] + 1e-12,
            "M25: bounds handed to IPM must be valid, got x_l={} > x_u={}",
            b.x_l[0],
            b.x_u[0]
        );
        assert!(
            (b.x_l[0] - 0.0).abs() < 1e-12 && (b.x_u[0] - 10.0).abs() < 1e-12,
            "box restored to the original [0, 10], got [{}, {}]",
            b.x_l[0],
            b.x_u[0]
        );
    }

    /// 1-variable, 2-constraint TNLP that drives FBBT into a partial
    /// tightening followed by a genuine infeasibility. Both constraints are
    /// `g = x` (default `NonLinear` linearity, so they are FBBT-handled, not
    /// turned into linear rows Phase 1 would catch first). Row 0 demands
    /// `x ∈ [0.3, 0.7]` (tightens the box `x ∈ [0, 1]`); row 1 demands
    /// `x = 5` (infeasible). FBBT applies row 0, then row 1 flags the box
    /// empty — leaving the bounds in the partially-tightened `[0.3, 0.7]`
    /// state the `FbbtReport` contract says must not be trusted.
    struct FbbtPartialThenInfeasible;
    impl TNLP for FbbtPartialThenInfeasible {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 2,
                nnz_jac_g: 2,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = 0.0;
            b.x_u[0] = 1.0;
            b.g_l[0] = 0.3;
            b.g_u[0] = 0.7;
            b.g_l[1] = 5.0;
            b.g_u[1] = 5.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 0.5;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0];
            g[1] = x[0];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 1]);
                    jcol.copy_from_slice(&[0, 0]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0]);
                }
            }
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// Returns the tape `x` (i.e. `Var(0)`) for both constraints.
    struct VarTapeProvider;
    impl ExpressionProvider for VarTapeProvider {
        fn constraint_expression(
            &self,
            i: usize,
        ) -> Option<pounce_nlp::expression_provider::FbbtTape> {
            use pounce_nlp::expression_provider::{FbbtOp, FbbtTape};
            if i < 2 {
                Some(FbbtTape {
                    ops: vec![FbbtOp::Var(0)],
                })
            } else {
                None
            }
        }
    }

    /// H12: when FBBT detects a genuine infeasibility, the partially
    /// tightened (and per its own contract "undefined") bounds must NOT
    /// reach the reduced problem. Presolve has no infeasibility channel, so
    /// it discards FBBT's bounds and lets the IPM run on the pre-FBBT box —
    /// here the original `[0, 1]`, never the corrupted `[0.3, 0.7]`.
    #[test]
    fn fbbt_infeasibility_discards_corrupted_bounds() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(FbbtPartialThenInfeasible));
        let provider: Rc<RefCell<dyn ExpressionProvider>> = Rc::new(RefCell::new(VarTapeProvider));
        let opts = PresolveOptions {
            enabled: true,
            fbbt: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped =
            PresolveTnlp::with_expression_provider(Rc::clone(&inner), Rc::clone(&provider), opts);
        let info = wrapped.get_nlp_info().expect("init ok");
        assert_eq!(info.n, 1);
        assert_eq!(info.m, 2, "nonlinear rows are not dropped by presolve");

        // FBBT must have reached the infeasible row.
        let rpt = wrapped.fbbt_report().expect("fbbt ran");
        // No Phase-0 elimination is in force here, so the witness cannot be an
        // artifact — it is certified, and carries the witnessing row.
        assert_eq!(
            wrapped.certified_infeasible(),
            Some(InfeasibilityProof::IntervalArithmetic { witness: 1 }),
            "an FBBT witness on an un-clamped box must be certified"
        );
        assert_eq!(
            rpt.infeasibility_witness,
            Some(1),
            "row 1 (`x = 5`) is the infeasibility witness"
        );

        // ...but its corrupted tightening is discarded: the reduced box is
        // the original [0, 1], not the partial [0.3, 0.7] FBBT left behind.
        let mut x_l = vec![0.0; info.n as usize];
        let mut x_u = vec![0.0; info.n as usize];
        let mut g_l = vec![0.0; info.m as usize];
        let mut g_u = vec![0.0; info.m as usize];
        assert!(wrapped.get_bounds_info(BoundsInfo {
            x_l: &mut x_l,
            x_u: &mut x_u,
            g_l: &mut g_l,
            g_u: &mut g_u,
        }));
        assert_eq!(
            (x_l[0], x_u[0]),
            (0.0, 1.0),
            "FBBT's undefined-on-infeasibility bounds must not reach the IPM (H12)"
        );
    }

    /// 3-variable, 3-row TNLP that pits a Phase-0 auxiliary clamp against a
    /// kept nonlinear row. Rows 0/1 are a square LINEAR equality block
    /// (`x0 + x1 = 3`, `x0 - x1 = 1` → `(x0, x1) = (2, 1)`) that Phase 0
    /// eliminates, clamping `x0 = 2`. Row 2 is tagged `NonLinear` (so it is
    /// FBBT-handled, never turned into a linear row Phase 1 would catch) and
    /// reads `x0 + x2 = 20` with `x2 ∈ [0, 1]`. Over the *aux-clamped* box
    /// `x0 ∈ [2, 2]` FBBT sees `x0 + x2 ∈ [2, 3]`, disjoint from the required
    /// `20`, and witnesses infeasibility on row 2 — an infeasibility that
    /// exists ONLY because the aux clamp pinned `x0`. On the un-clamped box
    /// `x0` is free and FBBT tightens it to `[19, 20]` with no witness.
    struct AuxClampBreaksKeptNonlinearRow;
    impl TNLP for AuxClampBreaksKeptNonlinearRow {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 3,
                m: 3,
                nnz_jac_g: 6,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = -1e19;
            b.x_u[0] = 1e19;
            b.x_l[1] = -1e19;
            b.x_u[1] = 1e19;
            b.x_l[2] = 0.0;
            b.x_u[2] = 1.0;
            b.g_l[0] = 3.0;
            b.g_u[0] = 3.0;
            b.g_l[1] = 1.0;
            b.g_u[1] = 1.0;
            b.g_l[2] = 20.0;
            b.g_u[2] = 20.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = 0.0;
                sp.x[1] = 0.0;
                sp.x[2] = 0.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for v in g.iter_mut() {
                *v = 0.0;
            }
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] + x[1];
            g[1] = x[0] - x[1];
            g[2] = x[0] + x[2];
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 0, 1, 1, 2, 2]);
                    jcol.copy_from_slice(&[0, 1, 0, 1, 0, 2]);
                }
                SparsityRequest::Values { values } => {
                    values.copy_from_slice(&[1.0, 1.0, 1.0, -1.0, 1.0, 1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::Linear;
            types[1] = Linearity::Linear;
            // Tagged NonLinear so it is FBBT-handled, not a Phase-1 linear row.
            types[2] = Linearity::NonLinear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// Supplies the FBBT tape `x0 + x2` for row 2 only (rows 0/1 are linear
    /// and need no tape).
    struct AuxBreakProvider;
    impl ExpressionProvider for AuxBreakProvider {
        fn constraint_expression(
            &self,
            i: usize,
        ) -> Option<pounce_nlp::expression_provider::FbbtTape> {
            use pounce_nlp::expression_provider::{FbbtOp, FbbtTape};
            if i == 2 {
                Some(FbbtTape {
                    ops: vec![FbbtOp::Var(0), FbbtOp::Var(2), FbbtOp::Add(0, 1)],
                })
            } else {
                None
            }
        }
    }

    /// Regression for F6 (follow-up to H12). An FBBT infeasibility witnessed
    /// while a Phase-0 auxiliary elimination is in force must roll back that
    /// elimination and re-run FBBT on the un-clamped box — mirroring the
    /// Phase-1 aux rollback (#53). Restoring only the pre-FBBT box leaves the
    /// aux clamps in place, handing the IPM a reduced problem whose
    /// infeasibility is a presolve artifact (here `x0` pinned to 2 makes the
    /// kept row `x0 + x2 = 20` infeasible). After the rollback the dropped
    /// rows are restored, FBBT finds no infeasibility on the free box, and the
    /// solver sees the full original problem.
    #[test]
    fn fbbt_infeasibility_with_aux_clamp_rolls_back_phase0() {
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(AuxClampBreaksKeptNonlinearRow));
        let provider: Rc<RefCell<dyn ExpressionProvider>> = Rc::new(RefCell::new(AuxBreakProvider));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            bound_tightening: true,
            fbbt: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped =
            PresolveTnlp::with_expression_provider(Rc::clone(&inner), Rc::clone(&provider), opts);
        let info = wrapped.get_nlp_info().expect("init ok");

        // Confirm the scenario actually fired: Phase 0 eliminated the 2×2
        // linear block (clamping x0, x1). Without this the test would pass
        // vacuously.
        assert_eq!(
            wrapped.auxiliary_diagnostics().vars_eliminated,
            2,
            "Phase 0 must have clamped the linear block — otherwise the \
             rollback path is never exercised"
        );

        // The headline fix: the elimination is rolled back, so every dropped
        // row is restored. Without the rollback Phase 0's two dropped rows
        // would leave `m == 1`.
        assert_eq!(
            info.m, 3,
            "FBBT infeasibility under an aux clamp must roll Phase 0 back, \
             restoring the dropped rows (got m={})",
            info.m
        );

        // FBBT, re-run on the un-clamped box, finds no infeasibility: the
        // witness was an artifact of the aux clamp.
        let rpt = wrapped.fbbt_report().expect("fbbt ran");
        assert!(
            rpt.infeasibility_witness.is_none(),
            "FBBT must not witness infeasibility on the un-clamped box, got {:?}",
            rpt.infeasibility_witness
        );

        // The soundness property for the certification channel, on the one
        // scenario that can break it: this model *is feasible*, and the only
        // infeasibility ever seen was manufactured by Phase 0's own clamp.
        // Certifying it would make POUNCE report "proved infeasible" for a
        // problem that has solutions — strictly worse than the numerical
        // verdict this replaces. The proof must be withheld.
        assert_eq!(
            wrapped.certified_infeasible(),
            None,
            "a presolve-manufactured infeasibility must never be certified — \
             this model is feasible"
        );

        // x0 is no longer pinned to the aux value 2; FBBT over `x0 + x2 = 20`
        // with `x2 ∈ [0, 1]` tightens it to `[19, 20]`. Bounds stay valid.
        let mut x_l = vec![0.0; info.n as usize];
        let mut x_u = vec![0.0; info.n as usize];
        let mut g_l = vec![0.0; info.m as usize];
        let mut g_u = vec![0.0; info.m as usize];
        assert!(wrapped.get_bounds_info(BoundsInfo {
            x_l: &mut x_l,
            x_u: &mut x_u,
            g_l: &mut g_l,
            g_u: &mut g_u,
        }));
        for i in 0..(info.n as usize) {
            assert!(
                x_l[i] <= x_u[i] + 1e-12,
                "bounds handed to IPM must be valid: x_l[{i}]={} > x_u[{i}]={}",
                x_l[i],
                x_u[i]
            );
        }
        assert!(
            x_l[0] > 2.0 + 1e-6,
            "x0 must no longer be clamped to the aux value 2 (got x_l[0]={})",
            x_l[0]
        );
    }

    /// 2-variable, 2-row TNLP whose ORIGINAL problem is FEASIBLE (at
    /// `(x, y) = (2, 1)`) but where Phase-0's nonlinear block solve picks the
    /// WRONG root of a multi-root block. Row 0 is the 1×1 nonlinear block
    /// `x² = 4` (roots ±2); the starting point `x = -1` drives the damped
    /// Newton to the root `x = -2`. Row 1 is the kept row `x + y = 3` with
    /// `y ∈ [0, 2]` — tagged `NonLinear` so it is FBBT-handled (a `Linear`
    /// tag would let Phase 1 catch the clamp first and exercise the Phase-1
    /// rollback instead of the F6 FBBT rollback under test). Under the clamp
    /// `x = -2` the kept row needs `y = 5 > 2`: infeasible. At the OTHER root
    /// `x = +2` it needs `y = 1 ∈ [0, 2]`: feasible. So FBBT's witness is
    /// purely a Phase-0 artifact of the wrong root, on a feasible original.
    ///
    /// Phase-0 mechanics that make the scenario fire: DM matches row 0 to
    /// `x` (its only column) and row 1 to `y`; BTF solves the block
    /// `{row 0, x}` first (Newton from the probe `x = -1` → `-2`, inside
    /// `x ∈ [-10, 10]`, residual 0, accepted and clamped) and then rejects
    /// the block `{row 1, y}` as `OutOfBounds` (`y = 3 − (−2) = 5 ∉ [0, 2]`),
    /// which keeps row 1 visible to FBBT.
    struct WrongRootClampBreaksFeasibleOriginal;
    impl TNLP for WrongRootClampBreaksFeasibleOriginal {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 2,
                nnz_jac_g: 3,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = -10.0;
            b.x_u[0] = 10.0;
            b.x_l[1] = 0.0;
            b.x_u[1] = 2.0;
            b.g_l[0] = 4.0;
            b.g_u[0] = 4.0;
            b.g_l[1] = 3.0;
            b.g_u[1] = 3.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                // The probe seeds the block Newton: from -1 it converges to
                // the wrong root -2 (a start of +1 would find +2 and the
                // scenario would not fire).
                sp.x[0] = -1.0;
                sp.x[1] = 0.0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for v in g.iter_mut() {
                *v = 0.0;
            }
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] * x[0];
            g[1] = x[0] + x[1];
            true
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 1, 1]);
                    jcol.copy_from_slice(&[0, 0, 1]);
                }
                SparsityRequest::Values { values } => {
                    // ∂(x²)/∂x = 2x — genuinely point-dependent so the block
                    // Newton sees the true (sign-carrying) derivative at each
                    // iterate. Values requests always carry `Some(x)` here
                    // (probe fetch and Phase-0 callback both pass it).
                    let x0 = x.map(|x| x[0]).unwrap_or(-1.0);
                    values.copy_from_slice(&[2.0 * x0, 1.0, 1.0]);
                }
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types[0] = Linearity::NonLinear;
            // Tagged NonLinear so it is FBBT-handled, not a Phase-1 linear row.
            types[1] = Linearity::NonLinear;
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    /// Supplies the FBBT tape `x + y` for the kept row 1 only (row 0 is the
    /// eliminated block; it needs no tape for this scenario).
    struct WrongRootKeptRowProvider;
    impl ExpressionProvider for WrongRootKeptRowProvider {
        fn constraint_expression(
            &self,
            i: usize,
        ) -> Option<pounce_nlp::expression_provider::FbbtTape> {
            use pounce_nlp::expression_provider::{FbbtOp, FbbtTape};
            if i == 1 {
                Some(FbbtTape {
                    ops: vec![FbbtOp::Var(0), FbbtOp::Var(1), FbbtOp::Add(0, 1)],
                })
            } else {
                None
            }
        }
    }

    /// F6's headline failure mode, end to end: a FEASIBLE original problem
    /// must not be declared infeasible because Phase 0 clamped a multi-root
    /// nonlinear block at the wrong root. The companion test above
    /// (`fbbt_infeasibility_with_aux_clamp_rolls_back_phase0`) exercises the
    /// rollback *mechanism* but its original model is itself infeasible; here
    /// the original is feasible at `(2, 1)`, Newton lands on `x = -2`
    /// (breaking the kept row `x + y = 3`, `y ∈ [0, 2]`), FBBT witnesses the
    /// artifact infeasibility, and the rollback must rescue the solve: full
    /// row set restored, clamp gone, FBBT witness-free on the re-run — so the
    /// IPM receives the feasible full problem instead of a wrong "infeasible"
    /// verdict.
    #[test]
    fn fbbt_rollback_rescues_feasible_original_from_wrong_root_clamp() {
        let inner: Rc<RefCell<dyn TNLP>> =
            Rc::new(RefCell::new(WrongRootClampBreaksFeasibleOriginal));
        let provider: Rc<RefCell<dyn ExpressionProvider>> =
            Rc::new(RefCell::new(WrongRootKeptRowProvider));
        let opts = PresolveOptions {
            enabled: true,
            auxiliary: true,
            auxiliary_coupling: AuxiliaryCouplingPolicy::Safe,
            bound_tightening: true,
            fbbt: true,
            ..PresolveOptions::defaults()
        };
        let mut wrapped =
            PresolveTnlp::with_expression_provider(Rc::clone(&inner), Rc::clone(&provider), opts);
        let info = wrapped.get_nlp_info().expect("init ok");

        // Pre-rollback, the scenario must actually have fired: Phase 0
        // eliminated the 1×1 block, clamping x at a Newton root. Without
        // this the test passes vacuously (diagnostics survive the rollback).
        assert_eq!(
            wrapped.auxiliary_diagnostics().vars_eliminated,
            1,
            "Phase 0 must have solved and clamped the x² = 4 block — \
             otherwise the wrong-root path is never exercised"
        );

        // The rollback restored the full model. `m == 2` proves both halves:
        // (a) FBBT witnessed the artifact infeasibility (no witness → no
        // rollback → row 0 stays dropped → m == 1), and (b) that witness can
        // only arise from the wrong root (-2): at +2 the kept row is
        // satisfiable with y = 1 ∈ [0, 2] and FBBT finds nothing.
        assert_eq!(
            info.m, 2,
            "wrong-root clamp on a feasible original must trigger the F6 \
             rollback and restore the dropped block row (got m={})",
            info.m
        );

        // The re-run on the un-clamped box is witness-free: the original is
        // feasible, so nothing survives the rollback.
        let rpt = wrapped.fbbt_report().expect("fbbt ran");
        assert!(
            rpt.infeasibility_witness.is_none(),
            "feasible original: FBBT must not witness infeasibility after \
             the rollback, got {:?}",
            rpt.infeasibility_witness
        );

        // The box handed to the IPM is valid, the wrong-root clamp is gone,
        // and the feasible point (x, y) = (2, 1) is still inside it.
        let mut x_l = vec![0.0; info.n as usize];
        let mut x_u = vec![0.0; info.n as usize];
        let mut g_l = vec![0.0; info.m as usize];
        let mut g_u = vec![0.0; info.m as usize];
        assert!(wrapped.get_bounds_info(BoundsInfo {
            x_l: &mut x_l,
            x_u: &mut x_u,
            g_l: &mut g_l,
            g_u: &mut g_u,
        }));
        for i in 0..(info.n as usize) {
            assert!(
                x_l[i] <= x_u[i] + 1e-12,
                "bounds handed to IPM must be valid: x_l[{i}]={} > x_u[{i}]={}",
                x_l[i],
                x_u[i]
            );
        }
        assert!(
            x_u[0] - x_l[0] > 1e-6,
            "x must no longer be clamped (got [{}, {}])",
            x_l[0],
            x_u[0]
        );
        assert!(
            x_l[0] <= 2.0 + 1e-9 && 2.0 <= x_u[0] + 1e-9,
            "the feasible root x = 2 must survive in the IPM's box, got \
             [{}, {}]",
            x_l[0],
            x_u[0]
        );
        assert!(
            x_l[1] <= 1.0 + 1e-9 && 1.0 <= x_u[1] + 1e-9,
            "the feasible y = 1 must survive in the IPM's box, got [{}, {}]",
            x_l[1],
            x_u[1]
        );
    }

    /// L46: under a row drop, a genuinely per-constraint metadata vector
    /// (length `m_in`) must be subset to the kept rows, and the inverse
    /// `expand` must restore it (dropped rows defaulted). This is the
    /// canonical, contract-conforming case — the only one any real inner
    /// TNLP produces (e.g. `nl_reader`'s per-row `con_names`).
    #[test]
    fn con_metadata_per_row_vector_round_trips_under_row_drop() {
        let m_in = 3;
        let rows_kept = vec![0usize, 2]; // row 1 dropped by presolve
        let mut inner = MetaData::default();
        inner.strings.insert(
            "names".to_string(),
            vec!["c0".to_string(), "c1".to_string(), "c2".to_string()],
        );
        inner.integers.insert("flags".to_string(), vec![10, 11, 12]);
        inner
            .numerics
            .insert("weights".to_string(), vec![1.0, 2.0, 3.0]);

        let reduced = project_con_metadata(&inner, &rows_kept, m_in);
        assert_eq!(
            reduced.strings["names"],
            vec!["c0".to_string(), "c2".to_string()]
        );
        assert_eq!(reduced.integers["flags"], vec![10, 12]);
        assert_eq!(reduced.numerics["weights"], vec![1.0, 3.0]);

        // finalize_metadata's inverse: reduced (m_out=2) back to m_in=3.
        let restored = expand_con_metadata(&reduced, &rows_kept, m_in);
        assert_eq!(
            restored.strings["names"],
            vec!["c0".to_string(), String::new(), "c2".to_string()]
        );
        assert_eq!(restored.integers["flags"], vec![10, 0, 12]);
        assert_eq!(restored.numerics["weights"], vec![1.0, 0.0, 3.0]);
    }

    /// L46 hazard, characterized: a *global* metadata vector that happens
    /// to have length `m_in` is indistinguishable from a per-row vector
    /// (the untyped `MetaData` API carries no arity tag), so projection
    /// silently subsets it. This pins the known-imperfect behavior so a
    /// future arity-aware fix has a regression anchor; the value pounce
    /// ships today is unaffected because no real inner TNLP emits a
    /// coincidentally-`m_in`-length global constraint vector.
    #[test]
    fn con_metadata_length_heuristic_misfires_on_coincidental_global() {
        let m_in = 3;
        let rows_kept = vec![0usize, 2]; // a row is dropped, so projection is non-identity
        let mut inner = MetaData::default();
        // Semantically global (e.g. a 3-entry config tuple) that just
        // happens to be length 3 == m_in.
        inner
            .integers
            .insert("global_triple".to_string(), vec![100, 200, 300]);

        let reduced = project_con_metadata(&inner, &rows_kept, m_in);
        // Wrongly treated as per-row and subset to [100, 300]; the
        // faithful value would be the untouched [100, 200, 300].
        assert_eq!(
            reduced.integers["global_triple"],
            vec![100, 300],
            "documents the L46 length-heuristic misfire on a coincidental global vector"
        );

        // A global vector whose length differs from m_in is correctly
        // passed through untouched (the heuristic only misfires on the
        // exact-length coincidence).
        let mut inner2 = MetaData::default();
        inner2
            .numerics
            .insert("two_globals".to_string(), vec![1.5, 2.5]);
        let reduced2 = project_con_metadata(&inner2, &rows_kept, m_in);
        assert_eq!(reduced2.numerics["two_globals"], vec![1.5, 2.5]);
    }

    /// gh#391 scratch. One variable in `[0, 1]`, `m` rows evaluated as
    /// `coeff[i] * x + offset[i]` against per-row declared bounds — enough to
    /// drive `witness_refutes_infeasibility` directly and compare the two
    /// [`WitnessRule`] forms on the same numbers.
    struct WitnessProbeRows {
        coeff: Vec<Number>,
        offset: Vec<Number>,
        g_l: Vec<Number>,
        g_u: Vec<Number>,
        x0: Number,
    }

    impl TNLP for WitnessProbeRows {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: self.coeff.len() as Index,
                nnz_jac_g: self.coeff.len() as Index,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = 0.0;
            b.x_u[0] = 1.0;
            b.g_l.copy_from_slice(&self.g_l);
            b.g_u.copy_from_slice(&self.g_u);
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            if sp.init_x {
                sp.x[0] = self.x0;
            }
            true
        }
        fn eval_f(&mut self, _x: &[Number], _new_x: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            for ((out, &c), &o) in g.iter_mut().zip(&self.coeff).zip(&self.offset) {
                *out = c * x[0] + o;
            }
            true
        }
        fn eval_jac_g(
            &mut self,
            _x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    for (i, r) in irow.iter_mut().enumerate() {
                        *r = i as Index;
                    }
                    jcol.fill(0);
                }
                SparsityRequest::Values { values } => values.copy_from_slice(&self.coeff),
            }
            true
        }
        fn get_constraints_linearity(&mut self, types: &mut [Linearity]) -> bool {
            types.fill(Linearity::Linear);
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    fn refutes(model: WitnessProbeRows, rule: WitnessRule) -> bool {
        let (g_l, g_u) = (model.g_l.clone(), model.g_u.clone());
        let m = g_l.len();
        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(model));
        witness_refutes_infeasibility(&inner, 1, m, &[0.0], &[1.0], &g_l, &g_u, 1e-8, rule)
    }

    /// The gh#391 mechanism in one assertion pair, on identical numbers.
    ///
    /// `s*x == 0.2*s` with `s*x == 0.8*s` over `x ∈ [0, 1]` is empty at every
    /// `s > 0` — the crossing is `0.6` regardless. At `s = 1e-12` the clamped
    /// accepting form reinstates the absolute floor (`tol * max(scale, 1)` =
    /// `1e-8`), every point of the box "satisfies" both rows, and the witness
    /// withdraws a scale-free proof. Measured against the row's declared
    /// magnitude with no clamp, the same point violates by `0.6` relative and
    /// refutes nothing.
    #[test]
    fn down_scaled_rows_only_defeat_the_clamped_witness_form() {
        let model = |_s: Number| WitnessProbeRows {
            coeff: vec![1e-12, 1e-12],
            offset: vec![0.0, 0.0],
            g_l: vec![0.2e-12, 0.8e-12],
            g_u: vec![0.2e-12, 0.8e-12],
            x0: 0.5,
        };
        assert!(
            refutes(model(1e-12), WitnessRule::SolverAcceptance),
            "the clamped form accepts every point of the box at this row scale \
             — this is the behavior gh#391 is about, pinned so the contrast \
             below stays meaningful"
        );
        assert!(
            !refutes(model(1e-12), WitnessRule::DeclaredRowRelative),
            "0.6 of the row's own declared magnitude is a violation at every \
             scale; the strict rule must not withdraw the proof"
        );
    }

    /// The same rows at unit scale: both forms agree, which is the point —
    /// the strict rule changes nothing where the clamp was never active.
    #[test]
    fn unit_scale_rows_refute_under_neither_form() {
        let model = || WitnessProbeRows {
            coeff: vec![1.0, 1.0],
            offset: vec![0.0, 0.0],
            g_l: vec![0.2, 0.8],
            g_u: vec![0.2, 0.8],
            x0: 0.5,
        };
        assert!(!refutes(model(), WitnessRule::SolverAcceptance));
        assert!(!refutes(model(), WitnessRule::DeclaredRowRelative));
    }

    /// The hazard the strict rule has to handle explicitly: a *homogeneous* row
    /// (`g_l = g_u = 0`) has no declared magnitude, so a pure relative test
    /// would be unsatisfiable by construction — the violation *is* the scale,
    /// and any float-noise residual reads as a full-magnitude violation. Such a
    /// row keeps the absolute floor, so a genuinely feasible point still
    /// refutes and the proof is still withheld. Fail closed.
    #[test]
    fn homogeneous_row_keeps_the_absolute_floor_under_the_strict_rule() {
        // `x - 0.5 == 0` evaluated at the modeller's `x0 = 0.5`, off by an ulp
        // of noise. Feasible; the witness must say so under both rules.
        let noisy = || WitnessProbeRows {
            coeff: vec![1.0],
            offset: vec![-0.5 + 1e-16],
            g_l: vec![0.0],
            g_u: vec![0.0],
            x0: 0.5,
        };
        assert!(refutes(noisy(), WitnessRule::SolverAcceptance));
        assert!(
            refutes(noisy(), WitnessRule::DeclaredRowRelative),
            "a homogeneous row has nothing to be relative to; float noise on it \
             must not be promoted to a violation"
        );

        // A real violation on a homogeneous row still refutes nothing: the
        // fallback is the absolute floor, not blanket acceptance.
        let violated = || WitnessProbeRows {
            coeff: vec![1.0],
            offset: vec![0.6],
            g_l: vec![0.0],
            g_u: vec![0.0],
            x0: 0.5,
        };
        assert!(!refutes(violated(), WitnessRule::SolverAcceptance));
        assert!(!refutes(violated(), WitnessRule::DeclaredRowRelative));
    }

    /// The strict rule is relative, not absolute: on a row declared near `1e30`
    /// a residual of `1e5` is fifteen relative digits of agreement and still
    /// refutes. This is the extreme-coefficient class the witness gate exists
    /// for, and it must survive the rule change.
    #[test]
    fn strict_rule_still_refutes_at_extreme_row_magnitude() {
        let model = || WitnessProbeRows {
            coeff: vec![2e30],
            offset: vec![1e5],
            g_l: vec![1e30],
            g_u: vec![1e30],
            x0: 0.5,
        };
        assert!(refutes(model(), WitnessRule::SolverAcceptance));
        assert!(refutes(model(), WitnessRule::DeclaredRowRelative));
    }

    /// End-to-end through the wrapper: the rule is a property of the call site,
    /// so a default wrapper (one that would be solved through) keeps the #380
    /// behavior at sub-tolerance row scale, and only a caller that has
    /// established the solve cannot run gets the proof.
    #[test]
    fn probing_without_a_solve_certifies_where_the_default_wrapper_withholds() {
        let opts = PresolveOptions {
            enabled: true,
            bound_tightening: true,
            auxiliary: false,
            ..PresolveOptions::defaults()
        };
        let model = || WitnessProbeRows {
            coeff: vec![1e-12, 1e-12],
            offset: vec![0.0, 0.0],
            g_l: vec![0.2e-12, 0.8e-12],
            g_u: vec![0.2e-12, 0.8e-12],
            x0: 0.5,
        };

        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(model()));
        let mut solved_through = PresolveTnlp::new(inner, opts);
        solved_through.get_nlp_info().expect("init ok");
        assert!(
            solved_through.tighten_report().infeasible,
            "bound propagation sees the contradiction at every scale"
        );
        assert_eq!(
            solved_through.certified_infeasible(),
            None,
            "#380: a wrapper that will be solved through must not claim a proof \
             the solver's own acceptance test would contradict"
        );

        let inner: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(model()));
        let mut probe = PresolveTnlp::new(inner, opts).probing_without_a_solve();
        probe.get_nlp_info().expect("init ok");
        assert_eq!(
            probe.certified_infeasible(),
            Some(InfeasibilityProof::BoundPropagation),
            "gh#391: with no solve to contradict, the scale-free crossing is a \
             proof at every row scale"
        );
    }
}