pounce-sensitivity 0.12.0

//! Post-solve activity classification (the covariance/information
//! roadmap's item 0, gh #362).
//!
//! Classifies every bounded variable and every finite-bounded inequality
//! row of a converged barrier solve into one of five statuses, keyed on
//! the ratio of barrier curvature to the model's own curvature:
//!
//! ```text
//! r = Σ / q,   Σ = z/s summed over the sides that exist,
//!              q = |H_ii|                        (variable)
//!                  |∇dⱼᵀ H ∇dⱼ| / ‖∇dⱼ‖⁴         (inequality row)
//! ```
//!
//! The row denominator carries the fourth power so that `r` is
//! invariant to rescaling the row: `d → c·d` sends `Σ → Σ/c²` while
//! the curvature along the unit normal is unchanged, and `‖∇d‖⁴`
//! restores the balance. Equivalently, the geometric barrier weight
//! `Σ‖∇d‖²` (distance to the surface is `d/‖∇d‖`, its conjugate
//! multiplier `v‖∇d‖`) is measured against the curvature along the
//! unit normal. Variable bounds are invariant as written. This also
//! absorbs the solver's own per-row `d_scale`.
//!
//! `H` is the exact Lagrangian Hessian, so constraint curvature
//! contributes to `q` alongside the objective's. For variables, `q`
//! reads the Hessian DIAGONAL only, so purely off-diagonal coupling is
//! invisible to it: `f = x₁x₂` with bounds on both variables reports
//! `unidentified` on every bound even though the bound directions have
//! well-defined curvature. Items 1-4 of the covariance roadmap inherit
//! these semantics where they consume the per-coordinate statuses;
//! their reduced-block classification is where coupling becomes
//! visible, folded into the reduced diagonal by elimination.
//!
//! # The diagonal is not the curvature that generates the multiplier
//!
//! The same DIAGONAL-only reading has a consequence sharper than a
//! missed `q`, and it is the one to know before building on a status
//! (gh#763). At a kink the multiplier is generated by the curvature
//! **reduced** along the coordinate — what is left after the other
//! free variables re-optimize — not by `H_ii`. Eliminating a free
//! partner `y` from `[[h, c], [c, m]]` leaves `h − c²/m`, and
//! `Σ = z/s` equals exactly that, so
//!
//! ```text
//! r = reduced / diagonal
//! ```
//!
//! which is `1` only where the coordinate is **decoupled**. Couple it
//! and a genuine kink drops out of the `[1e-1, 1e1]` band and reads
//! [`AMBIGUOUS`] — at any tolerance, because `r` there is
//! `μ`-independent, so re-solving tighter reports the same thing. On a
//! collocation model coupling between neighbouring coordinates is the
//! normal case, not a corner.
//!
//! So **[`AMBIGUOUS`] is not "probably not a kink"**, and the class
//! must not be used as a proxy for kink-ness. That inference is not
//! hypothetical: gh#756 made it and shipped a first-order wrong
//! derivative. [`reduced_activity`] is the accessor that answers the
//! question the class does not — one back-solve per coordinate,
//! normalizing by the reduced curvature — and on it the same kink
//! reads [`WEAKLY_ACTIVE`] at every coupling.
//!
//! The default stays the diagonal because the reduced normalizer is
//! the reciprocal diagonal of an *inverse*: there is no
//! diagonal-of-the-inverse shortcut, so classifying every bounded
//! variable that way is `n` back-solves, and at 62k variables that is
//! not a post-solve diagnostic any more. The refinement is on demand,
//! over the entries in question.
//!
//! The row path normalizes by a directional curvature rather than a
//! diagonal, and `∇dᵀH∇d/‖∇d‖²` is a genuine curvature along the
//! row's own gradient — strictly better than a bare `H_ii`, which is
//! why it was not the one gh#763 fixed. But it is not *reduced*
//! either: the other free coordinates still re-optimize. So a row's
//! `r` is `reduced/directional` by the same algebra, `1` only where
//! the row's direction is decoupled from the remaining free space,
//! and a coupled row kink reads [`AMBIGUOUS`] at any tolerance for
//! the same `μ`-independent reason (gh#804).
//! [`reduced_row_activity`] is the row half of the answer, one
//! back-solve per row: the row's own value IS a coordinate of the KKT
//! system — the slack the barrier acts on, tied to the model by
//! `dⱼ(x) = sⱼ` — so it is the same back-solve one block over.
//!
//! `r` is `O(μ)` when the bound is inactive, `O(1)` when weakly active
//! (slack and multiplier vanish together), and `O(1/μ)` when strongly
//! active, so one ratio separates the regimes at any `μ` where a fixed
//! threshold on the slack or the multiplier alone cannot: both are
//! `O(√μ)` at weak activity, so any constant tracks the solve rather
//! than the geometry.
//!
//! Everything read here is retained by the converged state the
//! backsolver already holds: the bound multipliers on the iterate, the
//! solver's own slacks, `Σ` through the backsolver's
//! `barrier_sigma_x` / `barrier_sigma_s` — `curr_sigma_x` /
//! `curr_sigma_s` unless the held iterate came from crossover, in which
//! case the declared-frame diagonal the factor is also built with
//! (gh#654) — the barrier parameter, and the exact Lagrangian Hessian,
//! so `H` is never recovered from the barrier-augmented factor.
//!
//! The report is indexed in **user space**: `var_*` arrays have the
//! user TNLP's full variable count and `row_*` arrays its full
//! constraint count. A variable removed internally by
//! `fixed_variable_treatment = make_parameter` (`lb == ub`, the
//! default) reports [`FIXED`] at its own user index, and an equality
//! constraint reports [`EQUALITY`], so user indices never shift.

use std::rc::Rc;

use pounce_common::types::{Index, Number};
use pounce_linalg::Matrix;
use pounce_linalg::dense_vector::{DenseVector, DenseVectorSpace};
use pounce_linalg::expansion_matrix::ExpansionMatrix;
use pounce_linalg::triplet::{GenTMatrix, SymTMatrix};

use crate::PdSensBacksolver;
use crate::backsolver::SensBacksolver;
use crate::vec_util::dense_to_vec;

// The status codes and the classification rule now live in
// `pounce-sens-core`, so the convex arm decides what a kink is with the same
// code rather than a parallel one. Re-exported here because
// `pounce_sensitivity::activity::WEAKLY_ACTIVE` (and its siblings) is the path
// `pounce-py` and four test files already use.
pub use pounce_sens_core::activity_kernel::{
    AMBIGUOUS, EQUALITY, FIXED, INACTIVE, STRONGLY_ACTIVE, UNBOUNDED, UNIDENTIFIED, WEAKLY_ACTIVE,
};
use pounce_sens_core::activity_kernel::{
    Entry, NOT_CLASSIFIED, classify_entry, off_path, sign_of, zero_gradient_row,
};

/// Per-variable and per-row classification of a converged solve.
///
/// All vectors are **user-space**: `var_*` have length `n_full_x` (the
/// user TNLP's `n`) and `row_*` length `n_full_g` (the user's `m`).
/// Entries with no finite bound hold [`UNBOUNDED`]; [`FIXED`]
/// variables and [`EQUALITY`] rows are placeholders for entries the
/// barrier never classified. All three carry `NaN` ratios.
pub struct ActivityReport {
    /// Barrier parameter of the converged iterate.
    pub mu: Number,
    /// Status per user variable (codes above).
    ///
    /// [`AMBIGUOUS`] here includes genuine kinks whose coordinate is
    /// coupled to a neighbour: `q` is the Hessian diagonal, not the
    /// reduced curvature that generates the multiplier, so the ratio
    /// is `reduced/diagonal` (gh#763). Do not read the class as an
    /// answer to "is this bound at a kink" — [`reduced_activity`]
    /// answers that, one back-solve per coordinate.
    pub var_status: Vec<i8>,
    /// `Σ_i / q_i` per user variable; `NaN` where not classified.
    /// For an [`UNIDENTIFIED`] entry the value is `Σ/floor`, a lower
    /// bound on any honest ratio rather than the ratio itself, since
    /// `q` is below the identification floor there.
    ///
    /// `q` is `|H_ii|`, so at a kink this ratio is
    /// `reduced/diagonal` — `1` only where the coordinate is
    /// decoupled, and `μ`-independent, so a tighter solve does not
    /// move it. See [`Self::var_status`].
    pub var_ratio: Vec<Number>,
    /// Sign of the signed curvature `H_ii` (−1, 0, +1); the absolute
    /// value goes into `q`, so an indefinite direction is reported
    /// rather than hidden.
    pub var_q_sign: Vec<i8>,
    /// `s·z` differs from `μ` by more than a factor of ten on some
    /// side: off the central path, or the bound was relaxed.
    pub var_off_central_path: Vec<bool>,
    /// Classified inactive yet `r` non-negligible: barrier curvature
    /// where none should be.
    pub var_contaminated: Vec<bool>,
    /// The barrier diagonal `Σ_i = z/s` itself per user variable, both
    /// sides summed; 0 where not classified. In **natural (unscaled)
    /// units**, the repo's sensitivity-output contract: classification
    /// runs on the solver's scaled quantities (the ratio is
    /// scale-invariant), the report does not. The covariance roadmap's
    /// item 1 subtracts exactly this from the factor's natural-units
    /// reduced Hessian.
    pub var_sigma: Vec<Number>,
    /// Status per user constraint row.
    ///
    /// [`AMBIGUOUS`] here includes genuine kinks whose direction is
    /// coupled to the remaining free space: `q` is the curvature
    /// along the row's own gradient, not the reduced curvature that
    /// generates the multiplier, so the ratio is
    /// `reduced/directional` (gh#804). Do not read the class as an
    /// answer to "is this row at a kink" — [`reduced_row_activity`]
    /// answers that, one back-solve per row.
    pub row_status: Vec<i8>,
    /// `Σ_j / q_j` per user row; `NaN` where not classified.
    /// [`UNIDENTIFIED`] entries hold `Σ/floor` as for variables.
    ///
    /// `q` is the directional curvature `|∇dᵀH∇d|/‖∇d‖²`, so at a
    /// kink this ratio is `reduced/directional` — `1` only where the
    /// row's direction is decoupled, and `μ`-independent, so a
    /// tighter solve does not move it. See [`Self::row_status`].
    pub row_ratio: Vec<Number>,
    /// Sign of the signed row curvature `∇dⱼᵀ H ∇dⱼ`.
    pub row_q_sign: Vec<i8>,
    /// Central-path check per row, as for variables.
    pub row_off_central_path: Vec<bool>,
    /// Contamination check per row, as for variables.
    pub row_contaminated: Vec<bool>,
    /// The row barrier diagonal `Σ_j = v/s` per user row, both sides
    /// summed; 0 where not classified. In **natural (unscaled) units**
    /// like [`Self::var_sigma`], and RAW rather than the geometric
    /// weight the classification uses: item 1 restricts the normal to
    /// its own fitted block and applies its own `‖a‖²` factor there.
    pub row_sigma: Vec<Number>,
}

/// Scatter a compressed (bounded-entries-only) vector to full length
/// through its expansion matrix. Entries without that bound stay 0.
fn expand(compressed: &[Number], px: &Rc<dyn Matrix>, n: usize) -> Vec<Number> {
    let em = px
        .as_any()
        .downcast_ref::<ExpansionMatrix>()
        .expect("bound projection is an ExpansionMatrix (orig_ipopt_nlp builds no other kind)");
    let idx = em.expanded_pos_indices();
    assert_eq!(
        idx.len(),
        compressed.len(),
        "compressed bound vector length disagrees with its expansion",
    );
    let mut full = vec![0.0; n];
    for (k, &pos) in idx.iter().enumerate() {
        full[pos as usize] = compressed[k];
    }
    full
}

/// Presence mask for a bound side, from the same expansion.
fn present(px: &Rc<dyn Matrix>, n: usize) -> Vec<bool> {
    let em = px
        .as_any()
        .downcast_ref::<ExpansionMatrix>()
        .expect("bound projection is an ExpansionMatrix (orig_ipopt_nlp builds no other kind)");
    let mut mask = vec![false; n];
    for &pos in em.expanded_pos_indices() {
        mask[pos as usize] = true;
    }
    mask
}

/// The exact Hessian diagonal: one pass over the triplet structure for
/// the type `eval_h` builds today. The mat-vec fallback keeps any
/// future non-triplet `SymMatrix` correct, at O(n·nnz) cost.
fn hessian_diagonal(hess: &Rc<dyn pounce_linalg::SymMatrix>, n: usize) -> Vec<Number> {
    let mut diag = vec![0.0; n];
    if let Some(t) = hess.as_any().downcast_ref::<SymTMatrix>() {
        // triplet indices are 1-based (the GenTMatrix convention);
        // duplicates accumulate, matching mult_vector
        for ((&i, &j), &v) in t.irows().iter().zip(t.jcols()).zip(t.values()) {
            if i == j {
                diag[(i - 1) as usize] += v;
            }
        }
        return diag;
    }
    let space = DenseVectorSpace::new(n as i32);
    let mut e = DenseVector::new(space.clone());
    let mut he = DenseVector::new(space);
    for (i, d) in diag.iter_mut().enumerate() {
        e.values_mut().fill(0.0);
        e.values_mut()[i] = 1.0;
        he.values_mut().fill(0.0);
        hess.mult_vector(1.0, &e, 0.0, &mut he);
        // values_mut, not values: a zero product may have left the
        // output homogeneous (empty backing slice); this materializes
        *d = he.values_mut()[i];
    }
    diag
}

/// The per-variable pieces both classifiers measure against, with the
/// `user-scaling` change of variables (gh#486 stage 3) divided out and
/// the objective scale `df` still in: `Σ̃·d²` and `H̃_ii·d² = df·H_ii`.
/// The ratio of the two is invariant to both, so classification runs
/// here; the identification `floor` is not (it is one number shared
/// across entries), which is why the change of variables comes out
/// before it is formed.
///
/// Extracted so [`compute`] and [`reduced_activity`] cannot drift on
/// the frame conversion — the dimension leg 1 of
/// `sens_invariance_legs.rs` exists for.
struct VarFrame {
    /// Barrier diagonal `Σ_i` per var-x column.
    sigma: Vec<Number>,
    /// Exact-Lagrangian-Hessian diagonal per var-x column.
    diag: Vec<Number>,
    /// Identification floor shared by every entry, relative to the
    /// largest curvature anywhere on the diagonal rather than just the
    /// bounded entries, so a row-only model still measures `q` against
    /// the model's own scale.
    floor: Number,
}

fn var_frame(bs: &PdSensBacksolver, hess: &Rc<dyn pounce_linalg::SymMatrix>, n: usize) -> VarFrame {
    let d_var = bs.variable_scaling();
    let dv = |i: usize| -> Number { d_var.map_or(1.0, |d| d[i]) };
    // `Σ̃_i = df·Σ_i/d_i²`: the `d_i²` comes out here, the `df` at the
    // caller's export boundary (it cancels in every ratio, so
    // classification never sees it).
    let sigma = dense_to_vec(bs.barrier_sigma_x().as_ref())
        .iter()
        .enumerate()
        .map(|(i, &s)| s * dv(i) * dv(i))
        .collect();
    let diag: Vec<Number> = hessian_diagonal(hess, n)
        .iter()
        .enumerate()
        .map(|(i, &h)| h * dv(i) * dv(i))
        .collect();
    let max_abs_diag = diag.iter().fold(0.0, |a: Number, d| a.max(d.abs()));
    let floor = Number::EPSILON.sqrt() * max_abs_diag.max(1.0);
    VarFrame { sigma, diag, floor }
}

pub(crate) fn compute(bs: &PdSensBacksolver) -> ActivityReport {
    let (data, cq, nlp) = bs.activity_handles();

    // scoped borrows: the Cq getters below re-borrow the NLP (mutably,
    // for lazy evaluation) and the data, so nothing here may hold
    // either across a Cq call
    let mu = bs.barrier_mu();
    let (mult_z_l, mult_z_u, mult_v_l, mult_v_u, n, m_d) = {
        let d = data.borrow();
        let curr = d.curr.as_ref().expect("converged state has an iterate");
        (
            Rc::clone(&curr.z_l),
            Rc::clone(&curr.z_u),
            Rc::clone(&curr.v_l),
            Rc::clone(&curr.v_u),
            curr.x.dim() as usize,
            curr.s.dim() as usize,
        )
    };
    let (px_l, px_u, pd_l, pd_u, obj_scale, d_scale) = {
        let nl = nlp.borrow();
        (
            nl.px_l(),
            nl.px_u(),
            nl.pd_l(),
            nl.pd_u(),
            nl.obj_scaling_factor(),
            nl.d_scale_vec(),
        )
    };
    let cq = cq.borrow();

    // Per-variable factors of a `user-scaling` change of variables
    // (gh#486 stage 3), in var-x space; 1.0 everywhere when none ran.
    // Every internal x-space quantity below is a `d`-transform of the
    // model's own — writing `a_j` for the gradient of inequality row
    // `j`, since `d` is spoken for here: `ã = a ⊘ d`,
    // `H̃ = H ⊘ (d ⊗ d)`, `Σ̃ = Σ · df ⊘ (d ⊙ d)`. Undoing that here rather than only on
    // the exported `Σ` is what keeps a status from depending on the
    // conditioning the user asked for: the per-entry ratio `Σ/q` is
    // invariant, but the identification `floor` is a single number
    // shared across entries, so a non-uniform `d` moves entries across
    // it. `1.0` multiplies are exact, so an unscaled solve is
    // bit-identical to the pre-#486 path.
    let d_var = bs.variable_scaling();
    let dv = |i: usize| -> Number { d_var.map_or(1.0, |d| d[i]) };

    // --- variables, in internal space ------------------------------------
    let has_l = present(&px_l, n);
    let has_u = present(&px_u, n);
    let z_l = expand(&dense_to_vec(mult_z_l.as_ref()), &px_l, n);
    let z_u = expand(&dense_to_vec(mult_z_u.as_ref()), &px_u, n);
    // The solver's own slacks, deliberately, even when `Σ` below comes
    // from the declared frame (gh#654): these feed `off_path` only, and
    // "is `s·z` near `μ`" is a question about the central path, which is
    // the barrier's geometry and therefore the barrier's slacks. A
    // crossed-over point is off that path by construction and reads so
    // under either frame.
    let s_l = expand(&dense_to_vec(cq.curr_slack_x_l().as_ref()), &px_l, n);
    let s_u = expand(&dense_to_vec(cq.curr_slack_x_u().as_ref()), &px_u, n);
    // `Σ` and the Hessian diagonal, with the change of variables
    // already divided out and the objective scale still in; the `df`
    // comes off on export below (it cancels in every ratio, so
    // classification never sees it).
    let hess = cq.curr_exact_hessian();
    let VarFrame {
        sigma: sigma_x,
        diag,
        floor,
    } = var_frame(bs, &hess, n);

    let mut vars = vec![NOT_CLASSIFIED; n];
    for i in 0..n {
        if !(has_l[i] || has_u[i]) {
            continue;
        }
        let mut e = classify_entry(sigma_x[i], diag[i], floor, mu);
        e.off_path = (has_l[i] && off_path(s_l[i], z_l[i], mu))
            || (has_u[i] && off_path(s_u[i], z_u[i], mu));
        // the ratio is scale-invariant, so classification ran in the
        // solver's own space up to the change of variables already
        // divided out of `sigma_x` / `diag` above; the REPORTED sigma
        // follows the repo's natural-units contract, and what is left
        // to undo is the objective scale the internal z carries
        e.sigma /= obj_scale;
        vars[i] = e;
    }

    // --- inequality rows, in internal space -------------------------------
    let rhas_l = present(&pd_l, m_d);
    let rhas_u = present(&pd_u, m_d);
    let v_l = expand(&dense_to_vec(mult_v_l.as_ref()), &pd_l, m_d);
    let v_u = expand(&dense_to_vec(mult_v_u.as_ref()), &pd_u, m_d);
    let rs_l = expand(&dense_to_vec(cq.curr_slack_s_l().as_ref()), &pd_l, m_d);
    let rs_u = expand(&dense_to_vec(cq.curr_slack_s_u().as_ref()), &pd_u, m_d);
    let sigma_s = dense_to_vec(bs.barrier_sigma_s().as_ref());

    let jac_d = cq.curr_jac_d();
    // One pass over the Jacobian triplets gathers every row's support
    // and one pass over the Hessian triplets builds an adjacency view,
    // so each row's curvature costs its own support times its
    // neighbours instead of a full mat-vec pair per row (second
    // review). The mat-vec loop below remains the fallback for any
    // future non-triplet matrix types.
    let mut rows = vec![NOT_CLASSIFIED; m_d];
    let fast = match (
        jac_d.as_any().downcast_ref::<GenTMatrix>(),
        hess.as_any().downcast_ref::<SymTMatrix>(),
    ) {
        (Some(jt), Some(ht)) => {
            // gather and merge each row's entries (triplet duplicates
            // sum, matching mult_vector; indices are 1-based)
            let mut support: Vec<Vec<(usize, Number)>> = vec![Vec::new(); m_d];
            for ((&r, &c), &v) in jt.irows().iter().zip(jt.jcols()).zip(jt.values()) {
                let col = (c - 1) as usize;
                // `a = ã ⊙ d`: the row's own scale stays (the ratio
                // divides it out), the change of variables does not.
                support[(r - 1) as usize].push((col, v * dv(col)));
            }
            for sup in &mut support {
                sup.sort_unstable_by_key(|&(c, _)| c);
                sup.dedup_by(|a, b| {
                    if a.0 == b.0 {
                        b.1 += a.1;
                        true
                    } else {
                        false
                    }
                });
            }
            let mut adj: Vec<Vec<(usize, Number)>> = vec![Vec::new(); n];
            for ((&i, &l), &v) in ht.irows().iter().zip(ht.jcols()).zip(ht.values()) {
                let (a, b) = ((i - 1) as usize, (l - 1) as usize);
                // `H = H̃ ⊙ (d ⊗ d)`, matching the `d²` already taken
                // out of `diag` (which sets the shared floor).
                let v = v * dv(a) * dv(b);
                adj[a].push((b, v));
                if a != b {
                    adj[b].push((a, v));
                }
            }
            let mut scratch = vec![0.0; n];
            for j in 0..m_d {
                if !(rhas_l[j] || rhas_u[j]) {
                    continue;
                }
                let sup = &support[j];
                let norm2: Number = sup.iter().map(|&(_, g)| g * g).sum();
                rows[j] = if norm2 <= 0.0 {
                    zero_gradient_row(sigma_s[j], floor)
                } else {
                    for &(k, g) in sup {
                        scratch[k] = g;
                    }
                    let mut ghg = 0.0;
                    for &(k, gk) in sup {
                        let mut acc = 0.0;
                        for &(l, v) in &adj[k] {
                            acc += v * scratch[l];
                        }
                        ghg += gk * acc;
                    }
                    for &(k, _) in sup {
                        scratch[k] = 0.0;
                    }
                    // Σ·‖∇d‖² against curvature along the unit
                    // normal: invariant to rescaling the row; the
                    // report keeps the raw Σ
                    let mut e = classify_entry(sigma_s[j] * norm2, ghg / norm2, floor, mu);
                    e.sigma = sigma_s[j];
                    e
                };
            }
            true
        }
        _ => false,
    };
    if !fast {
        let mspace = DenseVectorSpace::new(m_d as i32);
        let mut e_row = DenseVector::new(mspace);
        let nspace = DenseVectorSpace::new(n as i32);
        let mut grad = DenseVector::new(nspace.clone());
        let mut hgrad = DenseVector::new(nspace);
        for j in 0..m_d {
            if !(rhas_l[j] || rhas_u[j]) {
                continue;
            }
            // ∇dⱼ = Jdᵀ eⱼ, then the curvature along the normal;
            // values_mut throughout because a zero product may leave
            // the output homogeneous (empty backing slice)
            e_row.values_mut().fill(0.0);
            e_row.values_mut()[j] = 1.0;
            grad.values_mut().fill(0.0);
            jac_d.trans_mult_vector(1.0, &e_row, 0.0, &mut grad);
            // `a = ã ⊙ d`, then `aᵀHa = uᵀH̃u` with `u = d ⊙ a`
            // (since `H = H̃ ⊙ (d ⊗ d)`) — so the vector handed to the
            // internal Hessian carries `d²`, and `norm2` carries `d`.
            let norm2: Number = grad
                .values_mut()
                .iter()
                .enumerate()
                .map(|(i, g)| (*g * dv(i)) * (*g * dv(i)))
                .sum();
            rows[j] = if norm2 <= 0.0 {
                zero_gradient_row(sigma_s[j], floor)
            } else {
                for (i, g) in grad.values_mut().iter_mut().enumerate() {
                    *g *= dv(i) * dv(i);
                }
                hgrad.values_mut().fill(0.0);
                hess.mult_vector(1.0, &grad, 0.0, &mut hgrad);
                let ghg: Number = {
                    let h = hgrad.values_mut();
                    grad.values_mut()
                        .iter()
                        .zip(h.iter())
                        .map(|(g, h)| g * h)
                        .sum()
                };
                let mut e = classify_entry(sigma_s[j] * norm2, ghg / norm2, floor, mu);
                e.sigma = sigma_s[j];
                e
            };
        }
    }
    for j in 0..m_d {
        if !(rhas_l[j] || rhas_u[j]) {
            continue;
        }
        rows[j].off_path = (rhas_l[j] && off_path(rs_l[j], v_l[j], mu))
            || (rhas_u[j] && off_path(rs_u[j], v_u[j], mu));
        // natural-units report, as for variables: the scaled row
        // multiplier carries df/dg and the scaled slack dg, so
        // Sigma_nat = Sigma * dg^2 / df
        let dg = d_scale.as_ref().map_or(1.0, |v| v[j]);
        rows[j].sigma *= dg * dg / obj_scale;
    }

    // --- scatter to user space --------------------------------------------
    // all Cq evaluation is done, so borrowing the NLP again is safe
    let nl = nlp.borrow();
    let n_full_x = nl.n_full_x() as usize;
    let n_full_g = nl.n_full_g() as usize;

    let fixed_entry = Entry {
        status: FIXED,
        ..NOT_CLASSIFIED
    };
    let mut var_full = vec![fixed_entry; n_full_x];
    for (i, e) in vars.iter().enumerate() {
        var_full[nl.var_x_to_full_x(i as Index) as usize] = *e;
    }

    let equality_entry = Entry {
        status: EQUALITY,
        ..NOT_CLASSIFIED
    };
    let mut row_full = vec![equality_entry; n_full_g];
    // BoundClassification's d_map is one ascending scan over the
    // user's g, so the j-th full-g index outside the c-block is
    // internal inequality row j
    let mut d_pos = 0usize;
    for (full_idx, slot) in row_full.iter_mut().enumerate() {
        if nl.full_g_to_c_block(full_idx as Index).is_none() {
            *slot = rows[d_pos];
            d_pos += 1;
        }
    }
    assert_eq!(d_pos, m_d, "inequality count disagrees with the c/d split");

    ActivityReport {
        mu,
        var_status: var_full.iter().map(|e| e.status).collect(),
        var_ratio: var_full.iter().map(|e| e.ratio).collect(),
        var_q_sign: var_full.iter().map(|e| e.q_sign).collect(),
        var_off_central_path: var_full.iter().map(|e| e.off_path).collect(),
        var_contaminated: var_full.iter().map(|e| e.contaminated).collect(),
        var_sigma: var_full.iter().map(|e| e.sigma).collect(),
        row_status: row_full.iter().map(|e| e.status).collect(),
        row_ratio: row_full.iter().map(|e| e.ratio).collect(),
        row_q_sign: row_full.iter().map(|e| e.q_sign).collect(),
        row_off_central_path: row_full.iter().map(|e| e.off_path).collect(),
        row_contaminated: row_full.iter().map(|e| e.contaminated).collect(),
        row_sigma: row_full.iter().map(|e| e.sigma).collect(),
    }
}

/// One coordinate's activity re-measured against the curvature that
/// actually generates its multiplier, rather than against the Hessian
/// diagonal (gh#763).
///
/// Parallel arrays, one entry per requested variable, in the order
/// they were requested. See [`reduced_activity`] for what the
/// quantities mean.
#[derive(Debug, Clone)]
pub struct ReducedActivityReport {
    /// Barrier parameter of the converged iterate, as
    /// [`ActivityReport::mu`].
    pub mu: Number,
    /// The user-space variable index each entry answers about.
    pub var: Vec<usize>,
    /// Status from the same rule [`ActivityReport::var_status`] uses,
    /// applied to [`Self::ratio`].
    pub status: Vec<i8>,
    /// `Σ_i / |q_i^red|`. `NaN` where nothing was classified.
    pub ratio: Vec<Number>,
    /// The reduced curvature `q_i^red` itself, signed, in **natural
    /// (unscaled) units** like [`ActivityReport::var_sigma`]. `NaN`
    /// for a [`FIXED`] variable and where the factor offers no
    /// reduced curvature (see [`reduced_activity`]).
    pub q_reduced: Vec<Number>,
    /// Sign of [`Self::q_reduced`] (−1, 0, +1); the absolute value is
    /// what the ratio divides by.
    pub q_sign: Vec<i8>,
    /// `Σ_i`, natural units, identical to the same variable's entry in
    /// [`ActivityReport::var_sigma`].
    pub sigma: Vec<Number>,
}

/// Why a [`reduced_activity`] call could not be answered.
pub(crate) enum ReducedActivityError {
    /// A requested index is not a user variable. Carries the offending
    /// index and the user TNLP's `n`.
    OutOfRange { got: usize, n_full_x: usize },
    /// The back-solve against the held factor failed.
    Backsolve,
}

/// [`compute`]'s per-variable classification, re-normalized by the
/// **reduced** curvature along each requested coordinate instead of
/// the Hessian diagonal — one back-solve against the held factor per
/// variable asked about (gh#763).
///
/// # What is different
///
/// [`compute`] forms a variable's ratio as `Σ_i / |H_ii|`. At a kink
/// the multiplier is not generated by the diagonal: it is generated by
/// the curvature **reduced** along that coordinate, i.e. what is left
/// after the other free variables re-optimize. Eliminating a free
/// partner `y` from `[[h, c], [c, m]]` leaves `h − c²/m`, and
/// `Σ = z/s` equals that, so the diagonal ratio is `1` only when the
/// coordinate is **decoupled** (`c = 0`). Couple it and a genuine kink
/// reads [`AMBIGUOUS`] at any tolerance, because the ratio is
/// `μ`-independent — routine on a collocation model, where coupling
/// between neighbouring coordinates is the normal case. This function
/// is the answer to "is that AMBIGUOUS entry a kink?": on the reduced
/// normalizer the same kink reads [`WEAKLY_ACTIVE`] whatever it is
/// coupled to.
///
/// # How it is computed
///
/// The reduced curvature the barrier subproblem sees along coordinate
/// `i` is the reciprocal of the `i`-th diagonal entry of the inverse
/// of the barrier-augmented KKT matrix — one column of `K⁻¹`, so one
/// back-solve — with the barrier's own contribution at `i` taken back
/// off:
///
/// ```text
/// q_i^red = 1 / (K⁻¹)_ii − Σ_i
/// ```
///
/// Adding `diag(Σ)` to the free block shifts entry `i` by exactly
/// `Σ_i`, which the subtraction removes; it also shifts every OTHER
/// free coordinate, which is deliberately kept. That is what makes the
/// answer the curvature the multiplier is generated against: a
/// strongly active neighbour carries `Σ = O(1/μ)` and so does not
/// re-optimize, an inactive one carries `O(μ)` and re-optimizes
/// freely, and the elimination weights each accordingly. Eliminating
/// against the bare Hessian would let a pinned neighbour move.
///
/// `K` here is the natural-units KKT matrix, so `q_i^red` and `Σ_i`
/// are both in the model's own units and the ratio is invariant to
/// `nlp_scaling_method` and to a `user-scaling` change of variables,
/// exactly as [`compute`]'s is.
///
/// # Cost, and why this is not the default
///
/// The correct normalizer is the reciprocal diagonal of an *inverse*,
/// and there is no diagonal-of-the-inverse shortcut: it is one
/// back-solve per coordinate. Classifying every bounded variable this
/// way is `n` back-solves, which on a 62k-variable model is not a
/// post-solve diagnostic any more. So [`compute`] keeps the `O(nnz)`
/// diagonal and this is the on-demand refinement — the natural call is
/// over the [`AMBIGUOUS`] entries of a report, which is a handful on
/// the models where the question arises.
///
/// Pass those indices in ONE call: the back-solves batch, and the
/// per-call fixed cost — one pass over the Hessian for the shared
/// identification floor — is paid once rather than per index.
///
/// # Edge cases
///
/// * A [`FIXED`] variable has no column in the factor: status
///   [`FIXED`], `NaN` curvature and ratio, as in the report.
/// * A variable with no finite bound gets its `q_i^red` (it is
///   well-defined and already paid for) but status [`UNBOUNDED`] and a
///   `NaN` ratio: there is no bound question to answer.
/// * A coordinate the constraints determine outright has
///   `(K⁻¹)_ii = 0`: no direction left to reduce along, so `q_i^red`
///   is infinite and the ratio is `0`, i.e. [`INACTIVE`] — whatever
///   holds that coordinate, it is not its bound. A `(K⁻¹)_ii` at
///   roundoff level reads the same from either side of zero.
/// * A `NaN` from the back-solve reports [`UNIDENTIFIED`] with a `NaN`
///   curvature and ratio.
/// * `|q_i^red|` below the same identification floor [`compute`] uses
///   reports [`UNIDENTIFIED`], as there.
/// * Under `obj_scaling_factor < 0` — the documented way to maximize —
///   `Sigma` and `q_i^red` are both reported with the sign the
///   natural-units contract gives them (negative, as
///   [`ActivityReport::var_sigma`] is), but the CLASSIFICATION runs on
///   the objective-scale-positive orientation [`compute`] classifies
///   in, so a status here means the same thing at either sign of `df`.
///
/// Variable bounds only. The row path carries the same un-reduced
/// distinction — see [`reduced_row_activity`], which is this function
/// one KKT block over (gh#804).
pub(crate) fn reduced_activity(
    bs: &PdSensBacksolver,
    user_vars: &[usize],
) -> Result<ReducedActivityReport, ReducedActivityError> {
    let (data, cq, nlp) = bs.activity_handles();
    let mu = bs.barrier_mu();
    let n = {
        let d = data.borrow();
        d.curr
            .as_ref()
            .expect("converged state has an iterate")
            .x
            .dim() as usize
    };

    // full-x in, var-x rows out, through the same map the report
    // scatters through: reading the user index as a factor row returns
    // a NEIGHBORING variable's answer wherever a fixed variable
    // precedes it (the gh#450 hazard).
    let n_full_x = bs.n_full_x() as usize;
    let mut rows: Vec<Option<usize>> = Vec::with_capacity(user_vars.len());
    for &i in user_vars {
        if i >= n_full_x {
            return Err(ReducedActivityError::OutOfRange { got: i, n_full_x });
        }
        rows.push(bs.full_x_to_var_x(i as Index).map(|r| r as usize));
    }

    let (px_l, px_u, obj_scale) = {
        let nl = nlp.borrow();
        (nl.px_l(), nl.px_u(), nl.obj_scaling_factor())
    };
    let has_l = present(&px_l, n);
    let has_u = present(&px_u, n);
    // the NLP borrow above is dropped: the Cq getter re-borrows it
    // mutably for lazy evaluation
    let frame = {
        let cq = cq.borrow();
        let hess = cq.curr_exact_hessian();
        var_frame(bs, &hess, n)
    };
    // `Σ` and the floor come out of `var_frame` with the change of
    // variables divided out but the objective scale still in; the
    // reciprocal `(K⁻¹)_ii` below is natural units already, so both
    // sides of the subtraction meet there.
    let floor = frame.floor / obj_scale.abs();

    // One unit RHS per requested variable, batched against the held
    // factor. Chunked so the buffers stay bounded when a caller hands
    // over a long list rather than the report's ambiguous entries.
    const CHUNK: usize = 64;
    let dim = bs.dim();
    let solve_rows: Vec<usize> = rows.iter().flatten().copied().collect();
    let mut kinv: Vec<Number> = Vec::with_capacity(solve_rows.len());
    for chunk in solve_rows.chunks(CHUNK) {
        let k = chunk.len();
        let mut rhs = vec![0.0; k * dim];
        let mut lhs = vec![0.0; k * dim];
        for (c, &r) in chunk.iter().enumerate() {
            rhs[c * dim + r] = 1.0;
        }
        if !bs.solve_many(&rhs, &mut lhs, k) {
            return Err(ReducedActivityError::Backsolve);
        }
        for (c, &r) in chunk.iter().enumerate() {
            kinv.push(lhs[c * dim + r]);
        }
    }

    let mut out = ReducedActivityReport {
        mu,
        var: user_vars.to_vec(),
        status: Vec::with_capacity(user_vars.len()),
        ratio: Vec::with_capacity(user_vars.len()),
        q_reduced: Vec::with_capacity(user_vars.len()),
        q_sign: Vec::with_capacity(user_vars.len()),
        sigma: Vec::with_capacity(user_vars.len()),
    };
    let mut next = 0usize;
    for &row in &rows {
        let Some(row) = row else {
            out.status.push(FIXED);
            out.ratio.push(Number::NAN);
            out.q_reduced.push(Number::NAN);
            out.q_sign.push(0);
            out.sigma.push(0.0);
            continue;
        };
        let d = kinv[next];
        next += 1;
        let sigma = frame.sigma[row] / obj_scale;
        if !d.is_finite() {
            out.status.push(UNIDENTIFIED);
            out.ratio.push(Number::NAN);
            out.q_reduced.push(Number::NAN);
            out.q_sign.push(0);
            out.sigma.push(sigma);
            continue;
        }
        // `(K⁻¹)_ii = 0` — the constraints determine the coordinate
        // outright, so there is no direction left to reduce along —
        // sends `q` to an infinity the ratio divides to zero, i.e.
        // INACTIVE: whatever holds the coordinate there, it is not the
        // bound. A `(K⁻¹)_ii` at roundoff level lands in the same
        // class from either side of zero, which is why no guard
        // branches on its sign; a genuinely negative reduced curvature
        // is modest in magnitude and reports through `q_sign`, exactly
        // as an indefinite `H_ii` does in the report.
        let q = 1.0 / d - sigma;
        // classify in the same orientation `compute` does: it runs the
        // rule on the df-in `Sigma` (internal `z/s`, non-negative),
        // dividing the objective scale out only on export. Here both
        // sides are already natural, so a NEGATIVE df -- the
        // documented way to maximize -- would otherwise hand the rule
        // a negative ratio and read a pinned bound as INACTIVE.
        let sgn = if obj_scale < 0.0 { -1.0 } else { 1.0 };
        let e = classify_entry(sigma * sgn, q * sgn, floor, mu);
        // no finite bound: the curvature is still an answer, the
        // activity question is not
        let bounded = has_l[row] || has_u[row];
        out.status.push(if bounded { e.status } else { UNBOUNDED });
        out.ratio.push(if bounded { e.ratio } else { Number::NAN });
        out.q_reduced.push(q);
        out.q_sign.push(sign_of(q));
        out.sigma.push(sigma);
    }
    debug_assert_eq!(next, kinv.len(), "every free row consumed its solve");
    Ok(out)
}

/// One constraint row's activity re-measured against the curvature
/// that actually generates its multiplier, rather than against the
/// directional curvature along its own gradient (gh#804).
///
/// Parallel arrays, one entry per requested row, in the order they
/// were requested. See [`reduced_row_activity`] for what the
/// quantities mean.
#[derive(Debug, Clone)]
pub struct ReducedRowActivityReport {
    /// Barrier parameter of the converged iterate, as
    /// [`ActivityReport::mu`].
    pub mu: Number,
    /// The user-space constraint index each entry answers about.
    pub row: Vec<usize>,
    /// Status from the same rule [`ActivityReport::row_status`] uses,
    /// applied to [`Self::ratio`].
    pub status: Vec<i8>,
    /// `Σ_j‖∇dⱼ‖² / |q_j^red|`. `NaN` where nothing was classified.
    ///
    /// The numerator is the geometric barrier weight, not the raw
    /// `Σ_j` [`Self::sigma`] reports — the same pairing
    /// [`ActivityReport::row_ratio`] uses, so the two ratios are
    /// directly comparable and agree on a decoupled row.
    pub ratio: Vec<Number>,
    /// The reduced curvature `q_j^red` itself, signed, along the
    /// **unit** normal `∇dⱼ/‖∇dⱼ‖` and in **natural (unscaled)
    /// units** — the same quantity and units
    /// [`ActivityReport::row_q_sign`] takes the sign of, so it is what
    /// `|∇dⱼᵀH∇dⱼ|/‖∇dⱼ‖²` would have been had the curvature been
    /// reduced. `NaN` for an [`EQUALITY`] row and where the back-solve
    /// offers no reduced curvature (see [`reduced_row_activity`]).
    pub q_reduced: Vec<Number>,
    /// Sign of [`Self::q_reduced`] (−1, 0, +1); the absolute value is
    /// what the ratio divides by.
    pub q_sign: Vec<i8>,
    /// `Σ_j`, RAW (not the geometric weight the ratio uses) and in
    /// natural units, identical to the same row's entry in
    /// [`ActivityReport::row_sigma`].
    pub sigma: Vec<Number>,
}

/// Why a [`reduced_row_activity`] call could not be answered.
pub(crate) enum ReducedRowActivityError {
    /// A requested index is not a user constraint. Carries the
    /// offending index and the user TNLP's `m`.
    OutOfRange { got: usize, n_full_g: usize },
    /// The back-solve against the held factor failed.
    Backsolve,
}

/// `‖∇dⱼ‖²` for the requested internal inequality rows, in the frame
/// [`compute`] classifies in: `a = ã ⊙ d`, so the change of variables
/// is divided out and the row's own `d_scale` is still in. The caller
/// divides that `dg²` out to reach natural units.
///
/// One pass over the Jacobian triplets for the whole batch, with the
/// mat-vec loop kept as the fallback for any future non-triplet
/// matrix type — the same two paths, and the same duplicate-summing
/// convention, [`compute`] uses.
fn row_norm2(
    jac_d: &Rc<dyn Matrix>,
    d_var: Option<&[Number]>,
    wanted: &[Option<usize>],
    n: usize,
    m_d: usize,
    out_len: usize,
) -> Vec<Number> {
    let dv = |i: usize| -> Number { d_var.map_or(1.0, |d| d[i]) };
    let mut norm2 = vec![0.0; out_len];
    if let Some(jt) = jac_d.as_any().downcast_ref::<GenTMatrix>() {
        let mut support: Vec<Vec<(usize, Number)>> = vec![Vec::new(); out_len];
        for ((&r, &c), &v) in jt.irows().iter().zip(jt.jcols()).zip(jt.values()) {
            let Some(slot) = wanted[(r - 1) as usize] else {
                continue;
            };
            let col = (c - 1) as usize;
            support[slot].push((col, v * dv(col)));
        }
        for (slot, sup) in support.iter_mut().enumerate() {
            // triplet duplicates sum before the square, matching
            // `mult_vector` and `compute`'s own gather
            sup.sort_unstable_by_key(|&(c, _)| c);
            sup.dedup_by(|a, b| {
                if a.0 == b.0 {
                    b.1 += a.1;
                    true
                } else {
                    false
                }
            });
            norm2[slot] = sup.iter().map(|&(_, g)| g * g).sum();
        }
        return norm2;
    }
    let mspace = DenseVectorSpace::new(m_d as i32);
    let mut e_row = DenseVector::new(mspace);
    let nspace = DenseVectorSpace::new(n as i32);
    let mut grad = DenseVector::new(nspace);
    for (j, slot) in wanted.iter().enumerate() {
        let Some(slot) = *slot else { continue };
        // values_mut throughout: a zero product may leave the output
        // homogeneous (empty backing slice)
        e_row.values_mut().fill(0.0);
        e_row.values_mut()[j] = 1.0;
        grad.values_mut().fill(0.0);
        jac_d.trans_mult_vector(1.0, &e_row, 0.0, &mut grad);
        norm2[slot] = grad
            .values_mut()
            .iter()
            .enumerate()
            .map(|(i, g)| (*g * dv(i)) * (*g * dv(i)))
            .sum();
    }
    norm2
}

/// [`compute`]'s per-row classification, re-normalized by the
/// **reduced** curvature along each requested row's gradient instead
/// of the directional curvature `∇dᵀH∇d/‖∇d‖²` — one back-solve
/// against the held factor per row asked about (gh#804).
///
/// The row counterpart of [`reduced_activity`], and the same defect:
/// gh#763 for rows.
///
/// # What is different
///
/// [`compute`] forms a row's ratio as `Σ_j‖∇d‖⁴/|∇dᵀH∇d|`. That
/// denominator is a genuine curvature along the row's own gradient
/// direction — strictly better than the variable path's bare `H_ii`,
/// which is why gh#763 fixed the variables first — but it is not a
/// *reduced* curvature: it does not account for the other free
/// coordinates re-optimizing. The quantity that generates a row's
/// multiplier is what is left after that elimination, exactly as for
/// a variable, so a row's `r` is `reduced/directional` and equals `1`
/// only where the row's direction is decoupled from the remaining
/// free space. Couple it and a genuine row kink falls out of the
/// `[1e-1, 1e1]` band and reads [`AMBIGUOUS`] — at any tolerance,
/// because that ratio is `μ`-independent, so re-solving tighter
/// reports the same thing.
///
/// # How it is computed
///
/// The row's own value IS a coordinate of the KKT system: the slack
/// `s_j` the barrier acts on, tied to the model by `dⱼ(x) = s_j`. So
/// the reduced curvature the barrier subproblem sees along the row is
/// the reciprocal of that coordinate's diagonal entry of the inverse
/// — one back-solve against a unit right-hand side in the `s` block —
/// with the row's own barrier contribution taken back off:
///
/// ```text
/// q_j^raw = 1 / (K⁻¹)_{sⱼsⱼ} − Σ_j        (per unit of dⱼ)
/// q_j^red = q_j^raw · ‖∇dⱼ‖²              (per unit of x, reported)
/// ```
///
/// Driving the `x` block with the row's gradient `∇dⱼ` instead gives
/// the identical number — `∇dⱼᵀK⁻¹∇dⱼ = (K⁻¹)_{sⱼsⱼ}` for every `j`,
/// since `∇dⱼ` reaches the system only through the row it defines —
/// so the slack unit vector is used: it needs no gradient assembled
/// into the right-hand side, and it is the same call
/// [`reduced_activity`] makes one block over.
///
/// Everything the elimination weights is what makes this the
/// curvature the multiplier is generated against, exactly as in
/// [`reduced_activity`]: a strongly active neighbour carries
/// `Σ = O(1/μ)` and does not re-optimize, an inactive one carries
/// `O(μ)` and re-optimizes freely.
///
/// The `‖∇dⱼ‖²` puts the answer along the **unit** normal, which is
/// where [`compute`]'s `q` lives and what the shared identification
/// floor is scaled for. Paired with the geometric weight
/// `Σ_j‖∇dⱼ‖²` in the numerator — [`compute`]'s pairing — the ratio
/// is invariant to rescaling the row, so `d → c·d` does not change a
/// status here any more than it does there.
///
/// `K` is the natural-units KKT matrix, so the ratio is also
/// invariant to `nlp_scaling_method` and to a `user-scaling` change
/// of variables.
///
/// # Cost, and why this is not the default
///
/// One back-solve per row, for the same reason [`reduced_activity`]
/// costs one per variable: the correct normalizer is a diagonal entry
/// of an *inverse*, and there is no shortcut to it. So [`compute`]
/// keeps the `O(nnz)` directional curvature and this is the on-demand
/// refinement — the natural call is over the [`AMBIGUOUS`] rows of a
/// report. Pass them in ONE call: the back-solves batch, and the
/// per-call fixed cost (one pass over the Hessian for the shared
/// identification floor, one over the Jacobian for the norms) is paid
/// once rather than per row.
///
/// # Edge cases
///
/// * An equality row has no slack and no barrier multiplier pair:
///   status [`EQUALITY`], `NaN` curvature and ratio, as in the report.
/// * An inequality row with no finite bound gets its `q_j^red` but
///   status [`UNBOUNDED`] and a `NaN` ratio: no bound question to
///   answer.
/// * A row whose gradient vanishes at the iterate has no direction to
///   measure curvature along: [`UNIDENTIFIED`] with the raw
///   `Σ/floor` lower bound, exactly as [`compute`] reports it.
/// * A row the constraints determine outright has `(K⁻¹)_{sⱼsⱼ} = 0`:
///   no direction left to reduce along, so `q_j^red` is infinite and
///   the ratio is `0`, i.e. [`INACTIVE`] — whatever holds that row, it
///   is not its own bound.
/// * A `NaN` from the back-solve reports [`UNIDENTIFIED`] with a `NaN`
///   curvature and ratio.
/// * `|q_j^red|` below the same identification floor [`compute`] uses
///   reports [`UNIDENTIFIED`], as there.
/// * Under `obj_scaling_factor < 0` — the documented way to maximize —
///   `Σ` and `q_j^red` are both reported with the sign the
///   natural-units contract gives them, but the CLASSIFICATION runs on
///   the objective-scale-positive orientation [`compute`] classifies
///   in, so a status here means the same thing at either sign of `df`.
pub(crate) fn reduced_row_activity(
    bs: &PdSensBacksolver,
    user_rows: &[usize],
) -> Result<ReducedRowActivityReport, ReducedRowActivityError> {
    let (data, cq, nlp) = bs.activity_handles();
    let mu = bs.barrier_mu();
    let (n, m_d) = {
        let d = data.borrow();
        let curr = d.curr.as_ref().expect("converged state has an iterate");
        (curr.x.dim() as usize, curr.s.dim() as usize)
    };

    // full-g in, d-block rows out, through the NLP's own c/d map:
    // reading the user index as an inequality position returns a
    // NEIGHBORING row's answer wherever an equality precedes it (the
    // gh#450 hazard, one block over). It is the same map
    // `Solver::d_multiplier_rows` addresses the `y_d` block with
    // (gh#910), deliberately rather than a second ascending scan that
    // agrees with it today: the gate that accepts a strictly active
    // inequality's `dλ/dp` classifies with THIS function and then
    // reads THAT row, so a disagreement between the two would classify
    // one row and answer about another.
    let n_full_g = bs.n_full_g() as usize;
    let d_index: Vec<Option<usize>> = {
        let nl = nlp.borrow();
        let map: Vec<Option<usize>> = (0..n_full_g)
            .map(|g| nl.full_g_to_d_block(g as Index).map(|p| p as usize))
            .collect();
        // the same invariant [`compute`] asserts after its own scan,
        // restated here because every `sigma_s` / `d_scale` index
        // below rests on it
        assert_eq!(
            map.iter().filter(|p| p.is_some()).count(),
            m_d,
            "inequality count disagrees with the c/d split"
        );
        map
    };
    let mut rows: Vec<Option<usize>> = Vec::with_capacity(user_rows.len());
    for &j in user_rows {
        if j >= n_full_g {
            return Err(ReducedRowActivityError::OutOfRange { got: j, n_full_g });
        }
        rows.push(d_index[j]);
    }

    let (pd_l, pd_u, obj_scale, d_scale) = {
        let nl = nlp.borrow();
        (
            nl.pd_l(),
            nl.pd_u(),
            nl.obj_scaling_factor(),
            nl.d_scale_vec(),
        )
    };
    let rhas_l = present(&pd_l, m_d);
    let rhas_u = present(&pd_u, m_d);
    let sigma_s = dense_to_vec(bs.barrier_sigma_s().as_ref());

    // Which internal rows the batch asks about, and where each one's
    // norm lands. A row requested twice shares one slot and one gather.
    let mut slot_of: Vec<Option<usize>> = vec![None; m_d];
    let mut n_slots = 0usize;
    for row in rows.iter().flatten() {
        if slot_of[*row].is_none() {
            slot_of[*row] = Some(n_slots);
            n_slots += 1;
        }
    }

    // the NLP borrow above is dropped: the Cq getters re-borrow it
    // mutably for lazy evaluation
    let (floor, norm2) = {
        let cq = cq.borrow();
        let hess = cq.curr_exact_hessian();
        // only the shared identification floor is wanted from the
        // frame; it is the one number a per-entry ratio cannot supply
        let floor = var_frame(bs, &hess, n).floor;
        let jac_d = cq.curr_jac_d();
        let norm2 = row_norm2(&jac_d, bs.variable_scaling(), &slot_of, n, m_d, n_slots);
        (floor, norm2)
    };
    // the floor comes out of `var_frame` with the objective scale
    // still in; the reciprocal `(K⁻¹)_{ss}` below is natural units
    // already, so both sides of the comparison meet there.
    let floor = floor / obj_scale.abs();

    // One unit RHS per requested row, in the `s` block, batched
    // against the held factor. Chunked so the buffers stay bounded
    // when a caller hands over a long list rather than the report's
    // ambiguous rows.
    const CHUNK: usize = 64;
    let dim = bs.dim();
    let s_offset = bs.block_dims()[0];
    let solve_rows: Vec<usize> = rows.iter().flatten().map(|&r| s_offset + r).collect();
    let mut kinv: Vec<Number> = Vec::with_capacity(solve_rows.len());
    for chunk in solve_rows.chunks(CHUNK) {
        let k = chunk.len();
        let mut rhs = vec![0.0; k * dim];
        let mut lhs = vec![0.0; k * dim];
        for (c, &r) in chunk.iter().enumerate() {
            rhs[c * dim + r] = 1.0;
        }
        if !bs.solve_many(&rhs, &mut lhs, k) {
            return Err(ReducedRowActivityError::Backsolve);
        }
        for (c, &r) in chunk.iter().enumerate() {
            kinv.push(lhs[c * dim + r]);
        }
    }

    let mut out = ReducedRowActivityReport {
        mu,
        row: user_rows.to_vec(),
        status: Vec::with_capacity(user_rows.len()),
        ratio: Vec::with_capacity(user_rows.len()),
        q_reduced: Vec::with_capacity(user_rows.len()),
        q_sign: Vec::with_capacity(user_rows.len()),
        sigma: Vec::with_capacity(user_rows.len()),
    };
    // `compute` runs the rule on the df-in `Sigma` (internal `v/s`,
    // non-negative), dividing the objective scale out only on export.
    // Here every quantity is already natural, so a NEGATIVE df -- the
    // documented way to maximize -- would otherwise hand the rule a
    // negative ratio and read a pinned row as INACTIVE.
    let sgn = if obj_scale < 0.0 { -1.0 } else { 1.0 };
    let mut next = 0usize;
    for &row in &rows {
        let Some(row) = row else {
            out.status.push(EQUALITY);
            out.ratio.push(Number::NAN);
            out.q_reduced.push(Number::NAN);
            out.q_sign.push(0);
            out.sigma.push(0.0);
            continue;
        };
        let d = kinv[next];
        next += 1;
        // natural units, as `compute` exports them: the scaled row
        // multiplier carries df/dg and the scaled slack dg
        let dg = d_scale.as_ref().map_or(1.0, |v| v[row]);
        let sigma = sigma_s[row] * dg * dg / obj_scale;
        // `a = ã ⊙ d` leaves the row's own `dg` in, and the reported
        // curvature is natural, so it comes back out here
        let norm2 = norm2[slot_of[row].expect("every solved row has a norm slot")] / (dg * dg);
        if norm2 <= 0.0 {
            // no direction to measure curvature along, exactly as in
            // the report -- and the geometric weight is degenerate at
            // zero gradient, so the raw `Σ/floor` lower bound stands
            let e = zero_gradient_row(sigma * sgn, floor);
            out.status.push(e.status);
            out.ratio.push(e.ratio);
            out.q_reduced.push(Number::NAN);
            out.q_sign.push(0);
            out.sigma.push(sigma);
            continue;
        }
        if !d.is_finite() {
            out.status.push(UNIDENTIFIED);
            out.ratio.push(Number::NAN);
            out.q_reduced.push(Number::NAN);
            out.q_sign.push(0);
            out.sigma.push(sigma);
            continue;
        }
        // `(K⁻¹)_{ss} = 0` -- the rest of the model determines the
        // row's value outright, so there is no direction left to
        // reduce along -- sends `q` to an infinity the ratio divides
        // to zero, i.e. INACTIVE: whatever holds the row there, it is
        // not its own bound. A `(K⁻¹)_{ss}` at roundoff level lands in
        // the same class from either side of zero, which is why no
        // guard branches on its sign; a genuinely negative reduced
        // curvature is modest in magnitude and reports through
        // `q_sign`, exactly as an indefinite `∇dᵀH∇d` does in the
        // report.
        let q = (1.0 / d - sigma) * norm2;
        // the geometric weight against the curvature along the unit
        // normal: `compute`'s pairing, so a decoupled row's ratio here
        // IS the report's
        let e = classify_entry(sigma * norm2 * sgn, q * sgn, floor, mu);
        let bounded = rhas_l[row] || rhas_u[row];
        out.status.push(if bounded { e.status } else { UNBOUNDED });
        out.ratio.push(if bounded { e.ratio } else { Number::NAN });
        out.q_reduced.push(q);
        out.q_sign.push(sign_of(q));
        out.sigma.push(sigma);
    }
    debug_assert_eq!(next, kinv.len(), "every inequality row consumed its solve");
    Ok(out)
}

/// The gradient of one user constraint row at the converged iterate,
/// in user variable order (length `n_full_x`) and **natural (unscaled)
/// units**: the internal Jacobian row carries the solver's per-row
/// scale, which is divided out here per the sensitivity-output
/// contract. Works for equality and inequality rows alike; entries for
/// `make_parameter`-removed fixed variables are 0 because the solve
/// dropped their columns.
pub(crate) fn row_normal(bs: &PdSensBacksolver, user_row: usize) -> Result<Vec<Number>, usize> {
    let (data, cq, nlp) = bs.activity_handles();
    let n = {
        let d = data.borrow();
        d.curr
            .as_ref()
            .expect("converged state has an iterate")
            .x
            .dim() as usize
    };
    // position of the row within its own c/d block, by the same
    // ascending scan the report's scatter uses
    let c_pos = {
        let nl = nlp.borrow();
        if user_row >= nl.n_full_g() as usize {
            return Err(nl.n_full_g() as usize);
        }
        nl.full_g_to_c_block(user_row as Index)
    };
    let block_pos = match c_pos {
        Some(p) => p as usize,
        None => {
            let nl = nlp.borrow();
            (0..user_row)
                .filter(|&g| nl.full_g_to_c_block(g as Index).is_none())
                .count()
        }
    };

    let row_scale = {
        let nl = nlp.borrow();
        let sv = if c_pos.is_some() {
            nl.c_scale_vec()
        } else {
            nl.d_scale_vec()
        };
        sv.map_or(1.0, |v| v[block_pos])
    };
    let cq = cq.borrow();
    let jac = if c_pos.is_some() {
        cq.curr_jac_c()
    } else {
        cq.curr_jac_d()
    };
    let m_block = jac.n_rows() as usize;
    let mspace = DenseVectorSpace::new(m_block as i32);
    let mut e_row = DenseVector::new(mspace);
    let nspace = DenseVectorSpace::new(n as i32);
    let mut grad = DenseVector::new(nspace);
    e_row.values_mut().fill(0.0);
    e_row.values_mut()[block_pos] = 1.0;
    grad.values_mut().fill(0.0);
    jac.trans_mult_vector(1.0, &e_row, 0.0, &mut grad);

    let d_var = bs.variable_scaling();
    let nl = nlp.borrow();
    let n_full_x = nl.n_full_x() as usize;
    let mut full = vec![0.0; n_full_x];
    let g = grad.values_mut();
    for (i, slot) in g.iter().enumerate() {
        // `∇g̃ = (∇g ⊘ d) · row_scale`, so both come back out here
        // (gh#486 stage 3).
        let dx = d_var.map_or(1.0, |d| d[i]);
        full[nl.var_x_to_full_x(i as Index) as usize] = *slot * dx / row_scale;
    }
    Ok(full)
}

/// The exact Lagrangian Hessian times a user-space vector, in user
/// variable order and **natural (unscaled) units**: the internal
/// Hessian carries the objective scale, divided out here per the
/// sensitivity-output contract. Entries for `make_parameter`-removed
/// fixed variables are 0 in and out (their columns left the solve).
/// Serves the covariance roadmap's item 2: the tangent-recovered
/// reduced Hessian is `T^T (H T)`, one product per fitted column.
pub(crate) fn hessian_vec(bs: &PdSensBacksolver, v_full: &[Number]) -> Result<Vec<Number>, usize> {
    let (data, cq, nlp) = bs.activity_handles();
    let n = {
        let d = data.borrow();
        d.curr
            .as_ref()
            .expect("converged state has an iterate")
            .x
            .dim() as usize
    };
    let (n_full_x, obj_scale) = {
        let nl = nlp.borrow();
        (nl.n_full_x() as usize, nl.obj_scaling_factor())
    };
    if v_full.len() != n_full_x {
        return Err(n_full_x);
    }

    // `H = H̃ ⊙ (d ⊗ d)` under a change of variables (gh#486 stage 3),
    // so `H v = d ⊙ (H̃ (d ⊙ v))`: the factor goes in with the vector
    // and comes back out of the product.
    let d_var = bs.variable_scaling();
    let nspace = DenseVectorSpace::new(n as i32);
    let mut v_int = DenseVector::new(nspace.clone());
    let mut hv = DenseVector::new(nspace);
    {
        let nl = nlp.borrow();
        let vals = v_int.values_mut();
        vals.fill(0.0);
        for i in 0..n {
            let dx = d_var.map_or(1.0, |d| d[i]);
            vals[i] = v_full[nl.var_x_to_full_x(i as Index) as usize] * dx;
        }
    }
    let hess = {
        let cq = cq.borrow();
        cq.curr_exact_hessian()
    };
    hess.mult_vector(1.0, &v_int, 0.0, &mut hv);

    let nl = nlp.borrow();
    let mut out = vec![0.0; n_full_x];
    let h = hv.values_mut();
    for (i, slot) in h.iter().enumerate() {
        let dx = d_var.map_or(1.0, |d| d[i]);
        out[nl.var_x_to_full_x(i as Index) as usize] = *slot * dx / obj_scale;
    }
    Ok(out)
}