pounce-sens-core 0.12.0

//! The activity classification rule, independent of which solver produced the
//! numbers it reads.
//!
//! A converged interior-point solve leaves, for each bounded variable or row, a
//! barrier diagonal `Σ` and a curvature `q`. Their ratio `r = Σ/q` says what the
//! bound is doing: `O(μ)` and it is inactive, `O(1/μ)` and it holds the
//! coordinate, `O(1)` and slack and multiplier are vanishing together — a
//! *kink*, where the derivative is two-valued and the parametric step needs a
//! side.
//!
//! Only the decision lives here. Deriving `Σ`, `q` and `μ` is the caller's, and
//! is where the two arms genuinely differ: the NLP arm reads them off the
//! filter-IPM's eight-block iterate, the convex arm derives them from
//! `(problem, solution)`. Sharing the *rule* is what stops the two drifting into
//! disagreeing about what a kink is; sharing the plumbing would mean abstracting
//! `IpoptData` access behind another trait, which is a much larger change.
//!
//! # `AMBIGUOUS` is not "probably not a kink"
//!
//! The single most important thing to know before reading a status. A genuine
//! kink lands in [`AMBIGUOUS`] whenever its coordinate is *coupled*, because the
//! `q` a cheap classifier has to hand is the Hessian's **diagonal**, while the
//! multiplier is generated by the curvature *reduced* along that coordinate. The
//! ratio is then `reduced/diagonal`, which is 1 only when the coordinate is
//! decoupled — routine on a collocation model to be otherwise.
//!
//! That misreading shipped a defect once already (gh#763), so it is stated at
//! the rule rather than only at the call sites: **never use the activity class
//! as a proxy for kink-ness.** And note the failure is not one a tighter solve
//! fixes — `reduced/diagonal` is μ-independent, so re-solving does not separate
//! it. The refinement is a reduced-curvature normalizer, computed on demand over
//! the ambiguous entries because it costs a back-solve each.
//!
//! # The `floor` is the caller's
//!
//! [`classify_entry`] takes the noise floor below which `q` is not a curvature
//! at all as a parameter. That is the seam: the floor's derivation
//! (`√ε · max|diag|`, say) needs the whole operator, which is exactly the part
//! that does not generalize.

use pounce_common::types::Number;

/// No finite bound on this variable or row: nothing to classify.
pub const UNBOUNDED: i8 = -1;
/// `r = O(μ)`: the bound is not doing anything.
pub const INACTIVE: i8 = 0;
/// `r = O(1)`: slack and multiplier vanish together; kept, flagged.
pub const WEAKLY_ACTIVE: i8 = 1;
/// `r = O(1/μ)`: the bound holds the variable; projected out.
pub const STRONGLY_ACTIVE: i8 = 2;
/// `r` in a gap between the band and a `μ`-edge: undetermined at this `μ`.
///
/// **Not "probably not a kink"** — see the module documentation. A coupled
/// kink lands here by construction, and re-solving tighter does not move it.
pub const AMBIGUOUS: i8 = 3;
/// The curvature `q` is below noise scale: the bound question does not arise,
/// and the direction is poorly identified.
pub const UNIDENTIFIED: i8 = 4;
/// `lb == ub`: the variable was removed from the solve as a parameter, so there
/// is no barrier geometry to classify.
pub const FIXED: i8 = 5;
/// An equality constraint: always active by construction, with no slack or
/// multiplier pair on the barrier, so outside this classification.
pub const EQUALITY: i8 = 6;

/// The classification rule.
///
/// The band is fixed at `[1e-1, 1e1]` while the μ-edges `√μ` and `1/√μ` move
/// with the solve: they meet the band at `μ = 1e-2`, and a full decade separates
/// them from it at `μ = 1e-4`. Above `1e-4` that margin is what is thinning, so
/// only the two calls that stay clear are made and the middle is an honest
/// refusal rather than a guess.
pub fn classify(r: Number, mu: Number) -> i8 {
    if mu > 1e-4 {
        if r < 1e-1 {
            INACTIVE
        } else if r > 1e1 {
            STRONGLY_ACTIVE
        } else {
            AMBIGUOUS
        }
    } else if r < mu.sqrt() {
        INACTIVE
    } else if r > 1.0 / mu.sqrt() {
        STRONGLY_ACTIVE
    } else if (1e-1..=1e1).contains(&r) {
        WEAKLY_ACTIVE
    } else {
        AMBIGUOUS
    }
}

/// `-1`, `0` or `+1`.
pub fn sign_of(x: Number) -> i8 {
    if x > 0.0 {
        1
    } else if x < 0.0 {
        -1
    } else {
        0
    }
}

/// Central-path check for one side: `s·z` within a factor of ten of `μ`.
pub fn off_path(s: Number, z: Number, mu: Number) -> bool {
    let comp = s * z;
    comp > 10.0 * mu || comp < 0.1 * mu
}

/// Classified inactive yet `r` well above the `O(μ)` an inactive bound should
/// carry: barrier curvature where none should be.
///
/// The threshold is μ-relative because *inactive* MEANS `r = O(μ)`; a fixed
/// constant can never sit below the inactive edge `√μ` at any converged μ.
pub fn contaminated(status: i8, r: Number, mu: Number) -> bool {
    status == INACTIVE && r > 100.0 * mu
}

/// One classified entry, before any user-space scatter.
#[derive(Clone, Copy, Debug)]
pub struct Entry {
    pub status: i8,
    pub ratio: Number,
    pub q_sign: i8,
    pub off_path: bool,
    pub contaminated: bool,
    /// The RAW barrier diagonal, whatever weight classification used.
    pub sigma: Number,
}

/// The entry for something that was never classified.
pub const NOT_CLASSIFIED: Entry = Entry {
    status: UNBOUNDED,
    ratio: Number::NAN,
    q_sign: 0,
    off_path: false,
    contaminated: false,
    sigma: 0.0,
};

/// A bounded row whose gradient vanishes at the iterate has no direction to
/// measure curvature along: unidentified, exactly as a below-floor `q`, never
/// `UNBOUNDED` (the bounds are real). The ratio is the raw `Σ/floor` lower
/// bound; the geometric weight is degenerate at zero gradient.
pub fn zero_gradient_row(sigma: Number, floor: Number) -> Entry {
    Entry {
        status: UNIDENTIFIED,
        ratio: sigma / floor,
        q_sign: 0,
        off_path: false,
        contaminated: false,
        sigma,
    }
}

/// Classify one bounded variable or row from its `Σ` and signed `q`.
///
/// `off_path` is the caller's to fill: it reads the per-side slack and
/// multiplier, not the ratio.
pub fn classify_entry(sigma: Number, q_signed: Number, floor: Number, mu: Number) -> Entry {
    let q_sign = sign_of(q_signed);
    let q = q_signed.abs();
    if q < floor {
        return Entry {
            status: UNIDENTIFIED,
            ratio: sigma / floor,
            q_sign,
            off_path: false,
            contaminated: false,
            sigma,
        };
    }
    let r = sigma / q;
    let status = classify(r, mu);
    Entry {
        status,
        ratio: r,
        q_sign,
        off_path: false,
        contaminated: contaminated(status, r, mu),
        sigma,
    }
}

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

    #[test]
    fn tight_mu_walks_all_five_regions() {
        let mu = 1e-10; // edges at 1e-5 and 1e5
        assert_eq!(classify(0.9e-5, mu), INACTIVE);
        assert_eq!(classify(1.1e-5, mu), AMBIGUOUS); // gap: edge..band
        assert_eq!(classify(0.5, mu), WEAKLY_ACTIVE);
        assert_eq!(classify(50.0, mu), AMBIGUOUS); // gap: band..edge
        assert_eq!(classify(2e5, mu), STRONGLY_ACTIVE);
    }

    #[test]
    fn band_edges_are_inclusive_and_mu_edges_separate() {
        let mu = 1e-10;
        assert_eq!(classify(1e-1, mu), WEAKLY_ACTIVE);
        assert_eq!(classify(1e1, mu), WEAKLY_ACTIVE);
        // either side of each μ-edge (exactly-on is float-fragile:
        // √(1e-10) is not exactly 1e-5)
        assert_eq!(classify(0.99e-5, mu), INACTIVE);
        assert_eq!(classify(1.01e-5, mu), AMBIGUOUS);
        assert_eq!(classify(0.99e5, mu), AMBIGUOUS);
        assert_eq!(classify(1.01e5, mu), STRONGLY_ACTIVE);
    }

    #[test]
    fn loose_mu_refuses_the_weak_call() {
        // μ > 1e-4: three statuses only, the band reports ambiguous
        for mu in [1e-3, 1e-2, 1e-1] {
            assert_eq!(classify(0.05, mu), INACTIVE);
            assert_eq!(classify(1.0, mu), AMBIGUOUS);
            assert_eq!(classify(50.0, mu), STRONGLY_ACTIVE);
        }
        // at μ = 1e-4 exactly the μ-branch is not taken: the weak call
        // is available, with a decade of margin edge-to-band
        assert_eq!(classify(1.0, 1e-4), WEAKLY_ACTIVE);
    }

    #[test]
    fn off_path_is_a_factor_of_ten_both_ways() {
        let mu = 1e-2;
        assert!(!off_path(1.0, 1e-2, mu)); // s·z = μ exactly
        assert!(!off_path(0.5, 1e-2, mu)); // within 10×
        assert!(off_path(1.0, 0.2, mu)); // 20× above
        assert!(off_path(1.0, 5e-4, mu)); // 20× below
    }

    #[test]
    fn contamination_is_mu_relative_and_inactive_only() {
        let mu = 1e-10; // inactive edge at 1e-5, threshold at 1e-8
        assert!(contaminated(INACTIVE, 1e-6, mu));
        assert!(!contaminated(INACTIVE, 5e-9, mu));
        assert!(!contaminated(WEAKLY_ACTIVE, 1.0, mu));
        assert!(!contaminated(STRONGLY_ACTIVE, 1e5, mu));
        // the flag is reachable: 100μ sits below the inactive edge √μ
        // whenever μ < 1e-4, so an inactive r can exceed it
        assert!(100.0 * mu < mu.sqrt());
    }

    #[test]
    fn below_floor_reports_unidentified_with_the_sign() {
        let e = classify_entry(0.5, 1e-12, 1e-8, 1e-10);
        assert_eq!(e.status, UNIDENTIFIED);
        assert_eq!(e.q_sign, 1);
        let e = classify_entry(0.5, -1e-12, 1e-8, 1e-10);
        assert_eq!(e.status, UNIDENTIFIED);
        assert_eq!(e.q_sign, -1);
        // negative curvature above the floor classifies on |q| but
        // keeps its sign visible
        let e = classify_entry(1.0, -2.0, 1e-8, 1e-10);
        assert_eq!(e.status, WEAKLY_ACTIVE);
        assert_eq!(e.q_sign, -1);
        assert_eq!(e.ratio, 0.5);
    }
}