pounce-sens-core 0.12.0

//! `ReducedHessianCalculator` — port of upstream
//! [`SensReducedHessianCalculator.{hpp,cpp}`](../../../ref/Ipopt/contrib/sIPOPT/src/SensReducedHessianCalculator.cpp).
//!
//! # What this computes
//!
//! Given a converged KKT factor `K` and a parameter-row selector `B`
//! that picks out the **free variables** (post active-set
//! elimination), the reduced Hessian is
//!
//! ```text
//! H_R = obj_scal · B · K⁻¹ · Bᵀ
//! ```
//!
//! per upstream's sign + obj-scaling convention at
//! [`SensReducedHessianCalculator.cpp:90-97`](../../../ref/Ipopt/contrib/sIPOPT/src/SensReducedHessianCalculator.cpp):
//! the raw Schur output is `S = -B K⁻¹ Bᵀ` (with the leading minus
//! from the augmented-system reduction), which is then multiplied
//! by `-obj_scal` to produce the unscaled reduced Hessian.
//!
//! In pounce we default `obj_scal = 1.0` so the operation reduces to
//! `H_R = -S = B K⁻¹ Bᵀ`.
//!
//! # `B K⁻¹ Bᵀ` is `H_R` or `−H_R` depending on which rows `B` picks
//!
//! This function is agnostic about that; a caller must not be
//! (gh#937). Upstream — and pounce's CLI `red_hessian` suffix path,
//! which mirrors it — selects **free-variable (x-block) rows**, where
//! the corresponding block of `K⁻¹` is itself an inverse. The
//! `pounce-sensitivity` pin path instead selects **`y_c` multiplier
//! rows**, and for `K = [[H, Aᵀ], [A, 0]]` that block of `K⁻¹` is
//! `−(A H⁻¹ Aᵀ)⁻¹` — an inverse of an inverse, so the result is
//! `±H_R` itself rather than a submatrix of `K⁻¹`, and the `-S` above
//! leaves it at `−H_R`.
//!
//! The two blocks sit on opposite sides of one inversion, so the unit
//! test below — which feeds a synthetic dense `K` and selects two of
//! its rows — looks like it generalizes to the pin path and does not.
//! `pounce-sensitivity/tests/issue_937_reduced_hessian_sign.rs` pins
//! the pin path's orientation at the public API, on a model where
//! `±H_R` and `H_R⁻¹` differ in magnitude as well as sign.
//!
//! Unlike upstream, no NLP-side scaling needs
//! folding in here: since pounce#128 the live-factor backsolver
//! (`PdSensBacksolver`, in `pounce-sensitivity` — a downstream crate, so not
//! linkable from here) conjugates every back-solve by the
//! NLP scaling diagonal, so `K⁻¹` is already the natural-units KKT
//! inverse and `obj_scal` survives purely as a user-side extra
//! multiplier.
//!
//! Reference: Pirnay, López-Negrete & Biegler 2012, §5
//! (reduced-Hessian use case), DOI:
//! [10.1007/s12532-012-0043-2](https://doi.org/10.1007/s12532-012-0043-2).

use crate::p_calculator::PCalculator;
use crate::schur_data::SchurData;
use pounce_common::types::Number;

/// Compute the (column-major) reduced Hessian into the caller-supplied
/// buffer.
///
/// Mirrors `ReducedHessianCalculator::ComputeReducedHessian`
/// ([`SensReducedHessianCalculator.cpp:42-113`](../../../ref/Ipopt/contrib/sIPOPT/src/SensReducedHessianCalculator.cpp)).
///
/// # Arguments
///
/// - `pcalc`: a P-calculator that was built with `data_A = hess_data`
///   (the **same** SchurData on both sides — the reduced-Hessian
///   computation is the diagonal `B = A` case of `schur_matrix`).
///   Must have completed `compute_p` before this call; this function
///   calls `schur_matrix(hess_data, …)` which will run `compute_p`
///   lazily if needed.
/// - `hess_data`: the free-variable selector, in the IndexSchurData
///   format. Conceptually the same matrix as `pcalc.data_a()`.
/// - `obj_scal`: per-NLP objective scaling factor; pounce defaults to
///   `1.0` (no scaling). Mirrors upstream's `apply_obj_scaling(1.0)`
///   ([`SensReducedHessianCalculator.cpp:91`](../../../ref/Ipopt/contrib/sIPOPT/src/SensReducedHessianCalculator.cpp)).
/// - `out`: caller-allocated buffer of length `n_rows × n_rows` where
///   `n_rows = hess_data.nrows()`. Column-major. Overwritten on
///   success.
///
/// Returns `false` if the underlying `schur_matrix` call fails or
/// `out` is mis-sized.
///
/// # Note on scaling-induced warnings
///
/// Upstream prints a J_WARNING block when any of `x` / `c` / `d`
/// scaling is active because "a correct unscaled solution of the
/// reduced hessian cannot be guaranteed"
/// ([`SensReducedHessianCalculator.cpp:64-88`](../../../ref/Ipopt/contrib/sIPOPT/src/SensReducedHessianCalculator.cpp)).
/// Pounce's Phase-C surface takes `obj_scal` as an explicit argument
/// rather than reaching into an `NLP_scaling()` object; users
/// constructing a reduced Hessian outside a configured pounce IPM
/// own the responsibility to pass an `obj_scal` consistent with
/// whatever scaling their `K` factor encodes.
/// Convert a natural-units reduced Hessian (column-major `n×n`,
/// `n = dc.len()`) into the solver's internal scaled space in place:
/// `H̃_ij = (df / (dc_i·dc_j)) · H_ij`, with `df` the effective
/// objective scaling factor and `dc` the pin rows' constraint scaling
/// factors. This is the inverse of the natural-units correction the
/// live-factor backsolver applies (pounce#128) — i.e. the value
/// pounce returned before #128. Single home for the formula so the
/// `SensSolve` and `Solver` surfaces cannot drift.
///
/// Returns `false` when `hr` is mis-sized.
pub fn scale_to_solver_space(hr: &mut [Number], df: Number, dc: &[Number]) -> bool {
    let n = dc.len();
    if hr.len() != n * n {
        return false;
    }
    for j in 0..n {
        for i in 0..n {
            hr[j * n + i] *= df / (dc[i] * dc[j]);
        }
    }
    true
}

pub fn compute_reduced_hessian<P: PCalculator>(
    pcalc: &mut P,
    hess_data: &dyn SchurData,
    obj_scal: Number,
    out: &mut [Number],
) -> bool {
    let n = hess_data.nrows() as usize;
    if out.len() != n * n {
        return false;
    }
    // Step 1: build the Schur matrix S = -B K⁻¹ Bᵀ via the PCalculator.
    if !pcalc.schur_matrix(hess_data, out) {
        return false;
    }
    // Step 2: H_R = -obj_scal · S. Pounce default obj_scal = 1.0
    // yields H_R = B K⁻¹ Bᵀ — whose own sign depends on which rows
    // `hess_data` selects; see the module header (gh#937).
    let factor = -obj_scal;
    for v in out.iter_mut() {
        *v *= factor;
    }
    true
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::backsolver::DenseLuBacksolver;
    use crate::p_calculator::IndexPCalculator;
    use crate::schur_data::IndexSchurData;

    /// 3×3 SPD `K = [[2,-1,0], [-1,2,-1], [0,-1,2]]`. Inverse is
    /// `K⁻¹ = 1/4 · [[3, 2, 1], [2, 4, 2], [1, 2, 3]]`.
    ///
    /// Pick `hess_data` selecting rows/cols {0, 2}. The 2×2
    /// submatrix of `K⁻¹` at indices {0, 2} is
    /// `[[3/4, 1/4], [1/4, 3/4]]`. With `obj_scal = 1.0`, that's
    /// the reduced Hessian.
    #[test]
    fn reduced_hessian_matches_kinv_submatrix() {
        #[rustfmt::skip]
        let k = vec![
             2.0, -1.0,  0.0,
            -1.0,  2.0, -1.0,
             0.0, -1.0,  2.0,
        ];
        let backsolver = DenseLuBacksolver::from_dense(3, &k).unwrap();
        let hess_data = IndexSchurData::from_parts(vec![0, 2], vec![1, 1]).unwrap();
        let pcalc_a = IndexSchurData::from_parts(vec![0, 2], vec![1, 1]).unwrap();
        let mut pcalc = IndexPCalculator::new(backsolver, pcalc_a);

        let mut hr = vec![0.0; 2 * 2];
        assert!(compute_reduced_hessian(
            &mut pcalc, &hess_data, 1.0, &mut hr
        ));

        // Column-major: [j * n + i]
        // H_R[0,0] = 3/4
        // H_R[1,0] = 1/4
        // H_R[0,1] = 1/4
        // H_R[1,1] = 3/4
        assert!((hr[0] - 0.75).abs() < 1e-12, "H_R[0,0] = {}", hr[0]);
        assert!((hr[1] - 0.25).abs() < 1e-12, "H_R[1,0] = {}", hr[1]);
        assert!((hr[2] - 0.25).abs() < 1e-12, "H_R[0,1] = {}", hr[2]);
        assert!((hr[3] - 0.75).abs() < 1e-12, "H_R[1,1] = {}", hr[3]);
    }

    #[test]
    fn reduced_hessian_applies_obj_scal() {
        #[rustfmt::skip]
        let k = vec![
             2.0, -1.0,  0.0,
            -1.0,  2.0, -1.0,
             0.0, -1.0,  2.0,
        ];
        let backsolver = DenseLuBacksolver::from_dense(3, &k).unwrap();
        let hess_data = IndexSchurData::from_parts(vec![1], vec![1]).unwrap();
        let pcalc_a = IndexSchurData::from_parts(vec![1], vec![1]).unwrap();
        let mut pcalc = IndexPCalculator::new(backsolver, pcalc_a);

        // K⁻¹[1, 1] = 1, so H_R = obj_scal · 1.
        let mut hr = vec![0.0; 1];
        assert!(compute_reduced_hessian(
            &mut pcalc, &hess_data, 2.5, &mut hr
        ));
        assert!((hr[0] - 2.5).abs() < 1e-12, "H_R = {}", hr[0]);
    }

    #[test]
    fn reduced_hessian_rejects_wrong_buffer_size() {
        let backsolver = DenseLuBacksolver::from_dense(2, &[1.0, 0.0, 0.0, 1.0]).unwrap();
        let hd = IndexSchurData::from_parts(vec![0, 1], vec![1, 1]).unwrap();
        let pa = IndexSchurData::from_parts(vec![0, 1], vec![1, 1]).unwrap();
        let mut pc = IndexPCalculator::new(backsolver, pa);
        // hd is 2 rows so output should be 4 entries; pass 3.
        let mut wrong = vec![0.0; 3];
        assert!(!compute_reduced_hessian(&mut pc, &hd, 1.0, &mut wrong));
    }
}