pounce_nlp/

orig_ipopt_nlp.rs

//! `OrigIpoptNlp` — concrete `IpoptNlp` impl that wraps a [`TNLPAdapter`]
//! and an [`NlpScaling`] object. Port of
//! `Algorithm/IpOrigIpoptNLP.{hpp,cpp}` (Ipopt 3.14.19).
//!
//! # Design
//!
//! Upstream `OrigIpoptNLP` does four things:
//!
//! 1. Holds the equality / inequality-separated bound vectors
//!    (`x_L`, `x_U`, `d_L`, `d_U`) and the four expansion matrices
//!    (`Px_L`, `Px_U`, `Pd_L`, `Pd_U`).
//! 2. Routes `f / grad_f / c / d / jac_c / jac_d / h` evaluations down
//!    to the user's `NLP` (i.e. our `TNLP` via `TNLPAdapter`), splitting
//!    constraints into c/d and applying scaling.
//! 3. Caches each result keyed on the input vector tag (`CachedResults`
//!    upstream → [`pounce_common::cached::Cache`] here).
//! 4. Counts evaluations and forwards the unscaled solution to
//!    `TNLP::finalize_solution`.
//!
//! # Trait location
//!
//! The `Nlp` / `IpoptNlp` traits live in [`crate::ipopt_nlp`] and are
//! re-exported from `pounce_algorithm::ipopt_nlp` so the algorithm-side
//! code can keep its existing `crate::ipopt_nlp::IpoptNlp` import path.
//! We moved the traits down to `pounce-nlp` because the concrete impl
//! has to live alongside `TNLPAdapter` (its private dependency) and
//! `pounce-algorithm` already depends on `pounce-nlp` (the reverse
//! would cycle).
//!
//! # Phase scope
//!
//! Implemented for v1.0:
//! * `f / grad_f / c / d / jac_c / jac_d / h` with one-dependency
//!   tag-keyed caches.
//! * Bound vectors and 0/1 expansion matrices for `(x_L, x_U, d_L, d_U)`.
//! * Starting point retrieval + initial multiplier handling.
//! * Per-eval counters (used to populate `SolveStatistics`).
//! * `finalize_solution` plumbing.
//!
//! Deferred (phase numbers from
//! `we-are-going-to-polished-simon.md`):
//! * Phase 8: L-BFGS / SR1 quasi-Newton path
//!   (`hessian_approximation = limited-memory`). The `eval_h` path is
//!   wired here, but the `LowRankUpdateSymMatrix` h_space construction
//!   in `IpOrigIpoptNLP.cpp:251-278` is not.
//! * Phase 10: adaptive-mu's `objective_depends_on_mu` /
//!   `f(x, mu)` overload (CG-penalty objective).
//! * Bound relaxation (`bound_relax_factor`) and `honor_original_bounds`
//!   projection — these need `OptionsList` plumbing that lands later.
//! * `check_derivatives_for_naninf` — needs the journalist's NaN
//!   reporting, deferred with the option.
//! * Full `NLPScalingObject` integration (currently only `obj_scaling`
//!   is used; `apply_vector_scaling_*`, `apply_jac_*_scaling`,
//!   `apply_hessian_scaling` live behind future scaling-object API).
//! * Fixed-variable removal (`x_l == x_u`) — `TNLPAdapter` keeps fixed
//!   variables in `x_var` for now; the upstream
//!   `fixed_variable_treatment` knob lands when the option machinery
//!   does.

use crate::constant_derivatives::{
    ConstantDerivatives, DerivativeProof, DerivativeProofs, subsystem_proof,
};
use crate::ipopt_nlp::{IpoptNlp, Nlp, SplitNames};
use crate::tnlp::{IDX_NAMES, MetaData, NlpInfo, ScalingRequest, SparsityRequest, StartingPoint};
use crate::tnlp_adapter::{BoundClassification, TNLPAdapter};
use pounce_common::cached::Cache;
use pounce_common::timing::TimingStatistics;
use pounce_common::types::{Index, Number};
use pounce_linalg::{
    DenseVector, DenseVectorSpace, ExpansionMatrix, ExpansionMatrixSpace, GenTMatrix,
    GenTMatrixSpace, Matrix, SymMatrix, SymTMatrix, SymTMatrixSpace, Vector,
};
use std::cell::{Cell, RefCell};
use std::rc::Rc;

/// Opaque scaling-object handle. `OrigIpoptNlp` only consults this for
/// the *initial* objective scaling; the full per-row constraint /
/// Jacobian / Hessian scaling lives directly on `OrigIpoptNlp` (see
/// the `obj_scale_factor` / `c_scale` / `d_scale` fields and the
/// `determine_scaling_from_starting_point` method) so that the runtime
/// can compute gradient-based scaling without an upcall.
///
/// The trait is intentionally minimal and local to `pounce-nlp`: the
/// gradient-based scaling arithmetic lives on `OrigIpoptNlp` itself
/// (see `determine_scaling_from_starting_point`), so there is no
/// algorithm-layer scaling strategy object.
pub trait NlpScaling {
    /// Optional user-supplied multiplier on the objective scaling
    /// factor. Mirrors upstream's `obj_scaling_factor` option (default
    /// 1.0). Combined with the gradient-based factor in
    /// `OrigIpoptNlp::determine_scaling_from_starting_point`.
    fn obj_scaling(&self) -> Number {
        1.0
    }
}

/// No-op scaling — every factor is 1.0. Default for unit tests and
/// callers that have not configured a scaling strategy.
#[derive(Debug, Default, Clone, Copy)]
pub struct NoScaling;
impl NlpScaling for NoScaling {}

/// Constant objective scaling: carries the user's `obj_scaling_factor`
/// option value into [`OrigIpoptNlp`]. A negative factor flips the
/// optimization direction (the IPM minimizes `factor·f`, i.e.
/// maximizes `f`), matching upstream Ipopt's documented semantics.
#[derive(Debug, Clone, Copy)]
pub struct ConstObjScaling(pub Number);
impl NlpScaling for ConstObjScaling {
    fn obj_scaling(&self) -> Number {
        self.0
    }
}

/// Selector for [`OrigIpoptNlp::determine_scaling_from_starting_point`].
/// Mirrors upstream's `nlp_scaling_method` option.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ScalingMethod {
    /// No automatic scaling beyond the constant `obj_scaling_factor`.
    None,
    /// Gradient-based per `Algorithm/IpGradientScaling.cpp`. Default.
    GradientBased,
    /// User-supplied scaling via [`crate::tnlp::TNLP::get_scaling_parameters`].
    /// Port of upstream's `nlp_scaling_method=user-scaling`. The TNLP
    /// fills `obj_scaling` and the per-constraint `g_scaling`;
    /// `OrigIpoptNlp` models no variable-side rescale (the issue #61
    /// design), so a non-trivial per-variable `x_scaling` request is
    /// **refused** — see [`OrigIpoptNlp::user_x_scaling_rejected`].
    UserScaling,
}

/// The objective scale `df` gradient-based scaling picks from
/// `‖∇f(x₀)‖_∞`.
///
/// Free rather than inlined into [`OrigIpoptNlp::determine_scaling_from_starting_point`]
/// so a diagnostic can report the number the solver *will* use without
/// re-deriving it: `pounce check-x0` calls this, and a copy there would be
/// free to drift from the copy that runs (gh #703).
///
/// `obj_target_gradient > 0` (`nlp_scaling_obj_target_gradient`) pins the
/// scaled gradient's ∞-norm to that value and overrides both the
/// `max_gradient` cutoff and the 1.0 clamp; otherwise the objective is
/// scaled only when its gradient exceeds the cutoff. The result is floored
/// at `min_value`.
pub fn gradient_obj_scale(
    max_grad_f: Number,
    max_gradient: Number,
    min_value: Number,
    obj_target_gradient: Number,
) -> Number {
    let mut df = 1.0;
    if obj_target_gradient > 0.0 && max_grad_f > 0.0 {
        // Target overrides the cutoff (and the 1.0 clamp):
        // pin gradient ∞-norm to the requested value.
        df = obj_target_gradient / max_grad_f;
    } else if max_grad_f > max_gradient {
        df = max_gradient / max_grad_f;
    }
    if df < min_value {
        df = min_value;
    }
    df
}

/// The scale gradient-based scaling assigns one constraint row, given that
/// row's Jacobian ∞-norm at the starting point.
///
/// With `constr_target_gradient > 0` the user is asking for a *fixed*
/// gradient ∞-norm per row (overrides the cutoff and the 1.0 clamp).
/// Otherwise: scale only rows that exceed the cutoff, never amplify (clamp
/// at 1). Floored at `min_value`.
///
/// **A row whose Jacobian is entirely zero at the sampling point gets 1.0**
/// — `row_max` arrives as `f64::MIN_POSITIVE` (upstream's `dbl_min` seed),
/// the raw ratio overflows the clamp, and the row comes out unscaled. That
/// is the right answer for a row with no derivative at all, and the wrong
/// one for a row whose derivative merely *vanishes at x₀* — a
/// `½xᵀQx ≤ b` written about the origin and started from `x₀ = 0`, which
/// is how AMPL emits `qcqp1000-2c`. See `pounce check-x0`'s scaling
/// section and `dev-notes/quadratic-structure-exploitation.md` §8
/// (gh #703).
pub fn gradient_row_scale(
    row_max: Number,
    max_gradient: Number,
    min_value: Number,
    constr_target_gradient: Number,
) -> Number {
    let mut s = if constr_target_gradient > 0.0 {
        constr_target_gradient / row_max
    } else {
        let raw = max_gradient / row_max;
        if raw > 1.0 { 1.0 } else { raw }
    };
    if s < min_value {
        s = min_value;
    }
    s
}

/// Whether gradient-based scaling produces a scale vector at all for a
/// block of rows (the `c` equalities and the `d` inequalities are gated
/// separately). Unless some row exceeds the cutoff the whole block is left
/// unscaled, so a per-row scale below 1 elsewhere in the model does not
/// imply this block got one.
pub fn gradient_scaling_fires(
    row_max: &[Number],
    max_gradient: Number,
    constr_target_gradient: Number,
) -> bool {
    constr_target_gradient > 0.0 || row_max.iter().any(|&v| v > max_gradient)
}

/// Concrete `IpoptNlp` over a `TNLPAdapter`. Mirrors upstream
/// `Ipopt::OrigIpoptNLP`.
pub struct OrigIpoptNlp {
    /// Backing TNLP (and its bound classification).
    adapter: Rc<RefCell<TNLPAdapter>>,
    /// Constant objective-scaling multiplier supplied by the user
    /// (mirrors upstream's `obj_scaling_factor` option). The
    /// gradient-based factor is multiplied into [`obj_scale_factor`]
    /// after [`Self::determine_scaling_from_starting_point`] runs.
    scaling: Rc<dyn NlpScaling>,

    // ----- gradient-based scaling state (port of `IpGradientScaling.cpp`) -----
    /// Effective objective scaling factor `df_` (1.0 when no scaling).
    obj_scale_factor: Cell<Number>,
    /// The gradient-based factor alone (`df`), without the user's constant.
    /// See `IpoptNlp::computed_obj_scaling_factor`.
    computed_obj_scale: Cell<Number>,
    /// Per-row scaling for equality constraints (`dc_`). `None` ↔
    /// `IsValid(dc) == false` — i.e. row-max gradient is below the
    /// `nlp_scaling_max_gradient` cutoff so no scaling is applied.
    c_scale: RefCell<Option<Vec<Number>>>,
    /// Same as [`Self::c_scale`] but for inequality rows.
    d_scale: RefCell<Option<Vec<Number>>>,
    /// The *declared* compressed `d_L / d_U` — snapshotted (unscaled) at
    /// [`Self::relax_bounds`] time, before `bound_relax_factor` widens the
    /// live bounds and before any safe-slack adjustment moves them. The
    /// scale-relative feasibility measure keys row magnitudes off these: a
    /// relaxed zero bound reads as `~1e-8` on the live vector and would
    /// fabricate a magnitude for a row that has none. `None` until
    /// `relax_bounds` runs (direct drivers that skip it fall back to the
    /// live bounds).
    declared_d_l: RefCell<Option<Vec<Number>>>,
    declared_d_u: RefCell<Option<Vec<Number>>>,
    /// The *declared* compressed `x_L / x_U`, snapshotted (unrelaxed) at
    /// [`Self::relax_bounds`] time for the same reason as
    /// [`Self::declared_d_l`] — and used by
    /// [`Self::finalize_solution_x`] to project the final iterate back
    /// into the user's own box under `honor_original_bounds`.
    declared_x_l: RefCell<Option<Vec<Number>>>,
    declared_x_u: RefCell<Option<Vec<Number>>>,
    /// `honor_original_bounds` (upstream default `no`). When set, the
    /// final iterate handed to `TNLP::finalize_solution` is projected
    /// back into the declared bounds, undoing the `bound_relax_factor`
    /// widening. Registered but never read before gh#483's follow-up, so
    /// `honor_original_bounds=yes` was accepted and the reported solution
    /// still sat up to `min(bound_relax_factor·max(1,|b|), constr_viol_tol)`
    /// outside the box the user declared.
    honor_original_bounds: Cell<bool>,
    /// Set by [`Self::scale_user_supplied`] when the TNLP asked for
    /// **non-trivial** per-variable scaling (`use_x_scaling` with at
    /// least one factor `!= 1.0`). pounce models objective and
    /// constraint scaling only, so honoring the request is impossible
    /// and *ignoring* it silently hands back a differently-conditioned
    /// problem than the one the caller asked for (gh#483). The driver
    /// reads this via [`Self::user_x_scaling_rejected`] and fails the
    /// solve with `InvalidOption` instead.
    x_scaling_rejected: Cell<bool>,

    // ----- vector / matrix spaces (shared via Rc) -----
    x_space: Rc<DenseVectorSpace>,
    c_space: Rc<DenseVectorSpace>,
    d_space: Rc<DenseVectorSpace>,
    x_l_space: Rc<DenseVectorSpace>,
    x_u_space: Rc<DenseVectorSpace>,
    d_l_space: Rc<DenseVectorSpace>,
    d_u_space: Rc<DenseVectorSpace>,
    px_l_space: Rc<ExpansionMatrixSpace>,
    px_u_space: Rc<ExpansionMatrixSpace>,
    pd_l_space: Rc<ExpansionMatrixSpace>,
    pd_u_space: Rc<ExpansionMatrixSpace>,
    jac_c_space: Rc<GenTMatrixSpace>,
    jac_d_space: Rc<GenTMatrixSpace>,
    /// Hessian space; `None` when `eval_h` is not provided by the TNLP
    /// (the limited-memory quasi-Newton path lands in Phase 8).
    h_space: Option<Rc<SymTMatrixSpace>>,

    // ----- bound vectors (compressed-x sub-spaces) -----
    x_l: Rc<DenseVector>,
    x_u: Rc<DenseVector>,
    d_l: Rc<DenseVector>,
    d_u: Rc<DenseVector>,
    /// Constant equality right-hand side (upstream's `c_rhs`): for each
    /// equality row `i`, the bound `g_l[c_map[i]] == g_u[c_map[i]]`.
    /// Captured once at construction so [`Self::eval_c_internal`] forms the
    /// residual `g - c_rhs` without re-fetching all bounds (and the four
    /// full-size scratch allocations that requires) on every line-search
    /// trial. (Code review 2026-06 item M17.)
    c_rhs: Vec<Number>,
    /// Full TNLP start point fetched once for a warm start, then projected
    /// through the split x/y/z accessor calls without re-entering the TNLP.
    warm_start_snapshot: RefCell<Option<StartingPointSnapshot>>,

    // ----- expansion matrices (instances; spaces above) -----
    px_l: Rc<dyn Matrix>,
    px_u: Rc<dyn Matrix>,
    pd_l: Rc<dyn Matrix>,
    pd_u: Rc<dyn Matrix>,

    // ----- jacobian sparsity remap -----
    /// `jac_c_entry_in_g[k]` = position in the full TNLP jacobian's
    /// values array of the k-th equality-row entry.
    jac_c_entry_in_g: Vec<Index>,
    /// Same for inequality rows.
    jac_d_entry_in_g: Vec<Index>,
    /// Total nonzeros in the full (un-split) `eval_jac_g` triplet.
    nnz_jac_g_full: Index,

    // ----- hessian sparsity remap (fixed-var filtering) -----
    /// Total nonzeros the user's `eval_h` writes into. May exceed
    /// `h_space.nonzeros()` when fixed variables drop entries.
    nnz_h_lag_full: Index,
    /// `h_entry_in_full[k]` = position in the full TNLP hessian's
    /// values array of the k-th kept entry. Always has length
    /// `h_space.nonzeros()`; equals the identity `[0, 1, …, n-1]`
    /// when no fixed-var filtering dropped any entries.
    h_entry_in_full: Vec<Index>,

    // ----- caches (one entry; key = input vector tag) -----
    f_cache: RefCell<Cache<Number>>,
    grad_f_cache: RefCell<Cache<Rc<dyn Vector>>>,
    c_cache: RefCell<Cache<Rc<dyn Vector>>>,
    d_cache: RefCell<Cache<Rc<dyn Vector>>>,
    jac_c_cache: RefCell<Cache<Rc<dyn Matrix>>>,
    jac_d_cache: RefCell<Cache<Rc<dyn Matrix>>>,
    h_cache: RefCell<Cache<Rc<dyn SymMatrix>>>,
    /// Shared full-space buffers below the c/d split, so the dominant AD
    /// cost is paid once per iterate instead of twice. `eval_c`/`eval_d`
    /// both slice their rows out of one `eval_g` result (`full_g_cache`),
    /// and `eval_jac_c`/`eval_jac_d` both slice one `eval_jac_g` result
    /// (`full_jac_g_cache`). Keyed by the input vector's tag, like the
    /// per-subsystem caches; mirrors upstream's tagged `full_g_`/`jac_g_`
    /// buffers. (Code review 2026-06 item M16.)
    full_g_cache: RefCell<Cache<Rc<Vec<Number>>>>,
    full_jac_g_cache: RefCell<Cache<Rc<Vec<Number>>>>,

    // ----- evaluation counters -----
    f_evals: RefCell<Index>,
    grad_f_evals: RefCell<Index>,
    c_evals: RefCell<Index>,
    d_evals: RefCell<Index>,
    jac_c_evals: RefCell<Index>,
    jac_d_evals: RefCell<Index>,
    h_evals: RefCell<Index>,

    /// Cached `NlpInfo` (n, m, nnz_jac_g, nnz_h_lag, index_style) so we
    /// don't re-borrow the TNLP for dimension queries.
    info: NlpInfo,

    /// Shared per-subsystem timing accumulator. `None` until
    /// `IpoptApplication` installs the shared `Rc<TimingStatistics>` via
    /// [`Self::set_timing_stats`]; when `None`, all `eval_*` entry
    /// points skip the timing overhead.
    timing: RefCell<Option<Rc<TimingStatistics>>>,

    /// Which derivatives this solve reuses across iterates — the four
    /// upstream `*_constant` hints, resolved against what the model
    /// could prove (gh #588, Q6). Default: reuse nothing, which is the
    /// behaviour every path had before this field existed.
    ///
    /// A reused derivative is stored in its ordinary cache with an
    /// **empty** dependency list instead of the iterate's tag, so the
    /// lookup that already runs on every call finds it at every point.
    /// That is upstream's own mechanism (`IpOrigIpoptNLP.cpp` nulls the
    /// dependency); what pounce does differently is decide *when* it is
    /// allowed to (see [`crate::constant_derivatives`]).
    const_deriv: ConstantDerivatives,

    /// Which forms a *fixed* variable can actually reach — the guard on
    /// the `forget_variation` weakening in [`Self::derivative_proofs`].
    /// All-clear (and unread) when the model fixes nothing.
    fixed_reach: FixedVarReach,
}

/// Where `make_parameter`'s fixed variables are structurally reachable.
///
/// Fixing a variable can only *remove* variation from a derivative, so a
/// `Varying` proof does not survive it — that is why
/// [`DerivativeProof::forget_variation`] exists. But the weakening is
/// not free: `Varying` refuses a false user hint, and `Unknown` honours
/// it on trust (see the table in [`crate::constant_derivatives`]). Applied
/// to a form no fixed variable can reach, it would trade a refusal pounce
/// had earned for trust it has no basis for — on a model whose only sin
/// was fixing one unrelated variable. So it is applied per form, and this
/// records which forms qualify.
///
/// Both tests are read straight out of the sparsity the constructor
/// already queries for the `make_parameter` remap, so neither costs an
/// extra call into the user's TNLP.
#[derive(Clone, Debug, Default)]
struct FixedVarReach {
    /// `g_row_touches_fixed[i]` — row `i` of the full TNLP Jacobian has a
    /// declared nonzero in a fixed column. If `∂gᵢ/∂x_j` depends on a
    /// fixed `x_p` then `gᵢ` depends on `x_p`, so `∂gᵢ/∂x_p ≢ 0` and the
    /// entry is there: a row without one cannot change character.
    /// Indexed by the model's own row order, like `DerivativeProofs::jac`.
    g_row_touches_fixed: Vec<bool>,
    /// Whether a fixed variable is reachable in *second* order — set when
    /// the declared Hessian has an entry whose row or column is fixed, and
    /// (conservatively) whenever the TNLP declines to declare a Hessian at
    /// all, which leaves the question unanswerable.
    ///
    /// This is the test for both `∇f` and `∇²L`, because both are second
    /// order in the fixed variable: if `∂f/∂x_j` depends on a fixed `x_p`
    /// then `∂²f/∂x_j∂x_p ≢ 0`, and if `∂²L/∂x_j∂x_k` does then
    /// `∂²L/∂x_j∂x_p ≢ 0` — either way a Hessian entry in a fixed
    /// row/column is declared. A first-order-only appearance (`f` linear
    /// in `x_p`) changes no derivative w.r.t. the remaining variables.
    second_order_touches_fixed: bool,
}

#[derive(Clone)]
struct StartingPointSnapshot {
    x: Vec<Number>,
    z_l: Vec<Number>,
    z_u: Vec<Number>,
    lambda: Vec<Number>,
}

impl std::fmt::Debug for OrigIpoptNlp {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("OrigIpoptNlp")
            .field("info", &self.info)
            .field("f_evals", &*self.f_evals.borrow())
            .field("grad_f_evals", &*self.grad_f_evals.borrow())
            .field("c_evals", &*self.c_evals.borrow())
            .field("d_evals", &*self.d_evals.borrow())
            .field("jac_c_evals", &*self.jac_c_evals.borrow())
            .field("jac_d_evals", &*self.jac_d_evals.borrow())
            .field("h_evals", &*self.h_evals.borrow())
            .finish_non_exhaustive()
    }
}

impl OrigIpoptNlp {
    /// Construct an `OrigIpoptNlp` from a (already-classified) adapter
    /// and a scaling object. Mirrors
    /// `OrigIpoptNLP::OrigIpoptNLP` + `InitializeStructures`
    /// (`IpOrigIpoptNLP.cpp:22-457`) — the parts that don't need an
    /// `OptionsList` (those land with the option-machinery integration).
    pub fn new(
        adapter: Rc<RefCell<TNLPAdapter>>,
        scaling: Rc<dyn NlpScaling>,
    ) -> Result<Self, String> {
        // Snapshot dimensions / classification from the adapter.
        let (info, classification) = {
            let a = adapter.borrow();
            (*a.nlp_info(), a.classification().clone())
        };

        // ---- Vector spaces ----
        let n_x_var = classification.n_x_var();
        let x_space = DenseVectorSpace::new(n_x_var);
        let c_space = DenseVectorSpace::new(classification.n_c);
        let d_space = DenseVectorSpace::new(classification.n_d);
        let x_l_space = DenseVectorSpace::new(classification.n_x_l());
        let x_u_space = DenseVectorSpace::new(classification.n_x_u());
        let d_l_space = DenseVectorSpace::new(classification.n_d_l());
        let d_u_space = DenseVectorSpace::new(classification.n_d_u());

        // ---- Expansion matrix spaces (column-compressed → full x_var / d) ----
        let px_l_space =
            ExpansionMatrixSpace::new(n_x_var, classification.n_x_l(), &classification.x_l_map, 0);
        let px_u_space =
            ExpansionMatrixSpace::new(n_x_var, classification.n_x_u(), &classification.x_u_map, 0);
        let pd_l_space = ExpansionMatrixSpace::new(
            classification.n_d,
            classification.n_d_l(),
            &classification.d_l_map,
            0,
        );
        let pd_u_space = ExpansionMatrixSpace::new(
            classification.n_d,
            classification.n_d_u(),
            &classification.d_u_map,
            0,
        );
        let px_l: Rc<dyn Matrix> = Rc::new(ExpansionMatrix::new(Rc::clone(&px_l_space)));
        let px_u: Rc<dyn Matrix> = Rc::new(ExpansionMatrix::new(Rc::clone(&px_u_space)));
        let pd_l: Rc<dyn Matrix> = Rc::new(ExpansionMatrix::new(Rc::clone(&pd_l_space)));
        let pd_u: Rc<dyn Matrix> = Rc::new(ExpansionMatrix::new(Rc::clone(&pd_u_space)));

        // ---- Bound vectors. Pull the full `x_l/x_u/g_l/g_u` arrays from
        // the TNLP (the adapter discarded them after classification) and
        // pick out the entries pointed at by the `*_map` tables. -----
        let n_full_x = classification.n_full_x as usize;
        let n_full_g = classification.n_full_g as usize;
        let mut full_x_l = vec![0.0; n_full_x];
        let mut full_x_u = vec![0.0; n_full_x];
        let mut full_g_l = vec![0.0; n_full_g];
        let mut full_g_u = vec![0.0; n_full_g];
        {
            let a = adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            let ok = t.get_bounds_info(crate::tnlp::BoundsInfo {
                x_l: &mut full_x_l,
                x_u: &mut full_x_u,
                g_l: &mut full_g_l,
                g_u: &mut full_g_u,
            });
            if !ok {
                return Err("TNLP::get_bounds_info returned false on second call".into());
            }
        }

        let x_l = make_dense_from(&x_l_space, |i| {
            // x_l_map[i] is an index into x_var (== index into x_not_fixed_map).
            let var_idx = classification.x_l_map[i] as usize;
            let full_idx = classification.x_not_fixed_map[var_idx] as usize;
            full_x_l[full_idx]
        });
        let x_u = make_dense_from(&x_u_space, |i| {
            let var_idx = classification.x_u_map[i] as usize;
            let full_idx = classification.x_not_fixed_map[var_idx] as usize;
            full_x_u[full_idx]
        });
        let d_l = make_dense_from(&d_l_space, |i| {
            // d_l_map[i] is an index into d (== position in d_map).
            let d_idx = classification.d_l_map[i] as usize;
            let full_g_idx = classification.d_map[d_idx] as usize;
            full_g_l[full_g_idx]
        });
        let d_u = make_dense_from(&d_u_space, |i| {
            let d_idx = classification.d_u_map[i] as usize;
            let full_g_idx = classification.d_map[d_idx] as usize;
            full_g_u[full_g_idx]
        });

        // ---- Constant equality RHS (`c_rhs`). For an equality row the
        // bound satisfies `g_l == g_u`; capture it once here so the
        // line-search-hot `eval_c` subtracts a cached constant instead of
        // re-fetching every bound (and allocating four full-size scratch
        // vectors) per cache miss. (Code review 2026-06 item M17.) -----
        let c_rhs: Vec<Number> = classification
            .c_map
            .iter()
            .map(|&g_idx| full_g_l[g_idx as usize])
            .collect();

        // ---- Jacobian sparsity. Ask the TNLP for the full jacobian
        // structure (g rows × full-x cols), then split entries into
        // c-rows (equality) and d-rows (inequality). Within each split
        // the row index is remapped to the new dense indexing in
        // [0, n_c) / [0, n_d).
        //
        // We currently keep all original full-x columns (no fixed-var
        // removal); when fixed-var treatment lands, this is where the
        // column remap goes. -----
        let mut full_irow = vec![0 as Index; info.nnz_jac_g as usize];
        let mut full_jcol = vec![0 as Index; info.nnz_jac_g as usize];
        {
            let a = adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            let ok = t.eval_jac_g(
                None,
                false,
                SparsityRequest::Structure {
                    irow: &mut full_irow,
                    jcol: &mut full_jcol,
                },
            );
            if !ok {
                return Err("TNLP::eval_jac_g(Structure) returned false".into());
            }
        }

        // Build the inverse maps: g-row → c-row (or d-row).
        let mut g_to_c = vec![-1 as Index; n_full_g];
        for (c_idx, &g_idx) in classification.c_map.iter().enumerate() {
            g_to_c[g_idx as usize] = c_idx as Index;
        }
        let mut g_to_d = vec![-1 as Index; n_full_g];
        for (d_idx, &g_idx) in classification.d_map.iter().enumerate() {
            g_to_d[g_idx as usize] = d_idx as Index;
        }

        let style_offset = match info.index_style {
            crate::tnlp::IndexStyle::C => 0 as Index,
            crate::tnlp::IndexStyle::Fortran => 1 as Index,
        };

        let mut jac_c_irow_1based = Vec::new();
        let mut jac_c_jcol_1based = Vec::new();
        let mut jac_c_entry_in_g = Vec::new();
        let mut jac_d_irow_1based = Vec::new();
        let mut jac_d_jcol_1based = Vec::new();
        let mut jac_d_entry_in_g = Vec::new();

        // `make_parameter`: drop Jacobian entries in fixed-variable
        // columns. Their contribution to f and g is constant under the
        // active-x search so they don't appear in the KKT.
        let full_to_var = &classification.full_to_var;
        // Dropping an entry here *is* the observation "this row touches a
        // fixed variable" — see `FixedVarReach`. Recording it costs one
        // byte per row and no extra query.
        let mut g_row_touches_fixed = vec![false; info.m as usize];
        for k in 0..info.nnz_jac_g as usize {
            let g_row_0 = (full_irow[k] - style_offset) as usize;
            let x_col_0 = (full_jcol[k] - style_offset) as usize;
            let var_col = full_to_var[x_col_0];
            if var_col < 0 {
                g_row_touches_fixed[g_row_0] = true;
                continue;
            }
            // Triplet output is 1-based (matches `GenTMatrix` convention).
            let col_1based = var_col + 1;
            let c_row = g_to_c[g_row_0];
            if c_row >= 0 {
                jac_c_irow_1based.push(c_row + 1);
                jac_c_jcol_1based.push(col_1based);
                jac_c_entry_in_g.push(k as Index);
            } else {
                let d_row = g_to_d[g_row_0];
                debug_assert!(d_row >= 0, "g row {g_row_0} is neither in c_map nor d_map");
                jac_d_irow_1based.push(d_row + 1);
                jac_d_jcol_1based.push(col_1based);
                jac_d_entry_in_g.push(k as Index);
            }
        }

        let jac_c_space = GenTMatrixSpace::new(
            classification.n_c,
            n_x_var,
            jac_c_irow_1based,
            jac_c_jcol_1based,
        );
        let jac_d_space = GenTMatrixSpace::new(
            classification.n_d,
            n_x_var,
            jac_d_irow_1based,
            jac_d_jcol_1based,
        );

        // ---- Hessian sparsity (optional). If the TNLP doesn't
        // implement `eval_h`, we leave `h_space = None`. The Phase-8
        // limited-memory quasi-Newton path will populate it from
        // `LowRankUpdateSymMatrixSpace` instead. -----
        let nnz_h_lag_full = info.nnz_h_lag;
        let mut h_entry_in_full: Vec<Index> = Vec::new();
        // A TNLP with no `eval_h` leaves second-order reach unanswerable,
        // so it starts out assumed. A structurally zero Hessian answers it:
        // no second derivatives at all means no derivative w.r.t. the
        // remaining variables can depend on a fixed one.
        let mut second_order_touches_fixed = info.nnz_h_lag > 0;
        let h_space = if info.nnz_h_lag > 0 {
            let mut h_irow = vec![0 as Index; info.nnz_h_lag as usize];
            let mut h_jcol = vec![0 as Index; info.nnz_h_lag as usize];
            let supports_h = {
                let a = adapter.borrow();
                let mut t = a.tnlp().borrow_mut();
                t.eval_h(
                    None,
                    false,
                    1.0,
                    None,
                    false,
                    SparsityRequest::Structure {
                        irow: &mut h_irow,
                        jcol: &mut h_jcol,
                    },
                )
            };
            if supports_h {
                second_order_touches_fixed = false;
                // `make_parameter`: drop Hessian entries where row OR
                // column is fixed (the second derivatives w.r.t. a
                // parameter are not needed in the active-x KKT). The
                // surviving entries are remapped from full-x indices
                // to var-x indices via `full_to_var`.
                let mut h_irow_1: Vec<Index> = Vec::with_capacity(h_irow.len());
                let mut h_jcol_1: Vec<Index> = Vec::with_capacity(h_jcol.len());
                for k in 0..h_irow.len() {
                    let i_full = (h_irow[k] - style_offset) as usize;
                    let j_full = (h_jcol[k] - style_offset) as usize;
                    let i_var = full_to_var[i_full];
                    let j_var = full_to_var[j_full];
                    if i_var < 0 || j_var < 0 {
                        second_order_touches_fixed = true;
                        continue;
                    }
                    h_irow_1.push(i_var + 1);
                    h_jcol_1.push(j_var + 1);
                    h_entry_in_full.push(k as Index);
                }
                Some(SymTMatrixSpace::new(n_x_var, h_irow_1, h_jcol_1))
            } else {
                // TODO(Phase 8): wire the L-BFGS / SR1 path here.
                None
            }
        } else {
            // LPs and other problems with structurally zero Hessian: build an
            // empty SymTMatrixSpace so eval_h_internal returns a zero matrix
            // rather than panicking down the L-BFGS error path.
            Some(SymTMatrixSpace::new(n_x_var, Vec::new(), Vec::new()))
        };

        // Honor the scaling object's constant factor from construction
        // (a negative `obj_scaling_factor` means maximize). Callers
        // that run `determine_scaling_from_starting_point` overwrite
        // this with the combined automatic·user factor; callers that
        // don't (e.g. the SQP path) still get the user's constant.
        let initial_obj_scal = scaling.obj_scaling();
        Ok(Self {
            adapter,
            scaling,
            obj_scale_factor: Cell::new(initial_obj_scal),
            computed_obj_scale: Cell::new(1.0),
            c_scale: RefCell::new(None),
            d_scale: RefCell::new(None),
            declared_d_l: RefCell::new(None),
            declared_d_u: RefCell::new(None),
            declared_x_l: RefCell::new(None),
            declared_x_u: RefCell::new(None),
            honor_original_bounds: Cell::new(false),
            x_scaling_rejected: Cell::new(false),
            x_space,
            c_space,
            d_space,
            x_l_space,
            x_u_space,
            d_l_space,
            d_u_space,
            px_l_space,
            px_u_space,
            pd_l_space,
            pd_u_space,
            jac_c_space,
            jac_d_space,
            h_space,
            x_l: Rc::new(x_l),
            x_u: Rc::new(x_u),
            d_l: Rc::new(d_l),
            d_u: Rc::new(d_u),
            c_rhs,
            warm_start_snapshot: RefCell::new(None),
            px_l,
            px_u,
            pd_l,
            pd_u,
            jac_c_entry_in_g,
            jac_d_entry_in_g,
            nnz_jac_g_full: info.nnz_jac_g,
            nnz_h_lag_full,
            h_entry_in_full,
            f_cache: RefCell::new(Cache::new(1)),
            grad_f_cache: RefCell::new(Cache::new(1)),
            c_cache: RefCell::new(Cache::new(1)),
            d_cache: RefCell::new(Cache::new(1)),
            jac_c_cache: RefCell::new(Cache::new(1)),
            jac_d_cache: RefCell::new(Cache::new(1)),
            h_cache: RefCell::new(Cache::new(1)),
            full_g_cache: RefCell::new(Cache::new(1)),
            full_jac_g_cache: RefCell::new(Cache::new(1)),
            f_evals: RefCell::new(0),
            grad_f_evals: RefCell::new(0),
            c_evals: RefCell::new(0),
            d_evals: RefCell::new(0),
            jac_c_evals: RefCell::new(0),
            jac_d_evals: RefCell::new(0),
            h_evals: RefCell::new(0),
            info,
            timing: RefCell::new(None),
            const_deriv: ConstantDerivatives::default(),
            fixed_reach: FixedVarReach {
                g_row_touches_fixed,
                second_order_touches_fixed,
            },
        })
    }

    /// Install the shared timing accumulator. `IpoptApplication` calls
    /// this once per solve so each `eval_*` entrypoint records into the
    /// same `TimingStatistics` instance the algorithm reports at the
    /// end of the run. Calling with `None` (or never calling) leaves
    /// timing disabled.
    pub fn set_timing_stats(&self, t: Rc<TimingStatistics>) {
        *self.timing.borrow_mut() = Some(t);
    }

    /// Run `f` with two timers active for the duration of the call:
    /// `pick(&timing)` (e.g. `eval_obj`) and `total_function_evaluation_time`.
    /// When no `TimingStatistics` is installed, the closure is invoked
    /// directly with no overhead.
    fn timed_eval<R, F>(&self, pick: fn(&TimingStatistics) -> &pounce_common::TimedTask, f: F) -> R
    where
        F: FnOnce() -> R,
    {
        let guard = self.timing.borrow();
        match guard.as_deref() {
            Some(t) => {
                let task = pick(t);
                task.start();
                t.total_function_evaluation_time.start();
                let r = f();
                t.total_function_evaluation_time.end();
                task.end();
                r
            }
            None => {
                drop(guard);
                f()
            }
        }
    }

    // ---- constant-derivative hints (gh #588, Q6) ----

    /// Ask the model what it can prove about its own derivatives, and
    /// fold the per-row answers into the four hints the algorithm's
    /// evaluation sites are actually about.
    ///
    /// Returned in [`crate::constant_derivatives::HINT_OPTIONS`] order:
    /// `grad_f`, `hessian`, `jac_c`, `jac_d`. Two translations happen
    /// here and nowhere else, because this is the only layer that knows
    /// both the model's row order and the algorithm's:
    ///
    /// 1. **The c/d split.** The model proves things about *its* rows;
    ///    the algorithm evaluates an equality block and an inequality
    ///    block. `c_map` / `d_map` are the authoritative split — deriving
    ///    it a second time from `g_l == g_u` inside the model would be a
    ///    second rule that could disagree with this one.
    /// 2. **Fixed variables.** Under `fixed_variable_treatment =
    ///    make_parameter` a variable with `x_l == x_u` leaves the space
    ///    the algorithm differentiates in, so a row like `x·y` with `y`
    ///    fixed has a *constant* gradient here while the model correctly
    ///    proved it varies. A proof of variation therefore does not
    ///    survive fixing and is weakened to `Unknown`; a proof of
    ///    constancy does survive, since dropping variables cannot
    ///    introduce dependence. Getting this backwards would make pounce
    ///    refuse a hint that is true — the same kind of silent wrong
    ///    answer the phase exists to prevent, pointed the other way.
    pub fn derivative_proofs(&self) -> [DerivativeProof; 4] {
        let mut proofs: DerivativeProofs = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.derivative_proofs()
        };
        let cls = self.adapter.borrow().classification().clone();
        if cls.n_x_fixed > 0 {
            // A proof is stated about the model the *user* wrote; pounce
            // then solves a smaller one, with the fixed variables folded
            // into the forms as constants. That can only remove variation,
            // so `Varying` no longer holds — but only where a fixed
            // variable is reachable. Weakening a form it cannot reach
            // would convert a refusal into `Unknown`, which the hint layer
            // honours **on trust**, and a model that fixes one unrelated
            // variable would quietly lose the ability to refuse a false
            // `hessian_constant=yes` on the rest of itself. See
            // [`FixedVarReach`] for why each test is exactly right.
            //
            // Weakening happens per row, *before* the fold below, so a
            // subsystem is only weakened by the rows that earned it.
            let reach = &self.fixed_reach;
            for (row, proof) in proofs.jac.iter_mut().enumerate() {
                if reach.g_row_touches_fixed.get(row).copied().unwrap_or(true) {
                    *proof = proof.forget_variation();
                }
            }
            if reach.second_order_touches_fixed {
                proofs.grad_f = proofs.grad_f.forget_variation();
                proofs.hessian = proofs.hessian.forget_variation();
            }
        }
        [
            proofs.grad_f,
            proofs.hessian,
            subsystem_proof(&proofs, &cls.c_map),
            subsystem_proof(&proofs, &cls.d_map),
        ]
    }

    /// Install the resolved hints. Called once, before the first
    /// evaluation; calling it later would leave point-keyed entries in
    /// the caches, which is harmless (they simply never match again) but
    /// means the first reused value is whichever iterate got there
    /// first.
    pub fn set_constant_derivatives(&mut self, cd: ConstantDerivatives) {
        self.const_deriv = cd;
        // Anything cached against a point was computed before the
        // decision and must not be mistaken for the constant answer.
        self.invalidate_eval_caches();
    }

    /// The hints in force for this solve.
    pub fn constant_derivatives(&self) -> ConstantDerivatives {
        self.const_deriv
    }

    // ---- accessors used by the algorithm wiring layer ----

    pub fn nlp_info(&self) -> &NlpInfo {
        &self.info
    }
    pub fn classification_n_x_var(&self) -> Index {
        self.x_space.dim()
    }
    pub fn x_space(&self) -> &Rc<DenseVectorSpace> {
        &self.x_space
    }
    pub fn c_space(&self) -> &Rc<DenseVectorSpace> {
        &self.c_space
    }
    pub fn d_space(&self) -> &Rc<DenseVectorSpace> {
        &self.d_space
    }
    pub fn x_l_space(&self) -> &Rc<DenseVectorSpace> {
        &self.x_l_space
    }
    pub fn x_u_space(&self) -> &Rc<DenseVectorSpace> {
        &self.x_u_space
    }
    pub fn d_l_space(&self) -> &Rc<DenseVectorSpace> {
        &self.d_l_space
    }
    pub fn d_u_space(&self) -> &Rc<DenseVectorSpace> {
        &self.d_u_space
    }
    pub fn px_l_space(&self) -> &Rc<ExpansionMatrixSpace> {
        &self.px_l_space
    }
    pub fn px_u_space(&self) -> &Rc<ExpansionMatrixSpace> {
        &self.px_u_space
    }
    pub fn pd_l_space(&self) -> &Rc<ExpansionMatrixSpace> {
        &self.pd_l_space
    }
    pub fn pd_u_space(&self) -> &Rc<ExpansionMatrixSpace> {
        &self.pd_u_space
    }
    pub fn jac_c_space(&self) -> &Rc<GenTMatrixSpace> {
        &self.jac_c_space
    }
    pub fn jac_d_space(&self) -> &Rc<GenTMatrixSpace> {
        &self.jac_d_space
    }
    pub fn h_space(&self) -> Option<&Rc<SymTMatrixSpace>> {
        self.h_space.as_ref()
    }

    /// Effective objective scaling factor (`df_` upstream). 1.0 when
    /// no scaling has been determined.
    pub fn obj_scale_factor(&self) -> Number {
        self.obj_scale_factor.get()
    }

    /// Apply `bound_relax_factor` to the unscaled `x_L / x_U / d_L / d_U`
    /// in place. Mirrors `OrigIpoptNLP::relax_bounds`
    /// (`IpOrigIpoptNLP.cpp:343-358, 459-481`):
    ///
    /// `delta_i = min(constr_viol_tol, |relax| * max(|bound_i|, 1))`,
    /// then `x_L -= delta`, `x_U += delta`, `d_L -= delta`, `d_U += delta`.
    ///
    /// Must be called before `determine_scaling_from_starting_point`
    /// (which only reads bounds via cached evals — so order doesn't
    /// affect scaling — but the bounds themselves should be the
    /// post-relax values when they enter the algorithm).
    /// Record the bounds **as the caller declared them**, before anything
    /// widens them.
    ///
    /// Taken unconditionally — even when no relaxation will be applied — so
    /// [`Self::declared_d_bounds`] and [`Self::declared_box_violation`] have
    /// one authoritative answer per solve rather than an answer that exists
    /// only on widened runs. Safe-slack adjustments come later and only touch
    /// the live vectors, so this stays the user's own box.
    ///
    /// Split out of [`Self::relax_bounds`] for the active-set SQP arm, which
    /// applies no widening at all and so never called it. The cost of that was
    /// `Variable bound violation` reading `nan` on every SQP solve: the row is
    /// computed from this snapshot, and without it the accessor abstains. An
    /// arm that does not widen still has a declared box, and the distance the
    /// answer sits outside it is still the number that row exists to report.
    pub fn snapshot_declared_bounds(&mut self) {
        *self.declared_d_l.borrow_mut() = Some(self.d_l.expanded_values());
        *self.declared_d_u.borrow_mut() = Some(self.d_u.expanded_values());
        *self.declared_x_l.borrow_mut() = Some(self.x_l.expanded_values());
        *self.declared_x_u.borrow_mut() = Some(self.x_u.expanded_values());
    }

    pub fn relax_bounds(&mut self, bound_relax_factor: Number, constr_viol_tol: Number) {
        self.snapshot_declared_bounds();
        if bound_relax_factor <= 0.0 {
            return;
        }
        let relax = bound_relax_factor.abs();
        let cap = constr_viol_tol;
        let apply = |v: &mut DenseVector, sign: Number| {
            let xs = v.values_mut();
            for x in xs.iter_mut() {
                let delta = (relax * x.abs().max(1.0)).min(cap);
                *x += sign * delta;
            }
        };
        // Inequality-row bounds use a *scale-relative* delta (#385, Step 6):
        // `min(relax, cap) · |b|`, with the absolute `min(relax, cap)` kept
        // only for a declared-zero bound (`s·g >= 0` is the same row at every
        // `s`, so zero has no scale and the absolute form is already
        // invariant there). The upstream formula's `max(|b|, 1)` clamp is the
        // same absolute floor this migration removes everywhere else: it
        // relaxed a `2e-12`-bound row by `1e-8` — 5000× the bound — silently
        // erasing every down-scaled constraint before the solver saw it,
        // which is exactly how an infeasible model at row scale `1e-8` and
        // below reported `Solve_Succeeded`. In the other direction the
        // absolute `cap` pinned a `1e13`-bound row's relaxation at `1e-4` —
        // relative `1e-17`, i.e. no relaxation at all. Both directions are
        // now relative: `min(relax, cap)` is the *relative* width, identical
        // to the upstream formula on `1 <= |b| <= cap/relax`, which is where
        // the corpus lives. Variable bounds keep the upstream formula — row
        // scaling never touches them, and their relaxation is not this
        // change's business.
        let apply_d = |v: &mut DenseVector, sign: Number| {
            let rel_width = relax.min(cap);
            let xs = v.values_mut();
            for x in xs.iter_mut() {
                let scale = if *x == 0.0 { 1.0 } else { x.abs() };
                *x += sign * rel_width * scale;
            }
        };
        // The bound `Rc`s are uniquely owned (nothing clones them — the same
        // invariant `adjust_variable_bounds` relies on), so `get_mut` must
        // succeed. A shared `Rc` here would silently skip the relaxation,
        // leaving bounds tighter than `bound_relax_factor` requires; that is
        // a programming error, so fail loudly to match `adjust_variable_bounds`
        // rather than no-op.
        apply(
            Rc::get_mut(&mut self.x_l).expect("relax_bounds: x_l is uniquely owned"),
            -1.0,
        );
        apply(
            Rc::get_mut(&mut self.x_u).expect("relax_bounds: x_u is uniquely owned"),
            1.0,
        );
        apply_d(
            Rc::get_mut(&mut self.d_l).expect("relax_bounds: d_l is uniquely owned"),
            -1.0,
        );
        apply_d(
            Rc::get_mut(&mut self.d_u).expect("relax_bounds: d_u is uniquely owned"),
            1.0,
        );
    }

    /// Determine objective + per-constraint scaling from the starting
    /// point, per `Algorithm/IpGradientScaling.cpp::DetermineScalingParametersImpl`
    /// (and now also `nlp_scaling_method=user-scaling`). Should be called
    /// once, after construction and before the algorithm enters its main
    /// loop.
    ///
    /// Arguments:
    /// * `method` — `None` / `GradientBased` / `UserScaling`.
    /// * `max_gradient` — `nlp_scaling_max_gradient` (cutoff above which
    ///   gradient-based scaling fires; default 100).
    /// * `min_value` — `nlp_scaling_min_value` (floor on computed scale
    ///   factors; default 1e-8).
    /// * `obj_target_gradient` — `nlp_scaling_obj_target_gradient`
    ///   (default 0; when `> 0`, fixes `df = obj_target_gradient /
    ///   max_grad_f` unconditionally, overriding the cutoff).
    /// * `constr_target_gradient` — `nlp_scaling_constr_target_gradient`
    ///   (default 0; when `> 0`, fixes per-row scale to
    ///   `constr_target_gradient / row_max` unconditionally).
    ///
    /// Cache state is invalidated so subsequent eval calls produce
    /// scaled values.
    pub fn determine_scaling_from_starting_point(
        &mut self,
        method: ScalingMethod,
        max_gradient: Number,
        min_value: Number,
        obj_target_gradient: Number,
        constr_target_gradient: Number,
    ) {
        // Always pull the user's `obj_scaling_factor` constant first;
        // it multiplies whatever the automatic scheme computes.
        let user_obj_factor = self.scaling.obj_scaling();
        if matches!(method, ScalingMethod::None) {
            self.obj_scale_factor.set(user_obj_factor);
            *self.c_scale.borrow_mut() = None;
            *self.d_scale.borrow_mut() = None;
            self.invalidate_eval_caches();
            return;
        }

        // ---- Get starting x_full (needed by both gradient + user paths) ----
        let cls = self.adapter.borrow().classification().clone();
        let n_full_x = cls.n_full_x as usize;
        let n_full_g = cls.n_full_g as usize;
        let mut full_x = vec![0.0; n_full_x];
        let mut full_z_l = vec![0.0; n_full_x];
        let mut full_z_u = vec![0.0; n_full_x];
        let mut full_lambda = vec![0.0; n_full_g];
        let starting_ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.get_starting_point(StartingPoint {
                init_x: true,
                x: &mut full_x,
                init_z: false,
                z_l: &mut full_z_l,
                z_u: &mut full_z_u,
                init_lambda: false,
                lambda: &mut full_lambda,
            })
        };
        if !starting_ok {
            // Fall back to no automatic scaling.
            self.obj_scale_factor.set(user_obj_factor);
            *self.c_scale.borrow_mut() = None;
            *self.d_scale.borrow_mut() = None;
            self.invalidate_eval_caches();
            return;
        }

        // Lift fixed variables (x_l == x_u) to their fixed value before
        // sampling the gradient / Jacobian. Fixed vars never enter the
        // algorithm's compressed x; every algorithm-side eval re-inserts
        // their fixed value via `lift_x_to_full`, so scaling must be
        // computed at that same point. Upstream achieves this implicitly:
        // `TNLPAdapter::GetStartingPoint` projects the start onto the
        // (relaxed) bounds, pinning fixed vars to their value. A raw `x0`
        // that leaves them elsewhere can shift the objective gradient by
        // orders of magnitude (pounce: flosp2hm — 41 fixed vars sitting at
        // x0=0 instead of their fixed value 1 made ‖∇f‖∞ read 40 instead of
        // 2.4e5, so obj_scale_factor stayed 1.0 and the solve stalled at
        // max-iter while IPOPT, scaling correctly, converged in 5 iters).
        for (i, &full_idx) in cls.x_fixed_map.iter().enumerate() {
            full_x[full_idx as usize] = cls.x_fixed_vals[i];
        }

        match method {
            ScalingMethod::None => unreachable!("handled above"),
            ScalingMethod::GradientBased => {
                self.scale_gradient_based(
                    &cls,
                    &full_x,
                    user_obj_factor,
                    max_gradient,
                    min_value,
                    obj_target_gradient,
                    constr_target_gradient,
                );
            }
            ScalingMethod::UserScaling => {
                let applied = self.scale_user_supplied(&cls, user_obj_factor, min_value);
                if !applied {
                    // TNLP declined to supply scaling — fall through to
                    // no automatic scaling (matches upstream's behavior
                    // when `get_scaling_parameters` returns false).
                    self.obj_scale_factor.set(user_obj_factor);
                    *self.c_scale.borrow_mut() = None;
                    *self.d_scale.borrow_mut() = None;
                }
            }
        }

        // Apply the d-row scaling to the d_l/d_u bound vectors so
        // feasibility checks compare like with like (gh#54).
        self.apply_d_scale_to_bounds();

        // Drop any cached eval results computed before the scales were
        // set (their values would be wrong now).
        self.invalidate_eval_caches();
    }

    /// Gradient-based pathway: compute `df_`, `dc_`, `dd_` from the
    /// objective gradient and constraint Jacobian at the starting point.
    fn scale_gradient_based(
        &self,
        cls: &BoundClassification,
        full_x: &[Number],
        user_obj_factor: Number,
        max_gradient: Number,
        min_value: Number,
        obj_target_gradient: Number,
        constr_target_gradient: Number,
    ) {
        let n_full_x = cls.n_full_x as usize;
        let n_full_g = cls.n_full_g as usize;

        // ---- Objective gradient scale ----
        let mut full_grad_f = vec![0.0; n_full_x];
        let grad_ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_grad_f(full_x, true, &mut full_grad_f)
        };
        let mut df = 1.0;
        if grad_ok {
            // Amax over the *compressed* x_var space (matches upstream
            // which scales the algorithm-side gradient).
            let mut max_grad_f: Number = 0.0;
            for &full_idx in cls.x_not_fixed_map.iter() {
                let v = full_grad_f[full_idx as usize].abs();
                if v > max_grad_f {
                    max_grad_f = v;
                }
            }
            df = gradient_obj_scale(max_grad_f, max_gradient, min_value, obj_target_gradient);
        }
        self.computed_obj_scale.set(df);
        self.obj_scale_factor.set(df * user_obj_factor);

        // ---- Constraint Jacobian row-max scaling ----
        if cls.n_full_g == 0 {
            *self.c_scale.borrow_mut() = None;
            *self.d_scale.borrow_mut() = None;
            return;
        }
        // Evaluate full Jacobian once at x.
        let mut full_jac_vals = vec![0.0; self.nnz_jac_g_full as usize];
        let jac_ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_jac_g(
                Some(full_x),
                true,
                SparsityRequest::Values {
                    values: &mut full_jac_vals,
                },
            )
        };
        if !jac_ok {
            *self.c_scale.borrow_mut() = None;
            *self.d_scale.borrow_mut() = None;
            return;
        }
        // Recover row indices from the sparsity structure.
        let mut full_irow = vec![0 as Index; self.nnz_jac_g_full as usize];
        let mut full_jcol = vec![0 as Index; self.nnz_jac_g_full as usize];
        let _ = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_jac_g(
                None,
                false,
                SparsityRequest::Structure {
                    irow: &mut full_irow,
                    jcol: &mut full_jcol,
                },
            )
        };
        let style_offset: Index = match self.info.index_style {
            crate::tnlp::IndexStyle::C => 0,
            crate::tnlp::IndexStyle::Fortran => 1,
        };
        // Build inverse row maps to assign each entry to c or d.
        let mut g_to_c = vec![-1 as Index; n_full_g];
        for (c_idx, &g_idx) in cls.c_map.iter().enumerate() {
            g_to_c[g_idx as usize] = c_idx as Index;
        }
        let mut g_to_d = vec![-1 as Index; n_full_g];
        for (d_idx, &g_idx) in cls.d_map.iter().enumerate() {
            g_to_d[g_idx as usize] = d_idx as Index;
        }
        let n_c = cls.n_c as usize;
        let n_d = cls.n_d as usize;
        // Initialize row-max arrays to dbl_min as upstream does.
        let dbl_min = Number::MIN_POSITIVE;
        let mut c_row_max: Vec<Number> = vec![dbl_min; n_c];
        let mut d_row_max: Vec<Number> = vec![dbl_min; n_d];
        for k in 0..self.nnz_jac_g_full as usize {
            let g_row_0 = (full_irow[k] - style_offset) as usize;
            let v = full_jac_vals[k].abs();
            let cr = g_to_c[g_row_0];
            if cr >= 0 {
                let row = cr as usize;
                if v > c_row_max[row] {
                    c_row_max[row] = v;
                }
            } else {
                let dr = g_to_d[g_row_0];
                if dr >= 0 {
                    let row = dr as usize;
                    if v > d_row_max[row] {
                        d_row_max[row] = v;
                    }
                }
            }
        }

        let row_max_to_scale = |row_max: Number| -> Number {
            gradient_row_scale(row_max, max_gradient, min_value, constr_target_gradient)
        };
        let any_row_above = |rows: &[Number]| -> bool {
            gradient_scaling_fires(rows, max_gradient, constr_target_gradient)
        };

        if n_c > 0 && any_row_above(&c_row_max) {
            let dc: Vec<Number> = c_row_max.iter().map(|&v| row_max_to_scale(v)).collect();
            *self.c_scale.borrow_mut() = Some(dc);
        } else {
            *self.c_scale.borrow_mut() = None;
        }

        if n_d > 0 && any_row_above(&d_row_max) {
            let dd: Vec<Number> = d_row_max.iter().map(|&v| row_max_to_scale(v)).collect();
            *self.d_scale.borrow_mut() = Some(dd);
        } else {
            *self.d_scale.borrow_mut() = None;
        }
    }

    /// User-supplied scaling pathway: call `TNLP::get_scaling_parameters`
    /// and translate the user's `obj_scaling` and `g_scaling` arrays
    /// into the algorithm-side `obj_scale_factor`, `c_scale`, `d_scale`.
    /// Returns `true` if the TNLP supplied scaling (matches upstream's
    /// `GetScalingParameters` return-value contract).
    ///
    /// `OrigIpoptNlp` does not model per-variable rescaling (that would
    /// require transforming `eval_grad_f`, `eval_jac_*`, and `eval_h` in
    /// concert); issue #61's `nlp_scaling=user` design covers only
    /// `obj_scale` and `con_scale`. A **non-trivial** `x_scaling`
    /// request is therefore *rejected*, not dropped: it sets
    /// [`Self::x_scaling_rejected`] so the driver can fail the solve
    /// with a message. Quietly discarding it used to hand back a
    /// problem conditioned differently from the one the caller
    /// described, with nothing in the log to say so (gh#483). An
    /// all-ones request is a genuine no-op and passes through.
    fn scale_user_supplied(
        &self,
        cls: &BoundClassification,
        user_obj_factor: Number,
        min_value: Number,
    ) -> bool {
        let n_full_x = cls.n_full_x as usize;
        let n_full_g = cls.n_full_g as usize;
        let mut obj_scaling: Number = 1.0;
        let mut use_x_scaling = false;
        let mut x_scaling = vec![1.0; n_full_x];
        let mut use_g_scaling = false;
        let mut g_scaling = vec![1.0; n_full_g];
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.get_scaling_parameters(ScalingRequest {
                obj_scaling: &mut obj_scaling,
                use_x_scaling: &mut use_x_scaling,
                x_scaling: &mut x_scaling,
                use_g_scaling: &mut use_g_scaling,
                g_scaling: &mut g_scaling,
            })
        };
        if !ok {
            return false;
        }

        // Objective: user's obj_scaling combined with the constant
        // `obj_scaling_factor` (matches upstream's
        // `StandardScalingBase::DetermineScaling`).
        let mut df = obj_scaling;
        if df.abs() < min_value {
            // Defensively floor — a zero/near-zero obj scale would
            // make all duals divide-by-zero on the way out.
            df = df.signum().max(0.0).max(1.0) * min_value;
        }
        self.obj_scale_factor.set(df * user_obj_factor);

        // Constraint vector: split user g_scaling into c_scale / d_scale.
        if use_g_scaling && g_scaling.len() == n_full_g {
            let n_c = cls.n_c as usize;
            let n_d = cls.n_d as usize;
            let mut dc = vec![1.0; n_c];
            for (c_idx, &g_idx) in cls.c_map.iter().enumerate() {
                let s = g_scaling[g_idx as usize];
                dc[c_idx] = if s < min_value { min_value } else { s };
            }
            let mut dd = vec![1.0; n_d];
            for (d_idx, &g_idx) in cls.d_map.iter().enumerate() {
                let s = g_scaling[g_idx as usize];
                dd[d_idx] = if s < min_value { min_value } else { s };
            }
            // Only install the vectors when not all-ones (matches the
            // `Option::None ↔ identity` convention used elsewhere).
            let nontrivial_c = dc.iter().any(|&s| s != 1.0);
            *self.c_scale.borrow_mut() = if nontrivial_c && n_c > 0 {
                Some(dc)
            } else {
                None
            };
            let nontrivial_d = dd.iter().any(|&s| s != 1.0);
            *self.d_scale.borrow_mut() = if nontrivial_d && n_d > 0 {
                Some(dd)
            } else {
                None
            };
        } else {
            *self.c_scale.borrow_mut() = None;
            *self.d_scale.borrow_mut() = None;
        }
        // Per-variable factors are not modeled. Flag a request that
        // would actually change the problem so the driver can refuse
        // loudly; an all-ones vector asks for nothing and is accepted.
        if use_x_scaling && x_scaling.iter().any(|&s| s != 1.0) {
            self.x_scaling_rejected.set(true);
        }
        true
    }

    /// `true` when the last [`Self::determine_scaling_from_starting_point`]
    /// ran `user-scaling` and the TNLP asked for per-variable scaling
    /// factors that pounce cannot honor (see [`Self::scale_user_supplied`]).
    /// Drivers must turn this into a hard error rather than solve a
    /// problem the caller did not describe (gh#483).
    pub fn user_x_scaling_rejected(&self) -> bool {
        self.x_scaling_rejected.get()
    }

    /// Bring `d_l` / `d_u` into the scaled space so feasibility checks
    /// compare like with like (gh#54). Upstream's
    /// `OrigIpoptNLP::Initialize` does this via
    /// `Pd_L_->TransMultVector(scaling.apply_vec_d(...))`.
    fn apply_d_scale_to_bounds(&mut self) {
        let cls = self.adapter.borrow().classification().clone();
        if let Some(dd) = self.d_scale.borrow().as_ref() {
            if let Some(d_l) = Rc::get_mut(&mut self.d_l) {
                let xs = d_l.values_mut();
                for (i, slot) in xs.iter_mut().enumerate() {
                    let d_idx = cls.d_l_map[i] as usize;
                    *slot *= dd[d_idx];
                }
            }
            if let Some(d_u) = Rc::get_mut(&mut self.d_u) {
                let xs = d_u.values_mut();
                for (i, slot) in xs.iter_mut().enumerate() {
                    let d_idx = cls.d_u_map[i] as usize;
                    *slot *= dd[d_idx];
                }
            }
        }
    }

    fn invalidate_eval_caches(&self) {
        self.f_cache.borrow_mut().clear();
        self.grad_f_cache.borrow_mut().clear();
        self.c_cache.borrow_mut().clear();
        self.d_cache.borrow_mut().clear();
        self.jac_c_cache.borrow_mut().clear();
        self.jac_d_cache.borrow_mut().clear();
        self.h_cache.borrow_mut().clear();
    }

    pub fn f_evals(&self) -> Index {
        *self.f_evals.borrow()
    }
    pub fn grad_f_evals(&self) -> Index {
        *self.grad_f_evals.borrow()
    }
    pub fn c_evals(&self) -> Index {
        *self.c_evals.borrow()
    }
    pub fn d_evals(&self) -> Index {
        *self.d_evals.borrow()
    }
    pub fn jac_c_evals(&self) -> Index {
        *self.jac_c_evals.borrow()
    }
    pub fn jac_d_evals(&self) -> Index {
        *self.jac_d_evals.borrow()
    }
    pub fn h_evals(&self) -> Index {
        *self.h_evals.borrow()
    }

    /// Lift a compressed `x_var` (length `n_x_var`) up to the full TNLP
    /// `x` (length `n_full_x`), inserting `x_fixed_vals` at the
    /// `x_fixed_map` positions. Mirrors upstream
    /// `IpTNLPAdapter::ResortX` under `fixed_variable_treatment =
    /// make_parameter`.
    pub fn lift_x_to_full(&self, x: &dyn Vector) -> Vec<Number> {
        let Some(dx) = x.as_any().downcast_ref::<DenseVector>() else {
            panic!("OrigIpoptNlp expects DenseVector for x");
        };
        let a = self.adapter.borrow();
        let cls = a.classification();
        let mut full = vec![0.0; cls.n_full_x as usize];
        let vals = dx.expanded_values();
        for (var_idx, &full_idx) in cls.x_not_fixed_map.iter().enumerate() {
            full[full_idx as usize] = vals[var_idx];
        }
        for (i, &full_idx) in cls.x_fixed_map.iter().enumerate() {
            full[full_idx as usize] = cls.x_fixed_vals[i];
        }
        full
    }

    /// Select `honor_original_bounds`. Called once from the driver
    /// after [`Self::relax_bounds`], which is what captures the bounds
    /// to project back onto.
    pub fn set_honor_original_bounds(&self, on: bool) {
        self.honor_original_bounds.set(on);
    }

    /// The full-x handed to `TNLP::finalize_solution`: [`Self::lift_x_to_full`],
    /// then — under `honor_original_bounds` — clamped back into the
    /// bounds the user declared.
    ///
    /// `bound_relax_factor` (default `1e-8`) widens the box before the
    /// solve, so a solution pinned to a bound comes back *outside* it:
    /// `min (x−3)²` over `x ∈ [0, 1]` reports `x = 1.0000000094`. That
    /// is upstream's behavior too and is why upstream registers this
    /// option — but pounce registered it and never read it, so there was
    /// no way to turn the projection on (gh#483 follow-up). A value
    /// outside its declared domain is not a cosmetic difference: it
    /// breaks a downstream `sqrt(1 − x)`, a domain assertion, or a Pyomo
    /// `Var` whose bounds the value is loaded back into.
    ///
    /// Only the reported point moves. As upstream documents, the
    /// constraint-violation and complementarity numbers in the summary
    /// are for the **non-projected** point and are left alone.
    pub fn finalize_solution_x(&self, x: &dyn Vector) -> Vec<Number> {
        let mut full = self.lift_x_to_full(x);
        if !self.honor_original_bounds.get() {
            return full;
        }
        let cls = self.adapter.borrow().classification().clone();
        // Fixed variables are spliced in at their exact fixed value, so
        // only the free block can have drifted past a bound.
        if let Some(x_l) = self.declared_x_l.borrow().as_ref() {
            for (i, &var_idx) in cls.x_l_map.iter().enumerate() {
                let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
                if full[full_idx] < x_l[i] {
                    full[full_idx] = x_l[i];
                }
            }
        }
        if let Some(x_u) = self.declared_x_u.borrow().as_ref() {
            for (i, &var_idx) in cls.x_u_map.iter().enumerate() {
                let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
                if full[full_idx] > x_u[i] {
                    full[full_idx] = x_u[i];
                }
            }
        }
        full
    }

    /// Lift the algorithm-side `(y_c, y_d)` multipliers to the user
    /// TNLP's `lambda` array (length `m_full = n_c + n_d`, indexed
    /// by original constraint-row order). Result matches the user's
    /// **unscaled-Lagrangian** convention `min f + λ·g(x)` —
    /// i.e. without the obj_factor that the algorithm threads through
    /// `eval_h`. Mirror of upstream
    /// `IpOrigIpoptNLP::FinalizeSolution`'s `mult_g` packing
    /// (`lambda_user = c_scale * y_c / obj_scale_factor`,
    /// `mu_user = d_scale * y_d / obj_scale_factor`). Used by
    /// `application.rs::finalize_via_orig_nlp` to populate the
    /// `Solution.lambda` slot — pounce#11.
    pub fn finalize_solution_lambda(&self, y_c: &dyn Vector, y_d: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let mut lambda = self.pack_lambda_for_user(y_c, y_d, &cls);
        let obj_scal = self.obj_scale_factor.get();
        if obj_scal != 0.0 && obj_scal != 1.0 {
            let inv = 1.0 / obj_scal;
            for v in lambda.iter_mut() {
                *v *= inv;
            }
        }
        lambda
    }

    /// Lift the algorithm-side compressed `z_l` (length `n_x_l`,
    /// indexed via `x_l_map`) to the user's full-x bound multiplier
    /// array (length `n_full_x`). Slots without a finite lower bound
    /// — including fixed variables — are reported as `0.0`. Sign and
    /// scale match upstream Ipopt: `z_l ≥ 0` for active lower
    /// bounds, divided by `obj_scale_factor` so the user sees the
    /// unscaled-Lagrangian dual.
    pub fn finalize_solution_z_l(&self, z_l: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let n_full_x = cls.n_full_x as usize;
        let mut full_z_l = vec![0.0; n_full_x];
        let n_x_l = self.x_l.dim() as usize;
        if n_x_l == 0 {
            return full_z_l;
        }
        let Some(dz) = z_l.as_any().downcast_ref::<DenseVector>() else {
            panic!("OrigIpoptNlp::finalize_solution_z_l expects DenseVector");
        };
        let vals = dz.expanded_values();
        let obj_scal = self.obj_scale_factor.get();
        let inv = if obj_scal == 0.0 { 1.0 } else { 1.0 / obj_scal };
        for i in 0..n_x_l {
            let var_idx = cls.x_l_map[i] as usize;
            let full_idx = cls.x_not_fixed_map[var_idx] as usize;
            full_z_l[full_idx] = vals[i] * inv;
        }
        full_z_l
    }

    /// Mirror of [`Self::finalize_solution_z_l`] for the upper-bound
    /// duals. Indexed via `x_u_map`.
    pub fn finalize_solution_z_u(&self, z_u: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let n_full_x = cls.n_full_x as usize;
        let mut full_z_u = vec![0.0; n_full_x];
        let n_x_u = self.x_u.dim() as usize;
        if n_x_u == 0 {
            return full_z_u;
        }
        let Some(dz) = z_u.as_any().downcast_ref::<DenseVector>() else {
            panic!("OrigIpoptNlp::finalize_solution_z_u expects DenseVector");
        };
        let vals = dz.expanded_values();
        let obj_scal = self.obj_scale_factor.get();
        let inv = if obj_scal == 0.0 { 1.0 } else { 1.0 / obj_scal };
        for i in 0..n_x_u {
            let var_idx = cls.x_u_map[i] as usize;
            let full_idx = cls.x_not_fixed_map[var_idx] as usize;
            full_z_u[full_idx] = vals[i] * inv;
        }
        full_z_u
    }

    /// Clone the user-provided multipliers (already in the
    /// algorithm's eq/ineq-split form) into a single `lambda` array of
    /// length `m_full = n_c + n_d` ordered by original g-index. Used
    /// by `eval_h` and `finalize_solution`.
    /// Pack the algorithm-side `(y_c, y_d)` multipliers into the user
    /// TNLP's `lambda` array (full-g indexed), applying c/d scale
    /// factors so the result is in the user's unscaled-constraint
    /// multiplier space (`lambda_user_i = c_scale_i * y_c_i`). Used
    /// when invoking the user's `eval_h`.
    pub fn pack_lambda_for_user(
        &self,
        y_c: &dyn Vector,
        y_d: &dyn Vector,
        cls: &BoundClassification,
    ) -> Vec<Number> {
        let mut lambda = vec![0.0; cls.n_full_g as usize];
        if cls.n_c > 0 {
            let Some(dy) = y_c.as_any().downcast_ref::<DenseVector>() else {
                panic!("OrigIpoptNlp expects DenseVector for y_c");
            };
            let vals = dy.expanded_values();
            let cs = self.c_scale.borrow();
            for (i, &g_idx) in cls.c_map.iter().enumerate() {
                lambda[g_idx as usize] = match cs.as_ref() {
                    Some(v) => vals[i] * v[i],
                    None => vals[i],
                };
            }
        }
        if cls.n_d > 0 {
            let Some(dy) = y_d.as_any().downcast_ref::<DenseVector>() else {
                panic!("OrigIpoptNlp expects DenseVector for y_d");
            };
            let vals = dy.expanded_values();
            let ds = self.d_scale.borrow();
            for (i, &g_idx) in cls.d_map.iter().enumerate() {
                lambda[g_idx as usize] = match ds.as_ref() {
                    Some(v) => vals[i] * v[i],
                    None => vals[i],
                };
            }
        }
        lambda
    }

    // -------------------- Initialization --------------------

    fn fetch_warm_start_snapshot(&self) -> Option<StartingPointSnapshot> {
        let cls = self.adapter.borrow().classification().clone();
        // The bound-multiplier slots start *unseeded*, not at zero
        // (gh#622). `TNLP::get_starting_point` is asked for `init_z`
        // and is free to leave the blocks untouched — a caller warm
        // starting from a point alone does exactly that — and this
        // snapshot is what `get_starting_z` then hands the algorithm.
        // Zero is a legal multiplier value, so it sails past the "was
        // this seeded?" resolution in the warm-start initializer and
        // is merely floored at `warm_start_mult_bound_push`: 1e-3 by
        // default, 1e-9 under the tightened pushes `pounce.WarmStart`
        // ships. A start of z = 1e-9 declares every bound inactive and
        // breaks complementarity against mu before the first
        // iteration. NaN is the marker that initializer already
        // documents for "you decide", and resolves to
        // `bound_mult_init_val`.
        //
        // `lambda` deliberately keeps its zero fill. The equality
        // multipliers' unseeded resolution is *also* zero, so the
        // marker would buy nothing, and `any_dual_seeded` reads an
        // all-zero `y` as unseeded — which is what keeps gh#606's
        // reconstruction off a primal-only seed (measured there:
        // 1102 -> 1211 iterations across 27 parametric paths).
        let mut snapshot = StartingPointSnapshot {
            x: vec![0.0; cls.n_full_x as usize],
            z_l: vec![Number::NAN; cls.n_full_x as usize],
            z_u: vec![Number::NAN; cls.n_full_x as usize],
            lambda: vec![0.0; cls.n_full_g as usize],
        };
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.get_starting_point(StartingPoint {
                init_x: true,
                x: &mut snapshot.x,
                init_z: true,
                z_l: &mut snapshot.z_l,
                z_u: &mut snapshot.z_u,
                init_lambda: true,
                lambda: &mut snapshot.lambda,
            })
        };
        ok.then_some(snapshot)
    }

    /// Fill the algorithm's iterate slots with the TNLP's starting
    /// point. Mirrors the second half of upstream
    /// `InitializeStructures`. The caller passes already-allocated
    /// `DenseVector`s in the right spaces; we set them in place.
    ///
    /// Returns the four `init_*` flags so the caller can decide
    /// whether to overwrite zeros with the user's guess.
    #[allow(clippy::too_many_arguments)]
    pub fn initialize_starting_point(
        &mut self,
        x: &mut DenseVector,
        init_x: bool,
        y_c: &mut DenseVector,
        init_y_c: bool,
        y_d: &mut DenseVector,
        init_y_d: bool,
        z_l: &mut DenseVector,
        init_z_l: bool,
        z_u: &mut DenseVector,
        init_z_u: bool,
    ) -> bool {
        let n_full_x = self.adapter.borrow().classification().n_full_x as usize;
        let n_full_g = self.adapter.borrow().classification().n_full_g as usize;
        let n_x_l = self.x_l.dim() as usize;
        let n_x_u = self.x_u.dim() as usize;

        let mut full_x = vec![0.0; n_full_x];
        // Unseeded, not zero, for the reason spelled out in
        // `fetch_warm_start_snapshot` (gh#622): a block the TNLP
        // declines to write must not read as a supplied multiplier of
        // zero. `full_lambda` keeps its zero fill, also per that note.
        let mut full_z_l = vec![Number::NAN; n_full_x];
        let mut full_z_u = vec![Number::NAN; n_full_x];
        let mut full_lambda = vec![0.0; n_full_g];

        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.get_starting_point(StartingPoint {
                init_x,
                x: &mut full_x,
                init_z: init_z_l || init_z_u,
                z_l: &mut full_z_l,
                z_u: &mut full_z_u,
                init_lambda: init_y_c || init_y_d,
                lambda: &mut full_lambda,
            })
        };
        if !ok {
            return false;
        }

        let cls = self.adapter.borrow().classification().clone();
        let obj_scal = self.obj_scale_factor.get();
        let c_scale = self.c_scale.borrow();
        let d_scale = self.d_scale.borrow();

        // Compress full_x → x.
        if init_x {
            let xs = x.values_mut();
            for (var_idx, &full_idx) in cls.x_not_fixed_map.iter().enumerate() {
                xs[var_idx] = full_x[full_idx as usize];
            }
        }
        // Compress full_lambda → y_c, y_d. Upstream
        // (`IpOrigIpoptNLP.cpp:407-429`) divides the user multiplier
        // by the constraint scale (`unapply_vector_scaling_*`) and
        // multiplies by obj_scal so that the algorithm-side y_c sees
        // `(obj_scal / c_scale) * lambda_user`.
        if init_y_c && cls.n_c > 0 {
            let yc = y_c.values_mut();
            for (i, &g_idx) in cls.c_map.iter().enumerate() {
                let cs = c_scale.as_ref().map(|v| v[i]).unwrap_or(1.0);
                yc[i] = full_lambda[g_idx as usize] / cs * obj_scal;
            }
        }
        if init_y_d && cls.n_d > 0 {
            let yd = y_d.values_mut();
            for (i, &g_idx) in cls.d_map.iter().enumerate() {
                let ds = d_scale.as_ref().map(|v| v[i]).unwrap_or(1.0);
                yd[i] = full_lambda[g_idx as usize] / ds * obj_scal;
            }
        }
        // Compress full_z_l, full_z_u → z_l, z_u, indexed via x_l_map / x_u_map.
        if init_z_l && n_x_l > 0 {
            let zl = z_l.values_mut();
            for (i, slot) in zl.iter_mut().enumerate().take(n_x_l) {
                let var_idx = cls.x_l_map[i] as usize;
                let full_idx = cls.x_not_fixed_map[var_idx] as usize;
                *slot = full_z_l[full_idx] * obj_scal;
            }
        }
        if init_z_u && n_x_u > 0 {
            let zu = z_u.values_mut();
            for (i, slot) in zu.iter_mut().enumerate().take(n_x_u) {
                let var_idx = cls.x_u_map[i] as usize;
                let full_idx = cls.x_not_fixed_map[var_idx] as usize;
                *slot = full_z_u[full_idx] * obj_scal;
            }
        }
        true
    }

    // -------------------- Internal eval helpers --------------------

    fn eval_f_internal(&self, x: &dyn Vector) -> Number {
        if let Some(v) = self.f_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        *self.f_evals.borrow_mut() += 1;
        let full_x = self.lift_x_to_full(x);
        let unscaled = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            // A failed user eval (domain error, e.g. log of a negative) is
            // upstream Ipopt's `Eval_Error`. Return NaN so the line search's
            // non-finite-trial path backtracks the step, rather than aborting
            // — and a panic cannot unwind across the C FFI boundary anyway.
            t.eval_f(&full_x, true).unwrap_or(f64::NAN)
        };
        let scaled = unscaled * self.obj_scale_factor.get();
        self.f_cache.borrow_mut().add_1dep(scaled, x.as_tagged());
        scaled
    }

    fn eval_grad_f_internal(&self, x: &dyn Vector) -> Rc<dyn Vector> {
        // A reused derivative lives in the same cache under an empty
        // dependency list, so it matches at every point (gh #588, Q6).
        if self.const_deriv.grad_f
            && let Some(v) = self.grad_f_cache.borrow().get(&[], &[])
        {
            return v;
        }
        if let Some(v) = self.grad_f_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        *self.grad_f_evals.borrow_mut() += 1;
        let full_x = self.lift_x_to_full(x);
        let mut full_g = vec![0.0; full_x.len()];
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_grad_f(&full_x, true, &mut full_g)
        };
        // Eval failure → NaN-filled gradient, which propagates a non-finite
        // step the line search rejects (see `eval_f_internal`).
        if !ok {
            full_g.fill(f64::NAN);
        }
        // Compress full_g → grad in x_var-space, scale by obj_scal.
        let cls = self.adapter.borrow().classification().clone();
        let mut g_compressed = self.x_space.make_new_dense();
        let obj_scal = self.obj_scale_factor.get();
        {
            let gv = g_compressed.values_mut();
            for (var_idx, &full_idx) in cls.x_not_fixed_map.iter().enumerate() {
                gv[var_idx] = full_g[full_idx as usize] * obj_scal;
            }
        }
        // A failed evaluation NaN-fills; storing that as *the* constant
        // answer would poison every remaining iteration instead of
        // letting the line search back away from a bad point, so the
        // point-keyed entry is used and the next call tries again.
        let reuse = self.const_deriv.grad_f && all_finite(g_compressed.values());
        let result: Rc<dyn Vector> = Rc::new(g_compressed);
        if reuse {
            self.grad_f_cache
                .borrow_mut()
                .add(Rc::clone(&result), &[], &[]);
        } else {
            self.grad_f_cache
                .borrow_mut()
                .add_1dep(Rc::clone(&result), x.as_tagged());
        }
        result
    }

    /// Full-space constraint vector `g(x)` (length `n_full_g`), shared by
    /// `eval_c`/`eval_d` so the user `eval_g` runs once per iterate. On a
    /// cache hit no user evaluation occurs; on a failed eval the buffer is
    /// filled with NaN (so `theta_trial` goes non-finite and the line
    /// search backtracks), matching the per-subsystem paths.
    fn full_g(&self, x: &dyn Vector) -> Rc<Vec<Number>> {
        if let Some(v) = self.full_g_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        let n_full_g = self.adapter.borrow().classification().n_full_g as usize;
        let full_x = self.lift_x_to_full(x);
        let mut full_g = vec![0.0; n_full_g];
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_g(&full_x, true, &mut full_g)
        };
        if !ok {
            full_g.fill(f64::NAN);
        }
        let result = Rc::new(full_g);
        self.full_g_cache
            .borrow_mut()
            .add_1dep(Rc::clone(&result), x.as_tagged());
        result
    }

    /// Full-space Jacobian values (length `nnz_jac_g_full`, in the user's
    /// `eval_jac_g` order), shared by `eval_jac_c`/`eval_jac_d` so the user
    /// `eval_jac_g` runs once per iterate. Same NaN-on-failure contract as
    /// [`Self::full_g`].
    fn full_jac_g(&self, x: &dyn Vector) -> Rc<Vec<Number>> {
        if let Some(v) = self.full_jac_g_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        let mut full_vals = vec![0.0; self.nnz_jac_g_full as usize];
        let full_x = self.lift_x_to_full(x);
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_jac_g(
                Some(&full_x),
                true,
                SparsityRequest::Values {
                    values: &mut full_vals,
                },
            )
        };
        if !ok {
            full_vals.fill(f64::NAN);
        }
        let result = Rc::new(full_vals);
        self.full_jac_g_cache
            .borrow_mut()
            .add_1dep(Rc::clone(&result), x.as_tagged());
        result
    }

    fn eval_c_internal(&self, x: &dyn Vector) -> Rc<dyn Vector> {
        let cls = self.adapter.borrow().classification().clone();
        if cls.n_c == 0 {
            // Empty constraint vector — still cache so the tag is stable.
            if let Some(v) = self.c_cache.borrow().get_1dep(x.as_tagged()) {
                return v;
            }
            let v = self.c_space.make_new_dense();
            let result: Rc<dyn Vector> = Rc::new(v);
            self.c_cache
                .borrow_mut()
                .add_1dep(Rc::clone(&result), x.as_tagged());
            return result;
        }
        if let Some(v) = self.c_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        *self.c_evals.borrow_mut() += 1;
        // Shared full-space `g(x)` — computed once per iterate and reused
        // by `eval_d` (and vice versa). NaN-on-failure handled in `full_g`,
        // so `theta_trial` goes non-finite and the line search backtracks
        // (see `eval_f_internal`).
        let full_g = self.full_g(x);
        let mut c = self.c_space.make_new_dense();
        // c_i = g(g_idx) - c_rhs[i]  (since g_l == g_u for equalities,
        // upstream subtracts the bound to make it a residual). Matches
        // `OrigIpoptNLP::c` which calls `nlp_->Eval_c` after the adapter
        // subtracted the bound — TNLPAdapter doesn't subtract yet, so we
        // do it here. The RHS is the constant `g_l[g_idx]`, captured once
        // at construction (`self.c_rhs`, M17) — no per-iterate bounds
        // fetch or full-size scratch allocations in the line-search path.
        {
            let cv = c.values_mut();
            let cs = self.c_scale.borrow();
            for (i, &g_idx) in cls.c_map.iter().enumerate() {
                let raw = full_g[g_idx as usize] - self.c_rhs[i];
                cv[i] = match cs.as_ref() {
                    Some(v) => raw * v[i],
                    None => raw,
                };
            }
        }
        let result: Rc<dyn Vector> = Rc::new(c);
        self.c_cache
            .borrow_mut()
            .add_1dep(Rc::clone(&result), x.as_tagged());
        result
    }

    fn eval_d_internal(&self, x: &dyn Vector) -> Rc<dyn Vector> {
        let cls = self.adapter.borrow().classification().clone();
        if cls.n_d == 0 {
            if let Some(v) = self.d_cache.borrow().get_1dep(x.as_tagged()) {
                return v;
            }
            let v = self.d_space.make_new_dense();
            let result: Rc<dyn Vector> = Rc::new(v);
            self.d_cache
                .borrow_mut()
                .add_1dep(Rc::clone(&result), x.as_tagged());
            return result;
        }
        if let Some(v) = self.d_cache.borrow().get_1dep(x.as_tagged()) {
            return v;
        }
        *self.d_evals.borrow_mut() += 1;
        // Shared full-space `g(x)` — reused with `eval_c` (see `full_g`).
        let full_g = self.full_g(x);
        let mut d = self.d_space.make_new_dense();
        {
            let dv = d.values_mut();
            let ds = self.d_scale.borrow();
            for (i, &g_idx) in cls.d_map.iter().enumerate() {
                let raw = full_g[g_idx as usize];
                dv[i] = match ds.as_ref() {
                    Some(v) => raw * v[i],
                    None => raw,
                };
            }
        }
        let result: Rc<dyn Vector> = Rc::new(d);
        self.d_cache
            .borrow_mut()
            .add_1dep(Rc::clone(&result), x.as_tagged());
        result
    }

    fn eval_jac_c_internal(&self, x: &dyn Vector) -> Rc<dyn Matrix> {
        if self.const_deriv.jac_c
            && let Some(m) = self.jac_c_cache.borrow().get(&[], &[])
        {
            return m;
        }
        if let Some(m) = self.jac_c_cache.borrow().get_1dep(x.as_tagged()) {
            return m;
        }
        *self.jac_c_evals.borrow_mut() += 1;
        // Shared full-space Jacobian — computed once per iterate and reused
        // by `eval_jac_d` (and vice versa). NaN-on-failure handled in
        // `full_jac_g`.
        let full_vals = self.full_jac_g(x);
        let mut jac_c = GenTMatrix::new(Rc::clone(&self.jac_c_space));
        {
            let cs = self.c_scale.borrow();
            let irows = self.jac_c_space.irows().to_vec();
            let vs = jac_c.values_mut();
            for (k, &src) in self.jac_c_entry_in_g.iter().enumerate() {
                let raw = full_vals[src as usize];
                vs[k] = match cs.as_ref() {
                    // irows are 1-based.
                    Some(v) => raw * v[(irows[k] - 1) as usize],
                    None => raw,
                };
            }
        }
        let reuse = self.const_deriv.jac_c && all_finite(jac_c.values());
        let result: Rc<dyn Matrix> = Rc::new(jac_c);
        if reuse {
            self.jac_c_cache
                .borrow_mut()
                .add(Rc::clone(&result), &[], &[]);
        } else {
            self.jac_c_cache
                .borrow_mut()
                .add_1dep(Rc::clone(&result), x.as_tagged());
        }
        result
    }

    fn eval_jac_d_internal(&self, x: &dyn Vector) -> Rc<dyn Matrix> {
        if self.const_deriv.jac_d
            && let Some(m) = self.jac_d_cache.borrow().get(&[], &[])
        {
            return m;
        }
        if let Some(m) = self.jac_d_cache.borrow().get_1dep(x.as_tagged()) {
            return m;
        }
        *self.jac_d_evals.borrow_mut() += 1;
        // Shared full-space Jacobian — reused with `eval_jac_c` (see
        // `full_jac_g`).
        let full_vals = self.full_jac_g(x);
        let mut jac_d = GenTMatrix::new(Rc::clone(&self.jac_d_space));
        {
            let ds = self.d_scale.borrow();
            let irows = self.jac_d_space.irows().to_vec();
            let vs = jac_d.values_mut();
            for (k, &src) in self.jac_d_entry_in_g.iter().enumerate() {
                let raw = full_vals[src as usize];
                vs[k] = match ds.as_ref() {
                    Some(v) => raw * v[(irows[k] - 1) as usize],
                    None => raw,
                };
            }
        }
        let reuse = self.const_deriv.jac_d && all_finite(jac_d.values());
        let result: Rc<dyn Matrix> = Rc::new(jac_d);
        if reuse {
            self.jac_d_cache
                .borrow_mut()
                .add(Rc::clone(&result), &[], &[]);
        } else {
            self.jac_d_cache
                .borrow_mut()
                .add_1dep(Rc::clone(&result), x.as_tagged());
        }
        result
    }

    fn eval_h_internal(
        &self,
        x: &dyn Vector,
        obj_factor: Number,
        y_c: &dyn Vector,
        y_d: &dyn Vector,
    ) -> Rc<dyn SymMatrix> {
        // h_cache key: (x, y_c, y_d) tags + obj_factor scalar dep, as
        // upstream `IpOrigIpoptNLP.cpp:786`.
        //
        // A reused `∇²L` drops the three tags and keeps `obj_factor`:
        // the hint's premise is that every row is linear, so `λ` weights
        // nothing and `∇²L = σ·∇²f` — a function of `σ` alone. Upstream
        // drops `σ` too, which is safe only because the main algorithm
        // always passes 1.0; the restoration phase passes 0.0
        // (`resto_nlp.rs`), and keeping the scalar dependency costs one
        // float compare and makes that case right by construction rather
        // than by coincidence.
        if self.const_deriv.hessian
            && let Some(m) = self.h_cache.borrow().get(&[], &[obj_factor])
        {
            return m;
        }
        if let Some(m) = self.h_cache.borrow().get(
            &[x.as_tagged(), y_c.as_tagged(), y_d.as_tagged()],
            &[obj_factor],
        ) {
            return m;
        }
        *self.h_evals.borrow_mut() += 1;
        let Some(h_space) = self.h_space.as_ref() else {
            panic!(
                "OrigIpoptNlp::eval_h called but the TNLP did not provide \
                 eval_h sparsity. The L-BFGS path lands in Phase 8."
            );
        };
        let cls = self.adapter.borrow().classification().clone();
        let full_x = self.lift_x_to_full(x);
        // Upstream `IpOrigIpoptNLP.cpp:792-794` passes the user TNLP's
        // `eval_h` the multipliers in the user's unscaled-constraint
        // space, i.e. `lambda_user = c_scale * y_c` (and same for d).
        // The obj_factor is also scaled (`scaled_obj_factor = obj_scale
        // * obj_factor`). Together this gives the user-space Hessian
        // contribution that's already in the algorithm's scaled space
        // (no extra Hessian-side scaling because we don't scale x).
        let full_lambda = self.pack_lambda_for_user(y_c, y_d, &cls);
        let scaled_obj_factor = obj_factor * self.obj_scale_factor.get();

        // The user TNLP writes `nnz_h_lag_full` values; the kept
        // (var-x ⊗ var-x) subset has `h_space.nonzeros()` entries
        // selected via `h_entry_in_full`. They differ when fixed
        // variables drop entries.
        let mut full_vals = vec![0.0; self.nnz_h_lag_full as usize];
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.eval_h(
                Some(&full_x),
                true,
                scaled_obj_factor,
                Some(&full_lambda),
                true,
                SparsityRequest::Values {
                    values: &mut full_vals,
                },
            )
        };
        if !ok {
            full_vals.fill(f64::NAN);
        }
        let mut h = SymTMatrix::new(Rc::clone(h_space));
        let kept = h_space.nonzeros() as usize;
        let h_vals = h.values_mut();
        // `h_entry_in_full` always has length `kept` (identity when no
        // fixed-var filtering, sparse selection otherwise).
        debug_assert_eq!(kept, self.h_entry_in_full.len());
        for (k, &src) in self.h_entry_in_full.iter().enumerate() {
            h_vals[k] = full_vals[src as usize];
        }
        let reuse = self.const_deriv.hessian && all_finite(h.values());
        let result: Rc<dyn SymMatrix> = Rc::new(h);
        if reuse {
            self.h_cache
                .borrow_mut()
                .add(Rc::clone(&result), &[], &[obj_factor]);
        } else {
            self.h_cache.borrow_mut().add(
                Rc::clone(&result),
                &[x.as_tagged(), y_c.as_tagged(), y_d.as_tagged()],
                &[obj_factor],
            );
        }
        result
    }
}

// ---- helpers ----

/// Whether every value is finite. Guards the constant-derivative store:
/// a failed user evaluation NaN-fills its buffer, and a NaN written into
/// a cache entry that never expires would fail the rest of the solve
/// rather than the one trial point that caused it.
fn all_finite(v: &[Number]) -> bool {
    v.iter().all(|x| x.is_finite())
}

fn make_dense_from(
    space: &Rc<DenseVectorSpace>,
    mut f: impl FnMut(usize) -> Number,
) -> DenseVector {
    let mut v = space.make_new_dense();
    let dim = space.dim() as usize;
    if dim > 0 {
        let vs = v.values_mut();
        for (i, slot) in vs.iter_mut().enumerate().take(dim) {
            *slot = f(i);
        }
    }
    v
}

// -------------------- Trait impls --------------------

impl Nlp for OrigIpoptNlp {
    fn n(&self) -> Index {
        self.x_space.dim()
    }
    fn m_eq(&self) -> Index {
        self.c_space.dim()
    }
    fn m_ineq(&self) -> Index {
        self.d_space.dim()
    }

    fn eval_f(&mut self, x: &dyn Vector) -> Number {
        self.timed_eval(|t| &t.eval_obj, || self.eval_f_internal(x))
    }
    fn eval_grad_f(&mut self, x: &dyn Vector, g: &mut dyn Vector) {
        let result = self.timed_eval(|t| &t.eval_grad_obj, || self.eval_grad_f_internal(x));
        g.copy(&*result);
    }
    fn eval_c(&mut self, x: &dyn Vector, c: &mut dyn Vector) {
        let result = self.timed_eval(|t| &t.eval_constr, || self.eval_c_internal(x));
        c.copy(&*result);
    }
    fn eval_d(&mut self, x: &dyn Vector, d: &mut dyn Vector) {
        let result = self.timed_eval(|t| &t.eval_constr, || self.eval_d_internal(x));
        d.copy(&*result);
    }
    fn eval_jac_c(&mut self, x: &dyn Vector) -> Rc<dyn Matrix> {
        self.timed_eval(|t| &t.eval_constr_jac, || self.eval_jac_c_internal(x))
    }
    fn eval_jac_d(&mut self, x: &dyn Vector) -> Rc<dyn Matrix> {
        self.timed_eval(|t| &t.eval_constr_jac, || self.eval_jac_d_internal(x))
    }
    fn eval_h(
        &mut self,
        x: &dyn Vector,
        obj_factor: Number,
        y_c: &dyn Vector,
        y_d: &dyn Vector,
    ) -> Rc<dyn SymMatrix> {
        self.timed_eval(
            |t| &t.eval_lag_hess,
            || self.eval_h_internal(x, obj_factor, y_c, y_d),
        )
    }
}

impl IpoptNlp for OrigIpoptNlp {
    /// The Hessian's sparsity with zero values, built straight from
    /// `h_space` — no `eval_h`, no callback, no cache traffic. Falls
    /// back to a structurally empty block when the TNLP declared no
    /// Hessian sparsity at all (`nnz_h_lag == 0`), which is the same
    /// block the trait default would produce.
    fn uninitialized_h(&self) -> Rc<dyn SymMatrix> {
        match self.h_space.as_ref() {
            Some(space) => Rc::new(crate::ipopt_nlp::zeroed_sym_t(Rc::clone(space))),
            None => Rc::new(crate::ipopt_nlp::zeroed_sym_t(SymTMatrixSpace::new(
                self.x_space.dim(),
                Vec::new(),
                Vec::new(),
            ))),
        }
    }

    fn eval_counts(&self) -> [Index; 7] {
        [
            self.f_evals(),
            self.grad_f_evals(),
            self.c_evals(),
            self.d_evals(),
            self.jac_c_evals(),
            self.jac_d_evals(),
            self.h_evals(),
        ]
    }
    fn x_l(&self) -> &dyn Vector {
        &*self.x_l
    }
    fn x_u(&self) -> &dyn Vector {
        &*self.x_u
    }
    fn d_l(&self) -> &dyn Vector {
        &*self.d_l
    }
    fn d_u(&self) -> &dyn Vector {
        &*self.d_u
    }

    fn declared_d_bounds(&self) -> Option<(Vec<Number>, Vec<Number>)> {
        let mut dl = self.declared_d_l.borrow().clone()?;
        let mut du = self.declared_d_u.borrow().clone()?;
        // Return them in the live vectors' space: `apply_d_scale_to_bounds`
        // scaled `d_l`/`d_u` in place after the snapshot was taken, so the
        // same per-row factors apply here.
        if let Some(dd) = self.d_scale.borrow().as_ref() {
            let cls = self.adapter.borrow().classification().clone();
            for (i, slot) in dl.iter_mut().enumerate() {
                *slot *= dd[cls.d_l_map[i] as usize];
            }
            for (i, slot) in du.iter_mut().enumerate() {
                *slot *= dd[cls.d_u_map[i] as usize];
            }
        }
        Some((dl, du))
    }

    fn declared_box_violation(&self, x: &dyn Vector) -> Option<Number> {
        // Mirrors the `honor_original_bounds` projection just above
        // ([`Self::finalize_solution_x`]): same lift, same maps, same
        // declared bounds — it reports the distance instead of removing it.
        // Fixed variables are spliced in at their exact value, so only the
        // free block can have drifted past a bound.
        let x_l = self.declared_x_l.borrow();
        let x_u = self.declared_x_u.borrow();
        if x_l.is_none() && x_u.is_none() {
            return None;
        }
        let full = self.lift_x_to_full(x);
        let cls = self.adapter.borrow().classification().clone();
        let mut worst = 0.0_f64;
        if let Some(x_l) = x_l.as_ref() {
            for (i, &var_idx) in cls.x_l_map.iter().enumerate() {
                let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
                let viol = x_l[i] - full[full_idx];
                if viol.is_finite() && viol > worst {
                    worst = viol;
                }
            }
        }
        if let Some(x_u) = x_u.as_ref() {
            for (i, &var_idx) in cls.x_u_map.iter().enumerate() {
                let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
                let viol = full[full_idx] - x_u[i];
                if viol.is_finite() && viol > worst {
                    worst = viol;
                }
            }
        }
        Some(worst)
    }

    fn declared_x_bounds(&self) -> Option<(Vec<Number>, Vec<Number>)> {
        // No scaling to reapply, unlike `declared_d_bounds`: pounce models
        // objective and constraint scaling only, so nothing has touched the
        // variable box between the snapshot and now except the relaxation
        // this accessor exists to undo.
        let xl = self.declared_x_l.borrow().clone()?;
        let xu = self.declared_x_u.borrow().clone()?;
        Some((xl, xu))
    }

    fn declared_c_rhs(&self) -> Option<Vec<Number>> {
        // `c_rhs` is captured at construction from the user's `g_l` and never
        // touched afterwards — no relaxation applies to an equality row, so it
        // is already the declared value. Only the row scaling has to be
        // reapplied: `eval_c` emits `c_scale_i · (g_i(x) − b_i)`, so the RHS
        // must carry the same factor for the ratio to cancel it.
        let mut b = self.c_rhs.clone();
        if let Some(dc) = self.c_scale.borrow().as_ref() {
            for (i, slot) in b.iter_mut().enumerate() {
                *slot *= dc[i];
            }
        }
        Some(b)
    }

    fn px_l(&self) -> Rc<dyn Matrix> {
        Rc::clone(&self.px_l)
    }
    fn px_u(&self) -> Rc<dyn Matrix> {
        Rc::clone(&self.px_u)
    }
    fn pd_l(&self) -> Rc<dyn Matrix> {
        Rc::clone(&self.pd_l)
    }
    fn pd_u(&self) -> Rc<dyn Matrix> {
        Rc::clone(&self.pd_u)
    }

    /// Install moved bounds from the safe-slack mechanism. Mirrors
    /// `OrigIpoptNLP::AdjustVariableBounds` (`IpOrigIpoptNLP.cpp:990`):
    /// upstream simply swaps in the new bound vectors. We copy the values
    /// into the existing `Rc<DenseVector>` storage (falling back to a
    /// fresh allocation if the bound is somehow shared), which leaves the
    /// `Px_* / Pd_*` expansion matrices — keyed on the bound *spaces*,
    /// not values — untouched.
    fn adjust_variable_bounds(
        &mut self,
        new_x_l: &dyn Vector,
        new_x_u: &dyn Vector,
        new_d_l: &dyn Vector,
        new_d_u: &dyn Vector,
    ) {
        // The bound `Rc`s are uniquely owned (nothing clones them — same
        // invariant `relax_bounds` relies on), so `get_mut` always
        // succeeds and we copy the moved values into the existing storage.
        fn install(slot: &mut Rc<DenseVector>, new: &dyn Vector) {
            Rc::get_mut(slot)
                .expect("adjust_variable_bounds: bound vector is uniquely owned")
                .copy(new);
        }
        install(&mut self.x_l, new_x_l);
        install(&mut self.x_u, new_x_u);
        install(&mut self.d_l, new_d_l);
        install(&mut self.d_u, new_d_u);
    }

    fn obj_scaling_factor(&self) -> Number {
        self.obj_scale_factor.get()
    }

    fn computed_obj_scaling_factor(&self) -> Number {
        self.computed_obj_scale.get()
    }

    fn c_scale_vec(&self) -> Option<Vec<Number>> {
        self.c_scale.borrow().clone()
    }

    fn d_scale_vec(&self) -> Option<Vec<Number>> {
        self.d_scale.borrow().clone()
    }

    /// Project the underlying TNLP's `idx_names` metadata into the
    /// algorithm's split space. Variable names are gathered through the
    /// fixed-variable map (`x_not_fixed_map`), equality names through the
    /// c-block map (`c_map`), and inequality names through the d-block map
    /// (`d_map`) — exactly the permutations the adapter applied when it
    /// split the problem, so a residual at split index `k` is labeled with
    /// the equation the user actually wrote.
    ///
    /// Returns `None` when the TNLP exposes no names (e.g. presolve, which
    /// renumbers rows, declines `get_var_con_metadata`) so callers fall
    /// back to index labels rather than mislabeling permuted rows. This is
    /// the seam that turns "row 3" into `mass_balance` per Lee et al. (2024,
    /// <https://doi.org/10.69997/sct.147875>).
    fn split_space_names(&self) -> Option<SplitNames> {
        let a = self.adapter.borrow();
        let cls = a.classification();

        let mut var_meta = MetaData::default();
        let mut con_meta = MetaData::default();
        if !a
            .tnlp()
            .borrow_mut()
            .get_var_con_metadata(&mut var_meta, &mut con_meta)
        {
            return None;
        }

        // Full-space (original TNLP order) name pools. Either may be
        // absent — a model can name variables but not constraints, etc.
        let var_full = var_meta.strings.get(IDX_NAMES);
        let con_full = con_meta.strings.get(IDX_NAMES);
        if var_full.is_none() && con_full.is_none() {
            return None;
        }

        // Look a full-space name up safely; `None` for out-of-range or
        // empty entries so we degrade to an index label per slot.
        let pick = |pool: Option<&Vec<String>>, full_idx: Index| -> Option<String> {
            pool.and_then(|v| v.get(full_idx as usize))
                .filter(|s| !s.is_empty())
                .cloned()
        };

        let x_var = cls
            .x_not_fixed_map
            .iter()
            .map(|&full_idx| pick(var_full, full_idx))
            .collect();
        let eq = cls
            .c_map
            .iter()
            .map(|&full_idx| pick(con_full, full_idx))
            .collect();
        let ineq = cls
            .d_map
            .iter()
            .map(|&full_idx| pick(con_full, full_idx))
            .collect();

        let names = SplitNames { x_var, eq, ineq };
        names.any_present().then_some(names)
    }

    fn prepare_warm_start(&mut self) -> bool {
        let Some(snapshot) = self.fetch_warm_start_snapshot() else {
            return false;
        };
        *self.warm_start_snapshot.borrow_mut() = Some(snapshot);
        true
    }

    fn finish_warm_start(&mut self) {
        self.warm_start_snapshot.borrow_mut().take();
    }

    /// Populate `x` (length `n_x_var`) from the TNLP's starting point,
    /// compressed via `x_not_fixed_map`. Mirrors the `init_x` arm of
    /// upstream `IpOrigIpoptNLP::GetStartingPoint`.
    fn get_starting_x(&mut self, x: &mut dyn Vector) -> bool {
        let cls = self.adapter.borrow().classification().clone();
        if let Some(snapshot) = self.warm_start_snapshot.borrow().as_ref() {
            let Some(dx) = x.as_any_mut().downcast_mut::<DenseVector>() else {
                return false;
            };
            for (var_idx, &full_idx) in cls.x_not_fixed_map.iter().enumerate() {
                dx.values_mut()[var_idx] = snapshot.x[full_idx as usize];
            }
            return true;
        }
        let n_full_x = cls.n_full_x as usize;
        let n_full_g = cls.n_full_g as usize;
        let mut full_x = vec![0.0; n_full_x];
        let mut full_z_l = vec![0.0; n_full_x];
        let mut full_z_u = vec![0.0; n_full_x];
        let mut full_lambda = vec![0.0; n_full_g];
        let ok = {
            let a = self.adapter.borrow();
            let mut t = a.tnlp().borrow_mut();
            t.get_starting_point(StartingPoint {
                init_x: true,
                x: &mut full_x,
                init_z: false,
                z_l: &mut full_z_l,
                z_u: &mut full_z_u,
                init_lambda: false,
                lambda: &mut full_lambda,
            })
        };
        if !ok {
            return false;
        }
        let Some(dx) = x.as_any_mut().downcast_mut::<DenseVector>() else {
            return false;
        };
        let xs = dx.values_mut();
        for (var_idx, &full_idx) in cls.x_not_fixed_map.iter().enumerate() {
            xs[var_idx] = full_x[full_idx as usize];
        }
        true
    }

    fn get_starting_y(&mut self, y_c: &mut dyn Vector, y_d: &mut dyn Vector) -> bool {
        let Some(y_c) = y_c.as_any_mut().downcast_mut::<DenseVector>() else {
            return false;
        };
        let Some(y_d) = y_d.as_any_mut().downcast_mut::<DenseVector>() else {
            return false;
        };
        if let Some(snapshot) = self.warm_start_snapshot.borrow().as_ref() {
            let cls = self.adapter.borrow().classification().clone();
            let obj_scal = self.obj_scale_factor.get();
            let c_scale = self.c_scale.borrow();
            for (i, &g_idx) in cls.c_map.iter().enumerate() {
                let cs = c_scale.as_ref().map(|v| v[i]).unwrap_or(1.0);
                y_c.values_mut()[i] = snapshot.lambda[g_idx as usize] / cs * obj_scal;
            }
            let d_scale = self.d_scale.borrow();
            for (i, &g_idx) in cls.d_map.iter().enumerate() {
                let ds = d_scale.as_ref().map(|v| v[i]).unwrap_or(1.0);
                y_d.values_mut()[i] = snapshot.lambda[g_idx as usize] / ds * obj_scal;
            }
            return true;
        }
        let mut x = DenseVectorSpace::new(self.n()).make_new_dense();
        let mut z_l = DenseVectorSpace::new(self.x_l.dim()).make_new_dense();
        let mut z_u = DenseVectorSpace::new(self.x_u.dim()).make_new_dense();
        self.initialize_starting_point(
            &mut x, false, y_c, true, y_d, true, &mut z_l, false, &mut z_u, false,
        )
    }

    fn get_starting_z(
        &mut self,
        z_l: &mut dyn Vector,
        z_u: &mut dyn Vector,
        _v_l: &mut dyn Vector,
        _v_u: &mut dyn Vector,
    ) -> bool {
        // TNLP exposes only variable-bound multipliers.
        // Slack-bound v_l/v_u have no user-facing warm-start payload to forward.
        let Some(z_l) = z_l.as_any_mut().downcast_mut::<DenseVector>() else {
            return false;
        };
        let Some(z_u) = z_u.as_any_mut().downcast_mut::<DenseVector>() else {
            return false;
        };
        if let Some(snapshot) = self.warm_start_snapshot.borrow().as_ref() {
            let cls = self.adapter.borrow().classification().clone();
            let obj_scal = self.obj_scale_factor.get();
            for (i, slot) in z_l.values_mut().iter_mut().enumerate() {
                let var_idx = cls.x_l_map[i] as usize;
                let full_idx = cls.x_not_fixed_map[var_idx] as usize;
                *slot = snapshot.z_l[full_idx] * obj_scal;
            }
            for (i, slot) in z_u.values_mut().iter_mut().enumerate() {
                let var_idx = cls.x_u_map[i] as usize;
                let full_idx = cls.x_not_fixed_map[var_idx] as usize;
                *slot = snapshot.z_u[full_idx] * obj_scal;
            }
            return true;
        }
        let mut x = DenseVectorSpace::new(self.n()).make_new_dense();
        let mut y_c = DenseVectorSpace::new(self.m_eq()).make_new_dense();
        let mut y_d = DenseVectorSpace::new(self.m_ineq()).make_new_dense();
        self.initialize_starting_point(
            &mut x, false, &mut y_c, false, &mut y_d, false, z_l, true, z_u, true,
        )
    }

    fn lift_x_to_full(&self, x: &dyn Vector) -> Vec<Number> {
        OrigIpoptNlp::lift_x_to_full(self, x)
    }

    fn finalize_solution_x(&self, x: &dyn Vector) -> Vec<Number> {
        OrigIpoptNlp::finalize_solution_x(self, x)
    }

    fn n_full_x(&self) -> Index {
        self.adapter.borrow().classification().n_full_x
    }

    fn n_full_g(&self) -> Index {
        self.adapter.borrow().classification().n_full_g
    }

    fn pack_lambda_for_user(&self, y_c: &dyn Vector, y_d: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        OrigIpoptNlp::pack_lambda_for_user(self, y_c, y_d, &cls)
    }

    fn pack_g_for_user(&self, c: &dyn Vector, d: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let mut g = vec![0.0; cls.n_full_g as usize];
        if cls.n_c > 0 {
            let Some(dc) = c.as_any().downcast_ref::<DenseVector>() else {
                panic!("OrigIpoptNlp expects DenseVector for c");
            };
            let cs = self.c_scale.borrow();
            // Hoisted out of the loop on purpose. `expanded_values`
            // materializes a fresh `Vec` on every call (it has no cached
            // `expanded_values_` the way upstream's `DenseVector` does), so
            // calling it per row makes this scatter quadratic in the
            // constraint count: on 58k equality constraints that is 58k
            // allocations of 58k doubles, ~27 GB of memcpy and ~610 ms, once
            // per iteration for any caller that passes a non-NULL `g` to
            // `GetIpoptCurrentIterate`. The sibling `pack_z_*_for_user` /
            // `pack_lambda_for_user` scatters already lift it out. gh#698.
            let c_vals = dc.expanded_values();
            for (i, &g_idx) in cls.c_map.iter().enumerate() {
                let v = c_vals[i];
                g[g_idx as usize] = match cs.as_ref() {
                    Some(s) => v / s[i],
                    None => v,
                };
            }
        }
        if cls.n_d > 0 {
            let Some(dd) = d.as_any().downcast_ref::<DenseVector>() else {
                panic!("OrigIpoptNlp expects DenseVector for d");
            };
            let ds = self.d_scale.borrow();
            // Same hoist as the `c` block above.
            let d_vals = dd.expanded_values();
            for (i, &g_idx) in cls.d_map.iter().enumerate() {
                let v = d_vals[i];
                g[g_idx as usize] = match ds.as_ref() {
                    Some(s) => v / s[i],
                    None => v,
                };
            }
        }
        g
    }

    fn pack_z_l_for_user(&self, z_l: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let mut full = vec![0.0; cls.n_full_x as usize];
        if z_l.dim() == 0 {
            return full;
        }
        let Some(dz) = z_l.as_any().downcast_ref::<DenseVector>() else {
            panic!("OrigIpoptNlp expects DenseVector for z_l");
        };
        let vals = dz.expanded_values();
        for (k, &var_idx) in cls.x_l_map.iter().enumerate() {
            let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
            full[full_idx] = vals[k];
        }
        full
    }

    fn pack_z_u_for_user(&self, z_u: &dyn Vector) -> Vec<Number> {
        let cls = self.adapter.borrow().classification().clone();
        let mut full = vec![0.0; cls.n_full_x as usize];
        if z_u.dim() == 0 {
            return full;
        }
        let Some(dz) = z_u.as_any().downcast_ref::<DenseVector>() else {
            panic!("OrigIpoptNlp expects DenseVector for z_u");
        };
        let vals = dz.expanded_values();
        for (k, &var_idx) in cls.x_u_map.iter().enumerate() {
            let full_idx = cls.x_not_fixed_map[var_idx as usize] as usize;
            full[full_idx] = vals[k];
        }
        full
    }

    fn finalize_solution_lambda(&self, y_c: &dyn Vector, y_d: &dyn Vector) -> Vec<Number> {
        OrigIpoptNlp::finalize_solution_lambda(self, y_c, y_d)
    }

    fn finalize_solution_z_l(&self, z_l: &dyn Vector) -> Vec<Number> {
        OrigIpoptNlp::finalize_solution_z_l(self, z_l)
    }

    fn finalize_solution_z_u(&self, z_u: &dyn Vector) -> Vec<Number> {
        OrigIpoptNlp::finalize_solution_z_u(self, z_u)
    }

    fn variable_scaling(&self) -> Option<Vec<Number>> {
        // Forwarded, not stored: the substitution lives in the
        // `ScalingTnlp` the adapter wraps, and `TNLP::scaling_factors`
        // is the channel it reports through (gh#486). A transparent
        // decorator between the two forwards the inner answer, so one
        // hop off the adapter reaches whichever wrapper applied it.
        self.adapter.borrow().tnlp().borrow().scaling_factors()
    }

    fn full_x_to_var_x(&self, full_idx: Index) -> Option<Index> {
        let cls = self.adapter.borrow();
        let cls = cls.classification();
        let f = full_idx as usize;
        if f >= cls.full_to_var.len() {
            return None;
        }
        let v = cls.full_to_var[f];
        if v < 0 { None } else { Some(v) }
    }

    fn full_g_to_c_block(&self, full_idx: Index) -> Option<Index> {
        let cls = self.adapter.borrow();
        let cls = cls.classification();
        let f = full_idx as usize;
        if f >= cls.full_to_c.len() {
            return None;
        }
        let c = cls.full_to_c[f];
        if c < 0 { None } else { Some(c) }
    }

    fn full_g_to_d_block(&self, full_idx: Index) -> Option<Index> {
        let cls = self.adapter.borrow();
        let cls = cls.classification();
        let f = full_idx as usize;
        if f >= cls.full_to_d.len() {
            return None;
        }
        let d = cls.full_to_d[f];
        if d < 0 { None } else { Some(d) }
    }

    fn var_x_to_full_x(&self, var_idx: Index) -> Index {
        let cls = self.adapter.borrow();
        let cls = cls.classification();
        cls.x_not_fixed_map[var_idx as usize]
    }
}

// -------------------- Tests --------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::tnlp::{
        BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest,
        StartingPoint, TNLP,
    };

    /// HS071: min x[0]*x[3]*(x[0]+x[1]+x[2]) + x[2]
    /// s.t.   x[0]*x[1]*x[2]*x[3] >= 25                (inequality)
    ///        x[0]^2 + x[1]^2 + x[2]^2 + x[3]^2 == 40  (equality)
    ///        1 <= x[i] <= 5
    #[derive(Default)]
    struct Hs071 {
        eval_f_calls: usize,
        eval_grad_f_calls: usize,
        eval_g_calls: usize,
        eval_jac_g_value_calls: usize,
        eval_h_value_calls: usize,
        get_bounds_info_calls: usize,
        get_starting_point_calls: usize,
    }

    impl TNLP for Hs071 {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 4,
                m: 2,
                nnz_jac_g: 8,
                nnz_h_lag: 10,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            self.get_bounds_info_calls += 1;
            b.x_l.copy_from_slice(&[1.0; 4]);
            b.x_u.copy_from_slice(&[5.0; 4]);
            // Constraint 0: 25 <= g0 (inequality, finite lower only)
            // Constraint 1: g1 == 40                (equality)
            b.g_l.copy_from_slice(&[25.0, 40.0]);
            b.g_u.copy_from_slice(&[2.0e19, 40.0]);
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            self.get_starting_point_calls += 1;
            sp.x.copy_from_slice(&[1.0, 5.0, 5.0, 1.0]);
            if sp.init_z {
                sp.z_l.copy_from_slice(&[1.0, 2.0, 3.0, 4.0]);
                sp.z_u.copy_from_slice(&[5.0, 6.0, 7.0, 8.0]);
            }
            if sp.init_lambda {
                sp.lambda.copy_from_slice(&[11.0, 13.0]);
            }
            true
        }
        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
            self.eval_f_calls += 1;
            Some(x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2])
        }
        fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            self.eval_grad_f_calls += 1;
            // df/dx0 = x3*(2x0 + x1 + x2)
            // df/dx1 = x0*x3
            // df/dx2 = x0*x3 + 1
            // df/dx3 = x0*(x0 + x1 + x2)
            g[0] = x[3] * (2.0 * x[0] + x[1] + x[2]);
            g[1] = x[0] * x[3];
            g[2] = x[0] * x[3] + 1.0;
            g[3] = x[0] * (x[0] + x[1] + x[2]);
            true
        }
        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
            self.eval_g_calls += 1;
            // g0 = x0*x1*x2*x3 (>=25)
            // g1 = x0^2 + x1^2 + x2^2 + x3^2 (==40)
            g[0] = x[0] * x[1] * x[2] * x[3];
            g[1] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
            true
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            _new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    // Dense 2x4: row-major (g0 over x0..x3, then g1 over x0..x3).
                    irow.copy_from_slice(&[0, 0, 0, 0, 1, 1, 1, 1]);
                    jcol.copy_from_slice(&[0, 1, 2, 3, 0, 1, 2, 3]);
                }
                SparsityRequest::Values { values } => {
                    self.eval_jac_g_value_calls += 1;
                    let x = x.expect("eval_jac_g(Values) without x");
                    // d g0 / d x_j
                    values[0] = x[1] * x[2] * x[3];
                    values[1] = x[0] * x[2] * x[3];
                    values[2] = x[0] * x[1] * x[3];
                    values[3] = x[0] * x[1] * x[2];
                    // d g1 / d x_j
                    values[4] = 2.0 * x[0];
                    values[5] = 2.0 * x[1];
                    values[6] = 2.0 * x[2];
                    values[7] = 2.0 * x[3];
                }
            }
            true
        }
        fn eval_h(
            &mut self,
            x: Option<&[Number]>,
            _new_x: bool,
            obj_factor: Number,
            lambda: Option<&[Number]>,
            _new_lambda: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            // Dense lower triangle of 4x4 = 10 entries:
            // (0,0) (1,0) (1,1) (2,0) (2,1) (2,2) (3,0) (3,1) (3,2) (3,3)
            match mode {
                SparsityRequest::Structure { irow, jcol } => {
                    irow.copy_from_slice(&[0, 1, 1, 2, 2, 2, 3, 3, 3, 3]);
                    jcol.copy_from_slice(&[0, 0, 1, 0, 1, 2, 0, 1, 2, 3]);
                }
                SparsityRequest::Values { values } => {
                    self.eval_h_value_calls += 1;
                    let x = x.expect("eval_h(Values) without x");
                    let lam = lambda.expect("eval_h(Values) without lambda");
                    let of = obj_factor;
                    // Hessian of objective:
                    //   d2f/dx0^2 = 2*x3
                    //   d2f/dx0dx1 = x3,  d2f/dx0dx2 = x3,
                    //   d2f/dx0dx3 = 2*x0+x1+x2
                    //   d2f/dx1dx3 = x0,  d2f/dx2dx3 = x0
                    // Hessian of g0 = x0*x1*x2*x3:
                    //   d2/dx0dx1 = x2*x3, d2/dx0dx2 = x1*x3, d2/dx0dx3 = x1*x2
                    //   d2/dx1dx2 = x0*x3, d2/dx1dx3 = x0*x2, d2/dx2dx3 = x0*x1
                    // Hessian of g1 = sum x_i^2: 2*I.
                    let l0 = lam[0];
                    let l1 = lam[1];
                    values[0] = of * (2.0 * x[3]) + l1 * 2.0; // (0,0)
                    values[1] = of * x[3] + l0 * (x[2] * x[3]); // (1,0)
                    values[2] = l1 * 2.0; // (1,1)
                    values[3] = of * x[3] + l0 * (x[1] * x[3]); // (2,0)
                    values[4] = l0 * (x[0] * x[3]); // (2,1)
                    values[5] = l1 * 2.0; // (2,2)
                    values[6] = of * (2.0 * x[0] + x[1] + x[2]) + l0 * (x[1] * x[2]); // (3,0)
                    values[7] = of * x[0] + l0 * (x[0] * x[2]); // (3,1)
                    values[8] = of * x[0] + l0 * (x[0] * x[1]); // (3,2)
                    values[9] = l1 * 2.0; // (3,3)
                }
            }
            true
        }
        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
    }

    fn build_orig_nlp() -> (Rc<RefCell<TNLPAdapter>>, OrigIpoptNlp) {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Hs071::default()));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        (adapter, nlp)
    }

    fn dense_x(values: &[Number], space: &Rc<DenseVectorSpace>) -> DenseVector {
        let mut v = space.make_new_dense();
        v.values_mut().copy_from_slice(values);
        v
    }

    #[test]
    fn dimensions_match_classification() {
        let (_, nlp) = build_orig_nlp();
        // HS071: 4 vars (none fixed), 1 equality, 1 inequality.
        assert_eq!(nlp.n(), 4);
        assert_eq!(nlp.m_eq(), 1);
        assert_eq!(nlp.m_ineq(), 1);
        // 4 entries of jac_g go to c-row (g1), 4 go to d-row (g0).
        assert_eq!(nlp.jac_c_space().nonzeros(), 4);
        assert_eq!(nlp.jac_d_space().nonzeros(), 4);
        // Hessian sparsity comes through.
        assert_eq!(nlp.h_space().unwrap().nonzeros(), 10);
        // Bounds: all 4 x's bounded both sides; 1 ineq with finite lower only.
        assert_eq!(nlp.x_l().dim(), 4);
        assert_eq!(nlp.x_u().dim(), 4);
        assert_eq!(nlp.d_l().dim(), 1);
        assert_eq!(nlp.d_u().dim(), 0);
    }

    #[test]
    fn eval_f_at_starting_point() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        // f = 1*1*(1+5+5) + 5 = 11 + 5 = 16
        assert_eq!(nlp.eval_f(&x), 16.0);
        assert_eq!(nlp.f_evals(), 1);
    }

    #[test]
    fn eval_grad_f_at_starting_point() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut g = nlp.x_space().make_new_dense();
        nlp.eval_grad_f(&x, &mut g);
        // df/dx0 = 1*(2 + 5 + 5) = 12
        // df/dx1 = 1*1 = 1
        // df/dx2 = 1*1 + 1 = 2
        // df/dx3 = 1*(1 + 5 + 5) = 11
        assert_eq!(g.values(), &[12.0, 1.0, 2.0, 11.0]);
        assert_eq!(nlp.grad_f_evals(), 1);
    }

    #[test]
    fn eval_c_returns_equality_residual() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        // g1 = 1 + 25 + 25 + 1 = 52; residual = 52 - 40 = 12.
        assert_eq!(c.values(), &[12.0]);
        assert_eq!(nlp.c_evals(), 1);
    }

    #[test]
    fn eval_d_returns_inequality_value_unshifted() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut d = nlp.d_space().make_new_dense();
        nlp.eval_d(&x, &mut d);
        // g0 = 1*5*5*1 = 25.
        assert_eq!(d.values(), &[25.0]);
        assert_eq!(nlp.d_evals(), 1);
    }

    #[test]
    fn cache_returns_without_re_eval() {
        let (_, mut nlp) = build_orig_nlp();
        let mut x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let f1 = nlp.eval_f(&x);
        let f2 = nlp.eval_f(&x);
        assert_eq!(f1, f2);
        assert_eq!(nlp.f_evals(), 1, "second call must be served from cache");
        // Bumping x's tag (i.e. mutating it) should invalidate the cache.
        x.values_mut()[0] = 1.0; // values_mut bumps the cache.
        let _ = nlp.eval_f(&x);
        assert_eq!(nlp.f_evals(), 2);
    }

    #[test]
    fn jac_c_picks_only_equality_rows() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let m = nlp.eval_jac_c(&x);
        let g = m
            .as_any()
            .downcast_ref::<GenTMatrix>()
            .expect("jac_c is a GenTMatrix");
        // Equality is g1: dg1/dxj = 2*x_j.
        assert_eq!(g.values(), &[2.0, 10.0, 10.0, 2.0]);
        // 1-based row should all be 1 (the single equality row).
        assert_eq!(g.irows(), &[1, 1, 1, 1]);
        assert_eq!(g.jcols(), &[1, 2, 3, 4]);
    }

    #[test]
    fn jac_d_picks_only_inequality_rows() {
        let (_, mut nlp) = build_orig_nlp();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let m = nlp.eval_jac_d(&x);
        let g = m
            .as_any()
            .downcast_ref::<GenTMatrix>()
            .expect("jac_d is a GenTMatrix");
        // Inequality is g0: d/dxj of x0*x1*x2*x3 at (1,5,5,1).
        // d/dx0 = 5*5*1 = 25, d/dx1 = 1*5*1 = 5, d/dx2 = 1*5*1 = 5, d/dx3 = 1*5*5 = 25.
        assert_eq!(g.values(), &[25.0, 5.0, 5.0, 25.0]);
    }

    /// Build an `OrigIpoptNlp` over `Hs071` while retaining a typed handle
    /// to the underlying TNLP, so a test can read its `eval_g_calls` /
    /// `eval_jac_g_value_calls` counters (the adapter only exposes a
    /// `dyn TNLP`). Both `Rc`s alias the same allocation.
    fn build_orig_nlp_counting() -> (Rc<RefCell<Hs071>>, OrigIpoptNlp) {
        let concrete = Rc::new(RefCell::new(Hs071::default()));
        let tnlp: Rc<RefCell<dyn TNLP>> = concrete.clone();
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        (concrete, nlp)
    }

    #[test]
    fn eval_c_and_eval_d_share_one_eval_g_per_iterate() {
        // Code review 2026-06 item M16: `eval_c` and `eval_d` must slice
        // their rows out of ONE shared `g(x)`, not call the user `eval_g`
        // twice. Before the fix this asserted 2.
        let (tnlp, mut nlp) = build_orig_nlp_counting();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        let mut d = nlp.d_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        nlp.eval_d(&x, &mut d);
        assert_eq!(
            tnlp.borrow().eval_g_calls,
            1,
            "eval_c + eval_d at one iterate must share a single user eval_g"
        );
        // Per-subsystem counters still report one c and one d evaluation.
        assert_eq!(nlp.c_evals(), 1);
        assert_eq!(nlp.d_evals(), 1);
        // Values stay correct: c = g1 - 40 = 52 - 40 = 12, d = g0 = 25.
        assert_eq!(c.values(), &[12.0]);
        assert_eq!(d.values(), &[25.0]);

        // A genuinely new iterate (x mutated → tag bumped) costs exactly
        // one more eval_g shared across both subsystems.
        let mut x2 = x;
        x2.values_mut()[0] = 2.0;
        nlp.eval_c(&x2, &mut c);
        nlp.eval_d(&x2, &mut d);
        assert_eq!(
            tnlp.borrow().eval_g_calls,
            2,
            "a new iterate triggers exactly one more shared eval_g"
        );
    }

    #[test]
    fn eval_c_does_not_refetch_bounds_per_iterate() {
        // Code review 2026-06 item M17: the constant equality RHS is the
        // bound `g_l == g_u`, captured once at construction. `eval_c` must
        // NOT call the user's `get_bounds_info` on every (cache-missing)
        // iterate just to subtract that RHS. Before the fix each fresh
        // iterate re-fetched all bounds (and allocated four full-size
        // scratch vectors); this asserted the call count climbed with the
        // iterate count.
        let (tnlp, mut nlp) = build_orig_nlp_counting();
        // Construction fetches the bounds (once for classification, once in
        // `OrigIpoptNlp::new`). Snapshot whatever that baseline is.
        let baseline = tnlp.borrow().get_bounds_info_calls;

        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        // RHS is correct: c = g1 - 40 = (1+25+25+1) - 40 = 12.
        assert_eq!(c.values(), &[12.0]);

        // Several genuinely new iterates, each a cache miss.
        let mut x2 = x;
        for k in 0..5 {
            x2.values_mut()[0] = 2.0 + k as Number;
            nlp.eval_c(&x2, &mut c);
        }

        assert_eq!(
            tnlp.borrow().get_bounds_info_calls,
            baseline,
            "eval_c must reuse the captured c_rhs, not re-fetch bounds per iterate"
        );
    }

    #[test]
    fn eval_jac_c_and_eval_jac_d_share_one_eval_jac_g_per_iterate() {
        // Code review 2026-06 item M16: the full Jacobian is evaluated once
        // per iterate and sliced into jac_c / jac_d. Before the fix this
        // asserted 2.
        let (tnlp, mut nlp) = build_orig_nlp_counting();
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let _ = nlp.eval_jac_c(&x);
        let _ = nlp.eval_jac_d(&x);
        assert_eq!(
            tnlp.borrow().eval_jac_g_value_calls,
            1,
            "eval_jac_c + eval_jac_d at one iterate must share a single eval_jac_g"
        );
        assert_eq!(nlp.jac_c_evals(), 1);
        assert_eq!(nlp.jac_d_evals(), 1);
    }

    #[test]
    fn starting_point_is_compressed_into_x_var() {
        let (_, mut nlp) = build_orig_nlp();
        let mut x = nlp.x_space().make_new_dense();
        let mut yc = nlp.c_space().make_new_dense();
        let mut yd = nlp.d_space().make_new_dense();
        let mut zl = nlp.x_l_space().make_new_dense();
        let mut zu = nlp.x_u_space().make_new_dense();
        let ok = nlp.initialize_starting_point(
            &mut x, true, &mut yc, false, &mut yd, false, &mut zl, false, &mut zu, false,
        );
        assert!(ok);
        assert_eq!(x.values(), &[1.0, 5.0, 5.0, 1.0]);
    }

    #[test]
    fn warm_start_duals_are_forwarded_into_algorithm_vectors() {
        let (_, mut nlp) = build_orig_nlp();
        let mut y_c = nlp.c_space().make_new_dense();
        let mut y_d = nlp.d_space().make_new_dense();
        assert!(nlp.get_starting_y(&mut y_c, &mut y_d));
        assert_eq!(y_c.values(), &[13.0], "equality multiplier g1");
        assert_eq!(y_d.values(), &[11.0], "inequality multiplier g0");

        let mut z_l = nlp.x_l_space().make_new_dense();
        let mut z_u = nlp.x_u_space().make_new_dense();
        let mut v_l = nlp.d_l_space().make_new_dense();
        let mut v_u = nlp.d_u_space().make_new_dense();
        assert!(nlp.get_starting_z(&mut z_l, &mut z_u, &mut v_l, &mut v_u));
        assert_eq!(z_l.values(), &[1.0, 2.0, 3.0, 4.0]);
        assert_eq!(z_u.values(), &[5.0, 6.0, 7.0, 8.0]);
    }

    #[test]
    fn warm_start_prefetches_one_tnlp_snapshot_for_x_y_and_z() {
        let (tnlp, mut nlp) = build_orig_nlp_counting();
        assert!(nlp.prepare_warm_start());

        let mut x = nlp.x_space().make_new_dense();
        let mut y_c = nlp.c_space().make_new_dense();
        let mut y_d = nlp.d_space().make_new_dense();
        let mut z_l = nlp.x_l_space().make_new_dense();
        let mut z_u = nlp.x_u_space().make_new_dense();
        let mut v_l = nlp.d_l_space().make_new_dense();
        let mut v_u = nlp.d_u_space().make_new_dense();
        assert!(nlp.get_starting_x(&mut x));
        assert!(nlp.get_starting_y(&mut y_c, &mut y_d));
        assert!(nlp.get_starting_z(&mut z_l, &mut z_u, &mut v_l, &mut v_u));

        assert_eq!(tnlp.borrow().get_starting_point_calls, 1);
        assert_eq!(x.values(), &[1.0, 5.0, 5.0, 1.0]);
        assert_eq!(y_c.values(), &[13.0]);
        assert_eq!(y_d.values(), &[11.0]);
        assert_eq!(z_l.values(), &[1.0, 2.0, 3.0, 4.0]);
        assert_eq!(z_u.values(), &[5.0, 6.0, 7.0, 8.0]);

        nlp.finish_warm_start();
        let mut x_after_init = nlp.x_space().make_new_dense();
        assert!(nlp.get_starting_x(&mut x_after_init));
        assert_eq!(
            tnlp.borrow().get_starting_point_calls,
            2,
            "the snapshot must not affect later starting-point requests"
        );
    }

    /// Two-variable TNLP with `x[0]` fixed at 7.0 (`x_l == x_u`) and
    /// one equality on `x[1]`. Exercises the index-mapping methods on
    /// the `IpoptNlp` trait that are used by `pounce_sens` to support
    /// `.nl` files with fixed variables.
    struct OneFixedOneFree;
    impl TNLP for OneFixedOneFree {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                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] = 7.0;
            b.x_u[0] = 7.0; // fixed
            b.x_l[1] = -1.0e19;
            b.x_u[1] = 1.0e19;
            b.g_l[0] = 0.0;
            b.g_u[0] = 0.0; // equality
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 7.0;
            sp.x[1] = 0.5;
            true
        }
        fn eval_f(&mut self, x: &[Number], _: bool) -> Option<Number> {
            Some(x[1])
        }
        fn eval_grad_f(&mut self, _: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            g[1] = 1.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = x[1];
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 1;
                }
                SparsityRequest::Values { values } => values[0] = 1.0,
            }
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    /// Three-variable TNLP with `x[2]` fixed at 2.0, built so the fixed
    /// variable's reach is *different* for every form:
    ///
    /// ```text
    /// min  x0³                    ∇f, ∇²L vary; neither touches x2
    /// s.t. x0·x1 - 1 == 0         (equality)  ∇g0 varies, no x2
    ///      <row 1>       >= 0     (inequality) ∇g1 varies, touches x2
    /// ```
    ///
    /// Row 1 is the knob: `second_order = false` makes it `x1² + x2`, where
    /// x2 appears linearly and so reaches nothing but that row's Jacobian;
    /// `true` makes it `x1·x2`, where x2 is coupled to a free variable and
    /// so reaches ∇²L — and through it ∇f's and ∇g's character generally.
    /// Every declared proof is honest: all four forms really do vary.
    struct FixedVarReachModel {
        second_order: bool,
    }
    impl TNLP for FixedVarReachModel {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 3,
                m: 2,
                nnz_jac_g: 4,
                nnz_h_lag: 3,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l.copy_from_slice(&[-1.0e19, -1.0e19, 2.0]);
            b.x_u.copy_from_slice(&[1.0e19, 1.0e19, 2.0]); // x[2] fixed
            b.g_l.copy_from_slice(&[0.0, 0.0]);
            b.g_u.copy_from_slice(&[0.0, 1.0e19]); // g0 equality, g1 inequality
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x.copy_from_slice(&[1.0, 1.0, 2.0]);
            true
        }
        fn eval_f(&mut self, x: &[Number], _: bool) -> Option<Number> {
            Some(x[0] * x[0] * x[0])
        }
        fn eval_grad_f(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g.copy_from_slice(&[3.0 * x[0] * x[0], 0.0, 0.0]);
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = x[0] * x[1] - 1.0;
            g[1] = if self.second_order {
                x[1] * x[2]
            } else {
                x[1] * x[1] + x[2]
            };
            true
        }
        fn eval_jac_g(&mut self, x: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    // (0,0) (0,1) — row 0 is free-only; (1,1) (1,2) — row 1
                    // has a nonzero in the fixed column.
                    irow.copy_from_slice(&[0, 0, 1, 1]);
                    jcol.copy_from_slice(&[0, 1, 1, 2]);
                }
                SparsityRequest::Values { values } => {
                    let x = x.expect("values need x");
                    values[0] = x[1];
                    values[1] = x[0];
                    let (d1, d2) = if self.second_order {
                        (x[2], x[1])
                    } else {
                        (2.0 * x[1], 1.0)
                    };
                    values[2] = d1;
                    values[3] = d2;
                }
            }
            true
        }
        fn eval_h(
            &mut self,
            x: Option<&[Number]>,
            _: bool,
            sigma: Number,
            lambda: Option<&[Number]>,
            _: bool,
            m: SparsityRequest<'_>,
        ) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    // ∂²f/∂x0², ∂²g0/∂x1∂x0, and row 1's own second
                    // derivative — (1,1) when it is x1², (2,1) when it is
                    // x1·x2, which is the only entry in a fixed index.
                    irow.copy_from_slice(&[0, 1, if self.second_order { 2 } else { 1 }]);
                    jcol.copy_from_slice(&[0, 0, 1]);
                }
                SparsityRequest::Values { values } => {
                    let (x, l) = (x.expect("values need x"), lambda.expect("values need λ"));
                    values[0] = sigma * 6.0 * x[0];
                    values[1] = l[0];
                    values[2] = if self.second_order { l[1] } else { 2.0 * l[1] };
                }
            }
            true
        }
        fn derivative_proofs(&mut self) -> DerivativeProofs {
            DerivativeProofs {
                grad_f: DerivativeProof::Varying,
                hessian: DerivativeProof::Varying,
                jac: vec![DerivativeProof::Varying; 2],
            }
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    fn reach_proofs(second_order: bool) -> [DerivativeProof; 4] {
        let tnlp: Rc<RefCell<dyn TNLP>> =
            Rc::new(RefCell::new(FixedVarReachModel { second_order }));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(adapter, Rc::new(NoScaling)).unwrap();
        nlp.derivative_proofs()
    }

    /// Fixing a variable weakens the forms it reaches — and only those.
    ///
    /// `Varying` is what lets pounce *refuse* a false `*_constant=yes`;
    /// `Unknown` is honoured on trust. Weakening all four proofs because
    /// the model happened to fix one variable would hand a model back its
    /// own false hints on every form, which is the whole safety payload of
    /// gh#588 Q6. Row 1 is the only form the fixed x[2] reaches here, so
    /// it is the only one that may lose its refusal.
    #[test]
    fn only_the_forms_a_fixed_variable_reaches_lose_their_refusal() {
        use DerivativeProof::*;
        let [grad_f, hessian, jac_c, jac_d] = reach_proofs(false);
        assert_eq!(grad_f, Varying, "x[2] is nowhere in ∇f");
        assert_eq!(hessian, Varying, "x[2] appears linearly; ∇²L cannot see it");
        assert_eq!(jac_c, Varying, "row 0 has no nonzero in the fixed column");
        assert_eq!(jac_d, Unknown, "row 1 does, so fixing x[2] may flatten it");
    }

    /// The other direction: the weakening is not merely never applied.
    /// Couple the fixed variable to a free one and ∇²L — and with it ∇f,
    /// which is second-order in exactly the same sense — must give up its
    /// refusal too.
    #[test]
    fn a_second_order_coupling_to_a_fixed_variable_weakens_the_hessian() {
        use DerivativeProof::*;
        let [grad_f, hessian, jac_c, jac_d] = reach_proofs(true);
        assert_eq!(hessian, Unknown, "x1·x2 puts a fixed index in ∇²L");
        assert_eq!(grad_f, Unknown, "same test, and ∇f is reached the same way");
        assert_eq!(jac_c, Varying, "row 0 is still untouched by x[2]");
        assert_eq!(jac_d, Unknown);
    }

    #[test]
    fn ipopt_nlp_index_mapping_methods_handle_fixed_var() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneFixedOneFree));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();

        // Sanity: classification trimmed x[0] (fixed) from var-x space.
        assert_eq!(nlp.n_full_x(), 2);
        assert_eq!(nlp.n(), 1);

        // full_x_to_var_x: x[0] is fixed → None; x[1] → var idx 0.
        let nlp_dyn: &dyn crate::ipopt_nlp::IpoptNlp = &nlp;
        assert_eq!(nlp_dyn.full_x_to_var_x(0), None);
        assert_eq!(nlp_dyn.full_x_to_var_x(1), Some(0));

        // var_x_to_full_x: var 0 → full 1.
        assert_eq!(nlp_dyn.var_x_to_full_x(0), 1);

        // full_g_to_c_block: the one g is an equality → c-block 0.
        assert_eq!(nlp_dyn.full_g_to_c_block(0), Some(0));

        // full_g_to_d_block is its complement: an equality has no
        // d-block position, and an index past `m` is None in both
        // (gh#910).
        assert_eq!(nlp_dyn.full_g_to_d_block(0), None);
        assert_eq!(nlp_dyn.full_g_to_d_block(1), None);

        // lift_x_to_full inflates a compressed [v_0] back to [7.0, v_0].
        let mut x_var = nlp.x_space().make_new_dense();
        x_var.values_mut()[0] = 0.5;
        let lifted = nlp_dyn.lift_x_to_full(&x_var);
        assert_eq!(lifted, vec![7.0, 0.5]);
    }

    /// `OneFixedOneFree` plus `idx_names` metadata — used to check the
    /// split-space name projection threads names through the fixed-var
    /// and c/d-split permutations.
    struct NamedFixedOneFree;
    impl TNLP for NamedFixedOneFree {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            OneFixedOneFree.get_nlp_info()
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            OneFixedOneFree.get_bounds_info(b)
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            OneFixedOneFree.get_starting_point(sp)
        }
        fn eval_f(&mut self, x: &[Number], n: bool) -> Option<Number> {
            OneFixedOneFree.eval_f(x, n)
        }
        fn eval_grad_f(&mut self, x: &[Number], n: bool, g: &mut [Number]) -> bool {
            OneFixedOneFree.eval_grad_f(x, n, g)
        }
        fn eval_g(&mut self, x: &[Number], n: bool, g: &mut [Number]) -> bool {
            OneFixedOneFree.eval_g(x, n, g)
        }
        fn eval_jac_g(&mut self, x: Option<&[Number]>, n: bool, m: SparsityRequest<'_>) -> bool {
            OneFixedOneFree.eval_jac_g(x, n, m)
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
        fn get_var_con_metadata(&mut self, var: &mut MetaData, con: &mut MetaData) -> bool {
            var.strings.insert(
                IDX_NAMES.to_string(),
                vec!["fixed_x".to_string(), "free_x".to_string()],
            );
            con.strings
                .insert(IDX_NAMES.to_string(), vec!["balance".to_string()]);
            true
        }
    }

    #[test]
    fn split_space_names_threads_through_fixed_var_and_cd_split() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(NamedFixedOneFree));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();

        let names = nlp.split_space_names().expect("names present");
        // x[0] (fixed) dropped; var-x 0 is full-x 1 = "free_x".
        assert_eq!(names.x_var, vec![Some("free_x".to_string())]);
        // The single g is an equality → c-block 0 = "balance".
        assert_eq!(names.eq, vec![Some("balance".to_string())]);
        // No inequalities.
        assert!(names.ineq.is_empty());
        assert!(names.any_present());
    }

    #[test]
    fn split_space_names_none_when_tnlp_declines() {
        // OneFixedOneFree does not implement get_var_con_metadata.
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneFixedOneFree));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        assert!(nlp.split_space_names().is_none());
    }

    /// Regression: a TNLP with `x[0]` fixed and `nnz_h_lag = 1` whose
    /// only Hessian entry is (0,0). After fixed-var filtering kept = 0
    /// but `nnz_h_lag_full = 1`, which used to hit the broken
    /// `h_entry_in_full.is_empty()` fast path and panic in
    /// `copy_from_slice`.
    struct FixedOnlyHess;
    impl TNLP for FixedOnlyHess {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 1,
                nnz_jac_g: 1,
                nnz_h_lag: 1,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = 7.0;
            b.x_u[0] = 7.0; // fixed
            b.x_l[1] = -1.0e19;
            b.x_u[1] = 1.0e19;
            b.g_l[0] = 0.0;
            b.g_u[0] = 0.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 7.0;
            sp.x[1] = 0.5;
            true
        }
        fn eval_f(&mut self, x: &[Number], _: bool) -> Option<Number> {
            Some(0.5 * x[0] * x[0] + x[1])
        }
        fn eval_grad_f(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = x[0];
            g[1] = 1.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = x[1];
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 1;
                }
                SparsityRequest::Values { values } => values[0] = 1.0,
            }
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            obj_factor: Number,
            _: Option<&[Number]>,
            _: bool,
            m: SparsityRequest<'_>,
        ) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => values[0] = obj_factor,
            }
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    /// One variable with a single one-sided inequality whose Jacobian
    /// magnitude trips `nlp_scaling_max_gradient` (default 100): coeff
    /// 1000, bound `lo = 4e6`. After gradient-based scaling the
    /// `d_scale` for this row is `100/1000 = 0.1`, so the algorithm
    /// sees `d(x) = 0.1 * 1000 * x`. The bound must be scaled to
    /// `0.1 * 4e6 = 4e5` to match — otherwise the algorithm reads a
    /// 10x-too-large lower bound and reports phantom infeasibility
    /// (gh#54).
    struct OneIneqLargeOffset;
    impl TNLP for OneIneqLargeOffset {
        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] = -1.0e19;
            b.x_u[0] = 1.0e19;
            b.g_l[0] = 4.0e6;
            b.g_u[0] = 2.0e19;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 5000.0;
            true
        }
        fn eval_f(&mut self, _: &[Number], _: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 1000.0 * x[0];
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => values[0] = 1000.0,
            }
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    /// Equality twin of [`OneIneqLargeOffset`]: `1000·x == 4e6`, whose
    /// Jacobian magnitude likewise trips `nlp_scaling_max_gradient`, giving
    /// `c_scale = 100/1000 = 0.1`. Used to check that the declared equality
    /// RHS is reported in the same scaled space as `eval_c` (gh#390).
    struct OneEqLargeOffset;
    impl TNLP for OneEqLargeOffset {
        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] = -1.0e19;
            b.x_u[0] = 1.0e19;
            b.g_l[0] = 4.0e6;
            b.g_u[0] = 4.0e6;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 5000.0;
            true
        }
        fn eval_f(&mut self, _: &[Number], _: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 1000.0 * x[0];
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => values[0] = 1000.0,
            }
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    #[test]
    fn gradient_based_scaling_scales_d_l_and_d_u() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneIneqLargeOffset));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();

        // Pre-scaling: d_l carries the raw user bound 4e6.
        assert_eq!(nlp.d_l().dim(), 1);
        let pre = nlp
            .d_l()
            .as_any()
            .downcast_ref::<DenseVector>()
            .unwrap()
            .values()[0];
        assert_eq!(pre, 4.0e6);

        nlp.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            0.0,
            0.0,
        );

        // d_scale = 100 / 1000 = 0.1; bound must scale in step.
        let post = nlp
            .d_l()
            .as_any()
            .downcast_ref::<DenseVector>()
            .unwrap()
            .values()[0];
        assert!(
            (post - 4.0e5).abs() < 1e-9,
            "d_l should be scaled by d_scale=0.1; got {}",
            post
        );

        // And d(x) at the starting point must agree with the scaled
        // bound: d(5000) = 0.1 * 1000 * 5000 = 5e5 > 4e5, so feasible.
        let x = dense_x(&[5000.0], nlp.x_space());
        let mut d = nlp.d_space().make_new_dense();
        nlp.eval_d(&x, &mut d);
        assert!(
            (d.values()[0] - 5.0e5).abs() < 1e-6,
            "scaled d(x) mismatch; got {}",
            d.values()[0]
        );
        assert!(
            d.values()[0] >= post,
            "starting point must be feasible in scaled space"
        );
    }

    /// Same fixture as [`OneIneqLargeOffset`] but with a non-zero
    /// objective gradient (10), so we can verify that
    /// `nlp_scaling_obj_target_gradient` pins the scaled gradient
    /// ∞-norm exactly to the requested value (independent of the
    /// `max_gradient` cutoff).
    struct OneIneqWithObj;
    impl TNLP for OneIneqWithObj {
        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] = -1.0e19;
            b.x_u[0] = 1.0e19;
            b.g_l[0] = 4.0e6;
            b.g_u[0] = 2.0e19;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 5000.0;
            true
        }
        fn eval_f(&mut self, x: &[Number], _: bool) -> Option<Number> {
            Some(10.0 * x[0])
        }
        fn eval_grad_f(&mut self, _: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 10.0;
            true
        }
        fn eval_g(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 1000.0 * x[0];
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, m: SparsityRequest<'_>) -> bool {
            match m {
                SparsityRequest::Structure { irow, jcol } => {
                    irow[0] = 0;
                    jcol[0] = 0;
                }
                SparsityRequest::Values { values } => values[0] = 1000.0,
            }
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    #[test]
    fn obj_target_gradient_pins_obj_scale() {
        // grad_f = [10], so the default gradient-based path (max_grad=100,
        // 10 < cutoff) does NOT scale the objective: df = 1.
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneIneqWithObj));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        nlp.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            0.0, // no target → use cutoff path
            0.0,
        );
        assert!(
            (nlp.obj_scale_factor() - 1.0).abs() < 1e-12,
            "no-target path leaves df=1 when grad < cutoff; got {}",
            nlp.obj_scale_factor()
        );

        // With obj_target_gradient = 1.0 the scaled gradient ∞-norm
        // must be exactly 1, i.e. df = 1.0 / 10.0 = 0.1.
        let tnlp2: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneIneqWithObj));
        let adapter2 = Rc::new(RefCell::new(TNLPAdapter::new(tnlp2).unwrap()));
        let mut nlp2 = OrigIpoptNlp::new(Rc::clone(&adapter2), Rc::new(NoScaling)).unwrap();
        nlp2.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            1.0,
            0.0,
        );
        assert!(
            (nlp2.obj_scale_factor() - 0.1).abs() < 1e-12,
            "target_gradient=1, max_grad_f=10 → df=0.1; got {}",
            nlp2.obj_scale_factor()
        );
    }

    /// Regression (flosp2hm): gradient-based scaling must sample the
    /// objective gradient at the point the algorithm actually operates
    /// on — i.e. with fixed variables (`x_l == x_u`) lifted to their
    /// fixed value — not at the raw `x0` returned by `get_starting_point`.
    /// Here `x[1]` is fixed at 1000 but the starting point places it at 0,
    /// and the only free-variable gradient is `df/dx0 = x[1]`. Sampling at
    /// the raw `x0` gives `max_grad_f = 0` (df stays 1.0, no scaling);
    /// lifting `x[1]→1000` gives `max_grad_f = 1000`, so df = 100/1000 = 0.1.
    /// Pre-fix this left df=1 and stalled flosp2hm at max-iter.
    struct FixedVarShiftsObjGrad;
    impl TNLP for FixedVarShiftsObjGrad {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 2,
                m: 0,
                nnz_jac_g: 0,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = -1.0e19;
            b.x_u[0] = 1.0e19;
            b.x_l[1] = 1000.0;
            b.x_u[1] = 1000.0; // fixed at 1000
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 1.0;
            sp.x[1] = 0.0; // deliberately NOT the fixed value
            true
        }
        fn eval_f(&mut self, x: &[Number], _: bool) -> Option<Number> {
            Some(x[0] * x[1])
        }
        fn eval_grad_f(&mut self, x: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = x[1];
            g[1] = x[0];
            true
        }
        fn eval_g(&mut self, _: &[Number], _: bool, _: &mut [Number]) -> bool {
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, _: SparsityRequest<'_>) -> bool {
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    /// `min 0` over `x ∈ [-1, 2]`, one variable, no constraints. The
    /// objective is irrelevant: these tests drive the bound accessors
    /// directly rather than solving anything.
    struct BoxedVar;
    impl TNLP for BoxedVar {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Some(NlpInfo {
                n: 1,
                m: 0,
                nnz_jac_g: 0,
                nnz_h_lag: 0,
                index_style: IndexStyle::C,
            })
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            b.x_l[0] = -1.0;
            b.x_u[0] = 2.0;
            true
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            sp.x[0] = 0.0;
            true
        }
        fn eval_f(&mut self, _: &[Number], _: bool) -> Option<Number> {
            Some(0.0)
        }
        fn eval_grad_f(&mut self, _: &[Number], _: bool, g: &mut [Number]) -> bool {
            g[0] = 0.0;
            true
        }
        fn eval_g(&mut self, _: &[Number], _: bool, _: &mut [Number]) -> bool {
            true
        }
        fn eval_jac_g(&mut self, _: Option<&[Number]>, _: bool, _: SparsityRequest<'_>) -> bool {
            true
        }
        fn eval_h(
            &mut self,
            _: Option<&[Number]>,
            _: bool,
            _: Number,
            _: Option<&[Number]>,
            _: bool,
            _: SparsityRequest<'_>,
        ) -> bool {
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    fn boxed_var_nlp() -> OrigIpoptNlp {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(BoxedVar));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        OrigIpoptNlp::new(adapter, Rc::new(NoScaling)).unwrap()
    }

    fn at(x: Number) -> DenseVector {
        let space = DenseVectorSpace::new(1);
        let mut v = space.make_new_dense();
        v.set_values(&[x]);
        v
    }

    /// Without a snapshot the accessor **abstains**, and that is the whole
    /// reason the active-set SQP arm printed `Variable bound violation: nan`:
    /// it never called `relax_bounds`, because it applies no widening, so the
    /// declared box was never recorded and there was nothing to measure
    /// against.
    ///
    /// Pinned as its own case so the `nan` cannot come back by someone
    /// dropping the `snapshot_declared_bounds` call in
    /// `IpoptApplication::optimize_sqp_tnlp` — this states what that call is
    /// load-bearing for.
    #[test]
    fn the_box_violation_abstains_until_the_declared_bounds_are_snapshotted() {
        let nlp = boxed_var_nlp();
        assert!(
            nlp.declared_box_violation(&at(5.0)).is_none(),
            "no snapshot means no answer, not a fabricated zero"
        );
    }

    /// And after the snapshot it is a **measurement**: it tracks where the
    /// point is, on both sides of the box, and is zero only when the point is
    /// genuinely inside.
    ///
    /// This is what makes the SQP arm's zero trustworthy. That arm reports
    /// `0.0` on every fixture, which is the correct answer for it — it applies
    /// no `bound_relax_factor` widening, so its declared box IS the box it
    /// solves against, and the active-set QP keeps its iterate inside — but a
    /// hardcoded `0.0` would be indistinguishable from the outside. gh#900's
    /// subject was exactly that: "a zero that must not be fabricated". So the
    /// nonzero readings are pinned here, where the point can be placed by
    /// hand, rather than hoped for from a solve.
    #[test]
    fn the_box_violation_measures_the_distance_outside_the_declared_box() {
        let mut nlp = boxed_var_nlp();
        nlp.snapshot_declared_bounds();
        // Inside: zero, and it is measured rather than assumed.
        assert_eq!(nlp.declared_box_violation(&at(0.5)), Some(0.0));
        // Exactly on each bound is still inside.
        assert_eq!(nlp.declared_box_violation(&at(-1.0)), Some(0.0));
        assert_eq!(nlp.declared_box_violation(&at(2.0)), Some(0.0));
        // Past the upper bound, and past the lower one, by the distance.
        let over = nlp.declared_box_violation(&at(2.25)).unwrap();
        assert!((over - 0.25).abs() < 1e-15, "expected 0.25, got {over}");
        let under = nlp.declared_box_violation(&at(-1.5)).unwrap();
        assert!((under - 0.5).abs() < 1e-15, "expected 0.5, got {under}");
    }

    /// `relax_bounds` still takes the snapshot, so the interior-point arm is
    /// unaffected by the extraction — and it takes it BEFORE widening, which
    /// is what makes the row report the widening rather than zero.
    #[test]
    fn relaxing_the_bounds_still_snapshots_the_declared_ones_first() {
        let mut nlp = boxed_var_nlp();
        nlp.relax_bounds(1e-2, 1.0);
        // The live upper bound moved out to 2.02; the DECLARED one is still
        // 2.0, so a point at the widened bound reads 0.02 outside the model
        // the caller wrote. Reading the live bound here would report zero,
        // which is the defect the declared snapshot exists to prevent.
        let v = nlp.declared_box_violation(&at(2.02)).unwrap();
        assert!((v - 0.02).abs() < 1e-12, "expected 0.02, got {v}");
    }

    #[test]
    fn gradient_scaling_lifts_fixed_vars_to_their_value() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(FixedVarShiftsObjGrad));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();

        // Sanity: x[1] is fixed out of the var-x space, fixed value 1000.
        assert_eq!(nlp.n_full_x(), 2);
        assert_eq!(nlp.n(), 1);

        nlp.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            0.0,
            0.0,
        );

        // Lifted gradient ∞-norm over free vars is |df/dx0| = x[1] = 1000,
        // so df = 100/1000 = 0.1. Sampling at the raw x0 (x[1]=0) would
        // give 0 and leave df=1.0 (the pre-fix bug).
        assert!(
            (nlp.obj_scale_factor() - 0.1).abs() < 1e-12,
            "fixed var must be lifted before scaling; expected df=0.1, got {}",
            nlp.obj_scale_factor()
        );
    }

    #[test]
    fn constr_target_gradient_overrides_cutoff_and_clamp() {
        // Jacobian row max = 1000. Default gradient-based: cutoff 100
        // fires, dc = min(1, 100/1000) = 0.1. With
        // constr_target_gradient = 50 → dc = 50/1000 = 0.05 (no clamp
        // at 1, no cutoff check).
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneIneqLargeOffset));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        nlp.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            0.0,
            50.0,
        );
        let x = dense_x(&[5000.0], nlp.x_space());
        let mut d = nlp.d_space().make_new_dense();
        nlp.eval_d(&x, &mut d);
        // scaled d(x) = 0.05 * 1000 * 5000 = 2.5e5.
        assert!(
            (d.values()[0] - 2.5e5).abs() < 1e-6,
            "constr target=50 → dd=0.05; scaled d(5000)=2.5e5, got {}",
            d.values()[0]
        );
    }

    /// User-supplied TNLP that returns a per-constraint scaling vector
    /// via `get_scaling_parameters`. Constraint 0 is the equality (g1);
    /// constraint 1 is the inequality (g0). We reuse the HS071 fixture
    /// so the c/d split is well-defined.
    struct Hs071UserScaled;
    impl TNLP for Hs071UserScaled {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Hs071::default().get_nlp_info()
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            Hs071::default().get_bounds_info(b)
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            Hs071::default().get_starting_point(sp)
        }
        fn eval_f(&mut self, x: &[Number], new_x: bool) -> Option<Number> {
            Hs071::default().eval_f(x, new_x)
        }
        fn eval_grad_f(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_grad_f(x, new_x, g)
        }
        fn eval_g(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_g(x, new_x, g)
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_jac_g(x, new_x, mode)
        }
        fn eval_h(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            obj_factor: Number,
            lambda: Option<&[Number]>,
            new_lambda: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_h(x, new_x, obj_factor, lambda, new_lambda, mode)
        }
        fn get_scaling_parameters(&mut self, req: ScalingRequest<'_>) -> bool {
            *req.obj_scaling = 2.0;
            *req.use_x_scaling = false;
            *req.use_g_scaling = true;
            // HS071 g layout: g[0] = inequality, g[1] = equality.
            req.g_scaling[0] = 0.5;
            req.g_scaling[1] = 0.25;
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    #[test]
    fn user_scaling_dispatch_applies_obj_and_g_scaling() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Hs071UserScaled));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        nlp.determine_scaling_from_starting_point(
            ScalingMethod::UserScaling,
            100.0,
            1e-8,
            0.0,
            0.0,
        );

        // Objective scaling: 2.0 (no automatic floor needed since
        // user supplied a normal-sized factor).
        assert!(
            (nlp.obj_scale_factor() - 2.0).abs() < 1e-12,
            "user obj_scaling=2.0 should be installed; got {}",
            nlp.obj_scale_factor()
        );

        // Equality row (g1) gets g_scaling[1] = 0.25 → c-scaled
        // residual is 0.25× the unscaled one. Compute c at the
        // starting point.
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        // Unscaled: g1 = 1+25+25+1 = 52, residual = 52-40 = 12.
        // Scaled: 0.25 * 12 = 3.0.
        assert!(
            (c.values()[0] - 3.0).abs() < 1e-9,
            "user g_scaling=0.25 on equality → c=3.0; got {}",
            c.values()[0]
        );

        // Inequality row (g0) gets g_scaling[0] = 0.5 → d = 0.5 *
        // 1*5*5*1 = 12.5.
        let mut d = nlp.d_space().make_new_dense();
        nlp.eval_d(&x, &mut d);
        assert!(
            (d.values()[0] - 12.5).abs() < 1e-9,
            "user g_scaling=0.5 on inequality → d=12.5; got {}",
            d.values()[0]
        );

        // And d_l must have been brought along: the user lower bound
        // on g0 is 25 (HS071); scaled by 0.5 → 12.5.
        let post_d_l = nlp
            .d_l()
            .as_any()
            .downcast_ref::<DenseVector>()
            .unwrap()
            .values()[0];
        assert!(
            (post_d_l - 12.5).abs() < 1e-9,
            "d_l scaled in step: got {}",
            post_d_l
        );
    }

    /// TNLP whose `get_scaling_parameters` returns false — selecting
    /// `UserScaling` must fall back to no automatic scaling (matches
    /// upstream behavior).
    struct Hs071DeclinesScaling;
    impl TNLP for Hs071DeclinesScaling {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Hs071::default().get_nlp_info()
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            Hs071::default().get_bounds_info(b)
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            Hs071::default().get_starting_point(sp)
        }
        fn eval_f(&mut self, x: &[Number], new_x: bool) -> Option<Number> {
            Hs071::default().eval_f(x, new_x)
        }
        fn eval_grad_f(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_grad_f(x, new_x, g)
        }
        fn eval_g(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_g(x, new_x, g)
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_jac_g(x, new_x, mode)
        }
        fn eval_h(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            obj_factor: Number,
            lambda: Option<&[Number]>,
            new_lambda: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_h(x, new_x, obj_factor, lambda, new_lambda, mode)
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    #[test]
    fn user_scaling_falls_back_when_tnlp_declines() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Hs071DeclinesScaling));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        nlp.determine_scaling_from_starting_point(
            ScalingMethod::UserScaling,
            100.0,
            1e-8,
            0.0,
            0.0,
        );
        // No automatic scaling installed: obj_scale_factor = 1.0, c/d
        // unscaled.
        assert!((nlp.obj_scale_factor() - 1.0).abs() < 1e-12);
        let x = dense_x(&[1.0, 5.0, 5.0, 1.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        assert_eq!(c.values(), &[12.0], "unscaled equality residual");
    }

    /// HS071 whose `get_scaling_parameters` asks for per-variable
    /// factors — the channel `OrigIpoptNlp` cannot model (gh#483).
    struct Hs071XScaled(Vec<Number>);
    impl TNLP for Hs071XScaled {
        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
            Hs071::default().get_nlp_info()
        }
        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
            Hs071::default().get_bounds_info(b)
        }
        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
            Hs071::default().get_starting_point(sp)
        }
        fn eval_f(&mut self, x: &[Number], new_x: bool) -> Option<Number> {
            Hs071::default().eval_f(x, new_x)
        }
        fn eval_grad_f(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_grad_f(x, new_x, g)
        }
        fn eval_g(&mut self, x: &[Number], new_x: bool, g: &mut [Number]) -> bool {
            Hs071::default().eval_g(x, new_x, g)
        }
        fn eval_jac_g(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_jac_g(x, new_x, mode)
        }
        fn eval_h(
            &mut self,
            x: Option<&[Number]>,
            new_x: bool,
            obj_factor: Number,
            lambda: Option<&[Number]>,
            new_lambda: bool,
            mode: SparsityRequest<'_>,
        ) -> bool {
            Hs071::default().eval_h(x, new_x, obj_factor, lambda, new_lambda, mode)
        }
        fn get_scaling_parameters(&mut self, req: ScalingRequest<'_>) -> bool {
            *req.obj_scaling = 2.0;
            *req.use_x_scaling = true;
            req.x_scaling.copy_from_slice(&self.0);
            *req.use_g_scaling = false;
            true
        }
        fn finalize_solution(&mut self, _: Solution<'_>, _: &IpoptData, _: &IpoptCq) {}
    }

    fn user_x_scaling_run(factors: &[Number]) -> OrigIpoptNlp {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(Hs071XScaled(factors.to_vec())));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        nlp.determine_scaling_from_starting_point(
            ScalingMethod::UserScaling,
            100.0,
            1e-8,
            0.0,
            0.0,
        );
        nlp
    }

    /// gh#483: a per-variable scaling request pounce cannot honor is
    /// flagged for the driver instead of being discarded. Before this,
    /// `scale_user_supplied` ended in `let _ = use_x_scaling;` and the
    /// solve ran with the caller's variable scaling silently gone.
    #[test]
    fn user_x_scaling_request_is_flagged_not_discarded() {
        let nlp = user_x_scaling_run(&[1.0, 1e3, 1.0, 1.0]);
        assert!(
            nlp.user_x_scaling_rejected(),
            "a non-unit x_scaling must be refused, not dropped"
        );
        // The objective factor is still installed — the flag is the
        // driver's cue to abort, not a reason to skip the other axes.
        assert!((nlp.obj_scale_factor() - 2.0).abs() < 1e-12);
    }

    /// An all-ones request asks for nothing, so it is a genuine no-op
    /// and must not fail a solve that would otherwise run.
    #[test]
    fn unit_x_scaling_request_is_not_rejected() {
        let nlp = user_x_scaling_run(&[1.0, 1.0, 1.0, 1.0]);
        assert!(!nlp.user_x_scaling_rejected());
    }

    #[test]
    fn eval_h_with_all_entries_on_fixed_var_does_not_panic() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(FixedOnlyHess));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();

        // After filtering, the kept Hessian over var-x has 0 nonzeros,
        // while the user's full Hessian has 1.
        assert_eq!(nlp.h_space().unwrap().nonzeros(), 0);

        let x = dense_x(&[0.5], &nlp.x_space().clone());
        let yc = dense_x(&[0.0], &nlp.c_space().clone());
        let yd = nlp.d_space().make_new_dense();
        let h = nlp.eval_h(&x, 1.0, &yc, &yd);
        assert_eq!(h.n_rows(), 1);
    }

    #[test]
    fn relax_bounds_widens_uniquely_owned_bounds() {
        // Baseline: with uniquely-owned bound Rcs (the normal post-construction
        // state) relax_bounds loosens x_l downward and x_u upward.
        let (_adapter, mut nlp) = build_orig_nlp();
        let x_l_before = nlp.x_l.values().to_vec();
        let x_u_before = nlp.x_u.values().to_vec();
        nlp.relax_bounds(1e-2, 1.0);
        for (b, a) in x_l_before.iter().zip(nlp.x_l.values()) {
            assert!(a < b, "x_l should relax downward: {a} !< {b}");
        }
        for (b, a) in x_u_before.iter().zip(nlp.x_u.values()) {
            assert!(a > b, "x_u should relax upward: {a} !> {b}");
        }
    }

    /// #385 Step 6: inequality-row bounds relax by a *scale-relative* delta
    /// (`min(relax, cap) · |b|`), and the declared (pre-relax) bounds stay
    /// available for the scale-relative feasibility measure — the live vector
    /// alone cannot distinguish a declared `2e-12` bound from a relaxed zero.
    #[test]
    fn relax_bounds_is_scale_relative_on_d_and_snapshots_declared() {
        // HS071's inequality row is `x1*x2*x3*x4 >= 25`.
        let (_adapter, mut nlp) = build_orig_nlp();
        assert_eq!(nlp.d_l.values(), &[25.0]);
        assert_eq!(nlp.declared_d_bounds(), None, "no snapshot before relax");
        nlp.relax_bounds(1e-2, 1.0);
        // delta = min(1e-2, 1.0) * 25 = 0.25 — proportional to the bound, so
        // the same row written at any scaling relaxes to the same feasible
        // set. The upstream form `min(cap, relax*max(|b|,1))` coincides here;
        // where they differ (|b| < 1) the old absolute floor erased
        // down-scaled rows entirely (a 2e-12 bound relaxed by 1e-8).
        assert_eq!(nlp.d_l.values(), &[25.0 - 0.25]);
        let (dl, du) = nlp.declared_d_bounds().expect("snapshotted at relax");
        assert_eq!(dl, vec![25.0], "declared bound is the pre-relax value");
        assert!(du.is_empty() || du[0] >= 25.0); // HS071: no finite d upper
    }

    /// gh#612: the variable box the user declared survives the relaxation.
    ///
    /// Crossover pivots against these, not the live vector. HS071's box is
    /// `1 <= x_i <= 5`; a solution that sits exactly on `x = 1` is a full
    /// `delta` *inside* the relaxed `1 - delta`, so an activity test against
    /// the live bound reports the binding bound inactive — the one answer
    /// crossover exists to get right.
    #[test]
    fn declared_x_bounds_are_the_pre_relax_box() {
        let (_adapter, mut nlp) = build_orig_nlp();
        assert_eq!(nlp.declared_x_bounds(), None, "no snapshot before relax");
        let x_l_before = nlp.x_l.values().to_vec();
        let x_u_before = nlp.x_u.values().to_vec();
        nlp.relax_bounds(1e-2, 1.0);
        let (xl, xu) = nlp.declared_x_bounds().expect("snapshotted at relax");
        assert_eq!(xl, x_l_before, "declared lower box is the pre-relax value");
        assert_eq!(xu, x_u_before, "declared upper box is the pre-relax value");
        // And the live vectors did move, so the two are genuinely distinct
        // rather than the accessor happening to alias an unrelaxed bound.
        assert!(nlp.x_l.values()[0] < xl[0]);
        assert!(nlp.x_u.values()[0] > xu[0]);
    }

    /// gh#390: the equality RHS folded into `c(x) = 0` is plumbed back out, so
    /// the runtime feasibility measure has a magnitude to judge `|c_i|`
    /// against. Unlike an inequality bound it is never relaxed — an equality
    /// row has no bound to widen — so the captured value is already the
    /// declared one, at every point in the solve.
    #[test]
    fn declared_c_rhs_is_the_pre_fold_right_hand_side() {
        // HS071's equality row is `x1² + x2² + x3² + x4² == 40`.
        let (_adapter, mut nlp) = build_orig_nlp();
        assert_eq!(nlp.declared_c_rhs(), Some(vec![40.0]));
        nlp.relax_bounds(1e-2, 1.0);
        assert_eq!(
            nlp.declared_c_rhs(),
            Some(vec![40.0]),
            "bound relaxation must not reach the equality RHS"
        );
    }

    /// The RHS is reported in the same space as `eval_c`'s output, so the
    /// ratio `|c_i| / |b_i|` cancels the solver's own row scaling — the
    /// property that makes it a scale-*free* measure rather than one that
    /// merely moved which scale it depends on.
    #[test]
    fn declared_c_rhs_carries_the_row_scaling() {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(OneEqLargeOffset));
        let adapter = Rc::new(RefCell::new(TNLPAdapter::new(tnlp).unwrap()));
        let mut nlp = OrigIpoptNlp::new(Rc::clone(&adapter), Rc::new(NoScaling)).unwrap();
        assert_eq!(nlp.declared_c_rhs(), Some(vec![4.0e6]));

        nlp.determine_scaling_from_starting_point(
            ScalingMethod::GradientBased,
            100.0,
            1e-8,
            0.0,
            0.0,
        );
        // c_scale = 100 / 1000 = 0.1.
        let rhs = nlp.declared_c_rhs().unwrap();
        assert!(
            (rhs[0] - 4.0e5).abs() < 1e-9,
            "declared RHS should carry c_scale=0.1; got {}",
            rhs[0]
        );

        // At x = 5000: unscaled residual 5e6 - 4e6 = 1e6 over an unscaled RHS
        // of 4e6 is 0.25 — and the scaled pair reads the same 0.25.
        let x = dense_x(&[5000.0], nlp.x_space());
        let mut c = nlp.c_space().make_new_dense();
        nlp.eval_c(&x, &mut c);
        assert!((c.values()[0] / rhs[0] - 0.25).abs() < 1e-12);
    }

    #[test]
    #[should_panic(expected = "x_l is uniquely owned")]
    fn relax_bounds_panics_on_shared_bound_rc() {
        // Code review L33: a shared bound Rc used to make relax_bounds silently
        // skip the relaxation, leaving bounds tighter than bound_relax_factor
        // requires. The unique-ownership invariant is now enforced loudly,
        // matching adjust_variable_bounds' `expect`.
        let (_adapter, mut nlp) = build_orig_nlp();
        let _shared = Rc::clone(&nlp.x_l); // bump strong_count so get_mut fails
        nlp.relax_bounds(1e-2, 1.0);
    }
}