pounce_sens_core/

activity.rs

//! 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 comes from the engine through
//! [`BarrierIterate::activity_data`]: the bound multipliers on the
//! iterate, the engine's own slacks, the barrier diagonal `Σ`, the
//! barrier parameter, and the exact Lagrangian Hessian, so `H` is never
//! recovered from the barrier-augmented factor. pounce's interior-point
//! solver reports the declared-frame diagonal there when the held
//! iterate came from crossover (gh#654). The reduced classifications
//! add one back-solve per entry through [`crate::backsolver::SensBacksolver::solve_many`].
//!
//! 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::Number;
use pounce_linalg::Matrix;
use pounce_linalg::dense_vector::{DenseVector, DenseVectorSpace};
use pounce_linalg::triplet::{GenTMatrix, SymTMatrix};

use crate::barrier::{ActivityData, BarrierIterate, EngineScaling, ModelMaps};

// The status codes and the classification rule live in `activity_kernel`,
// so the convex arm decides what a kink is with the same code rather than a
// parallel one. Re-exported here so a caller of this module has them at hand.
pub use crate::activity_kernel::{
    AMBIGUOUS, EQUALITY, FIXED, INACTIVE, STRONGLY_ACTIVE, UNBOUNDED, UNIDENTIFIED, WEAKLY_ACTIVE,
};
use crate::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 the positions of its bounds. Entries without that bound stay 0.
fn expand(compressed: &[Number], pos: &[usize], n: usize) -> Vec<Number> {
    let mut full = vec![0.0; n];
    for (k, &p) in pos.iter().enumerate() {
        full[p] = compressed[k];
    }
    full
}

/// Presence mask for a bound side, from the same positions.
fn present(pos: &[usize], n: usize) -> Vec<bool> {
    let mut mask = vec![false; n];
    for &p in pos {
        mask[p] = true;
    }
    mask
}

/// What every classification reads from the engine: the converged
/// quantities, the model maps and the scaling, each checked against the
/// engine's block lengths before anything indexes with them.
struct EngineRead {
    data: ActivityData,
    maps: ModelMaps,
    scaling: EngineScaling,
}

fn read_engine<E: BarrierIterate + ?Sized>(bs: &E) -> Result<EngineRead, String> {
    let data = bs.activity_data()?;
    data.check()?;
    let dims = bs.block_dims();
    if data.sigma_x.len() != dims[0] || data.sigma_s.len() != dims[1] {
        return Err(format!(
            "activity data: {} variables and {} inequality rows, but the KKT blocks hold {} and {}",
            data.sigma_x.len(),
            data.sigma_s.len(),
            dims[0],
            dims[1]
        ));
    }
    let maps = bs.model_maps()?;
    maps.check(&dims)?;
    let scaling = bs.scaling();
    scaling.check(&dims)?;
    Ok(EngineRead {
        data,
        maps,
        scaling,
    })
}

/// 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(data: &ActivityData, d_var: Option<&[Number]>, n: usize) -> VarFrame {
    let dv = |i: usize| -> Number { d_var.map_or(1.0, |d| d[i]) };
    let hess = &data.hessian;
    // `Σ̃_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 = data
        .sigma_x
        .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 }
}

/// Classifies every bounded variable and every finite-bounded inequality row
/// of the engine's converged point. See the module documentation for the
/// rule.
///
/// The slacks have to be measured against the bounds the model declares,
/// which means the engine must not have relaxed them (`bound_relax_factor =
/// 0` in Ipopt and pounce, where the default is `1e-8`). A relaxed bound
/// leaves an active slack near the relaxation rather than near `μ/z`, so the
/// complementarity products no longer track `μ` and the ratios land in the
/// wrong band. pounce's `Solver::classify_activity` refuses such a solve. The
/// trait cannot report a relaxation, so an engine outside pounce checks it
/// before calling this.
pub fn compute<E: BarrierIterate + ?Sized>(bs: &E) -> Result<ActivityReport, String> {
    let EngineRead {
        data,
        maps,
        scaling,
    } = read_engine(bs)?;
    let mu = bs.barrier_mu();
    let (n, m_d) = (data.sigma_x.len(), data.sigma_s.len());
    let (obj_scale, d_scale) = (scaling.objective, scaling.inequality_rows.as_deref());

    // 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 = scaling.variables.as_deref();
    let dv = |i: usize| -> Number { d_var.map_or(1.0, |d| d[i]) };

    // --- variables, in internal space ------------------------------------
    let has_l = present(&data.x_l.pos, n);
    let has_u = present(&data.x_u.pos, n);
    let z_l = expand(&data.x_l.mult, &data.x_l.pos, n);
    let z_u = expand(&data.x_u.mult, &data.x_u.pos, 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(&data.x_l.slack, &data.x_l.pos, n);
    let s_u = expand(&data.x_u.slack, &data.x_u.pos, 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 = &data.hessian;
    let VarFrame {
        sigma: sigma_x,
        diag,
        floor,
    } = var_frame(&data, d_var, 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(&data.d_l.pos, m_d);
    let rhas_u = present(&data.d_u.pos, m_d);
    let v_l = expand(&data.d_l.mult, &data.d_l.pos, m_d);
    let v_u = expand(&data.d_u.mult, &data.d_u.pos, m_d);
    let rs_l = expand(&data.d_l.slack, &data.d_l.pos, m_d);
    let rs_u = expand(&data.d_u.slack, &data.d_u.pos, m_d);
    let sigma_s = &data.sigma_s;

    let jac_d = &data.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.map_or(1.0, |v| v[j]);
        rows[j].sigma *= dg * dg / obj_scale;
    }

    // --- scatter to user space --------------------------------------------
    let n_full_x = maps.n_full_x;
    let n_full_g = maps.n_full_g();

    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[maps.var_x_to_full_x[i]] = *e;
    }

    let equality_entry = Entry {
        status: EQUALITY,
        ..NOT_CLASSIFIED
    };
    let mut row_full = vec![equality_entry; n_full_g];
    for (slot, d_pos) in row_full.iter_mut().zip(&maps.full_g_to_d) {
        if let Some(d) = d_pos {
            *slot = rows[*d];
        }
    }

    Ok(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.
#[derive(Debug, Clone, PartialEq)]
pub 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,
    /// The engine could not supply the converged quantities, or supplied
    /// inconsistent ones. The message says which.
    Data(String),
}

/// [`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 fn reduced_activity<E: BarrierIterate + ?Sized>(
    bs: &E,
    user_vars: &[usize],
) -> Result<ReducedActivityReport, ReducedActivityError> {
    let EngineRead {
        data,
        maps,
        scaling,
    } = read_engine(bs).map_err(ReducedActivityError::Data)?;
    let mu = bs.barrier_mu();
    let n = data.sigma_x.len();
    let full_x_to_var_x = maps.full_x_to_var_x();

    // 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 = maps.n_full_x;
    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(full_x_to_var_x[i]);
    }

    let obj_scale = scaling.objective;
    let has_l = present(&data.x_l.pos, n);
    let has_u = present(&data.x_u.pos, n);
    let frame = var_frame(&data, scaling.variables.as_deref(), 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.
#[derive(Debug, Clone, PartialEq)]
pub 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,
    /// The engine could not supply the converged quantities, or supplied
    /// inconsistent ones. The message says which.
    Data(String),
}

/// `‖∇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 fn reduced_row_activity<E: BarrierIterate + ?Sized>(
    bs: &E,
    user_rows: &[usize],
) -> Result<ReducedRowActivityReport, ReducedRowActivityError> {
    let EngineRead {
        data,
        maps,
        scaling,
    } = read_engine(bs).map_err(ReducedRowActivityError::Data)?;
    let mu = bs.barrier_mu();
    let (n, m_d) = (data.sigma_x.len(), data.sigma_s.len());

    // 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 = maps.n_full_g();
    let d_index = &maps.full_g_to_d;
    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 (obj_scale, d_scale) = (scaling.objective, scaling.inequality_rows.as_deref());
    let rhas_l = present(&data.d_l.pos, m_d);
    let rhas_u = present(&data.d_u.pos, m_d);
    let sigma_s = &data.sigma_s;

    // 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 = var_frame(&data, scaling.variables.as_deref(), n).floor;
    let norm2 = row_norm2(
        &data.jac_d,
        scaling.variables.as_deref(),
        &slot_of,
        n,
        m_d,
        n_slots,
    );
    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.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)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::backsolver::SensBacksolver;
    use crate::barrier::{BarrierIterate, BoundSide, KktOperator};
    use pounce_linalg::triplet::{GenTMatrix, GenTMatrixSpace, SymTMatrix, SymTMatrixSpace};
    use std::rc::Rc;

    const MU: Number = 1e-6;

    /// An engine that supplies the classification's data and nothing else.
    /// Two variables and two inequality rows, each with a lower bound at a
    /// unit slack. The first bound's multiplier puts `s·z` at `μ`, on the
    /// central path, and the second's at `100μ`, off it.
    struct TwoBoundsPerSide;

    fn side() -> BoundSide {
        BoundSide {
            pos: vec![0, 1],
            mult: vec![MU, 100.0 * MU],
            slack: vec![1.0, 1.0],
        }
    }

    impl SensBacksolver for TwoBoundsPerSide {
        fn dim(&self) -> usize {
            self.block_dims().iter().sum()
        }
        fn solve(&self, _rhs: &[Number], _lhs: &mut [Number]) -> bool {
            false
        }
    }

    impl BarrierIterate for TwoBoundsPerSide {
        fn block_dims(&self) -> [usize; 8] {
            [2, 2, 0, 2, 2, 0, 2, 0]
        }
        fn barrier_mu(&self) -> Number {
            MU
        }
        fn converged_iterate(&self, _out: &mut [Number]) -> Result<(), String> {
            Err("not used by the classification".into())
        }
        fn barrier_residual(&self, _point: &[Number], _out: &mut [Number]) -> Result<(), String> {
            Err("not used by the classification".into())
        }
        fn corrector_operator<'a>(
            &'a self,
            _predicted: &[Number],
            _pinned: &[(usize, Number)],
            _exact_hessian: bool,
        ) -> Result<Box<dyn KktOperator + 'a>, String> {
            Err("not used by the classification".into())
        }
        fn activity_data(&self) -> Result<ActivityData, String> {
            let mut hess = SymTMatrix::new(SymTMatrixSpace::new(2, vec![1, 2], vec![1, 2]));
            hess.set_values(&[1.0, 1.0]);
            let mut jac = GenTMatrix::new(GenTMatrixSpace::new(2, 2, vec![1, 2], vec![1, 2]));
            jac.set_values(&[1.0, 1.0]);
            Ok(ActivityData {
                x_l: side(),
                x_u: BoundSide::default(),
                d_l: side(),
                d_u: BoundSide::default(),
                sigma_x: vec![MU, 100.0 * MU],
                sigma_s: vec![MU, 100.0 * MU],
                hessian: Rc::new(hess),
                jac_c: Rc::new(GenTMatrix::new(GenTMatrixSpace::new(0, 2, vec![], vec![]))),
                jac_d: Rc::new(jac),
            })
        }
        fn model_maps(&self) -> Result<ModelMaps, String> {
            Ok(ModelMaps {
                n_full_x: 2,
                var_x_to_full_x: vec![0, 1],
                full_g_to_c: vec![None, None],
                full_g_to_d: vec![Some(0), Some(1)],
            })
        }
        fn scaling(&self) -> EngineScaling {
            EngineScaling::none()
        }
        fn primal_bounds(&self) -> Result<(Vec<Number>, Vec<Number>), String> {
            Err("not used by the classification".into())
        }
    }

    /// Each bound's central-path check reads its own multiplier and slack,
    /// not those of another bound on the same side.
    #[test]
    fn each_bound_is_checked_against_the_central_path_with_its_own_product() {
        let report = compute(&TwoBoundsPerSide).expect("the classification");
        assert_eq!(report.var_off_central_path, vec![false, true]);
        assert_eq!(report.row_off_central_path, vec![false, true]);
    }
}