pounce_algorithm/mu/oracle/
probing.rs

1//! Probing (Mehrotra) mu oracle — port of
2//! `IpProbingMuOracle.{hpp,cpp}`. Phase 10.
3//!
4//! Given the affine-step complementarity `mu_aff` and current `mu`,
5//! the probing rule chooses
6//!
7//! ```text
8//!   sigma = min((mu_aff / mu)^3, sigma_max)
9//!   mu_new = sigma * mu
10//! ```
11//!
12//! per `IpProbingMuOracle.cpp:117-121` and the standard Mehrotra
13//! predictor-corrector formula.
14
15use crate::ipopt_cq::IpoptCqHandle;
16use crate::ipopt_data::IpoptDataHandle;
17use crate::ipopt_nlp::IpoptNlp;
18use crate::iterates_vector::IteratesVector;
19use crate::kkt::pd_search_dir_calc::PdSearchDirCalc;
20use crate::mu::oracle::r#trait::MuOracle;
21use pounce_common::types::Number;
22use std::cell::RefCell;
23use std::rc::Rc;
24
25pub struct ProbingMuOracle {
26    pub sigma_max: Number,
27    pub mu_min: Number,
28    pub mu_max: Number,
29    pub mu_curr: Number,
30    pub mu_aff: Number,
31}
32
33impl Default for ProbingMuOracle {
34    fn default() -> Self {
35        Self {
36            // Default `sigma_max = 100` per `RegisterOptions`.
37            sigma_max: 100.0,
38            mu_min: 1e-11,
39            mu_max: 1e5,
40            mu_curr: 1.0,
41            mu_aff: 1.0,
42        }
43    }
44}
45
46impl ProbingMuOracle {
47    pub fn new() -> Self {
48        Self::default()
49    }
50
51    /// Pure-arithmetic probing formula from `IpProbingMuOracle.cpp:117`.
52    pub fn probing_mu(mu_curr: Number, mu_aff: Number, sigma_max: Number) -> Number {
53        let sigma = (mu_aff / mu_curr).powi(3).min(sigma_max);
54        sigma * mu_curr
55    }
56
57    /// Drive the affine (predictor) step through `pd_search_dir`,
58    /// snapshot `μ_aff` from the projected complementarity, then return
59    /// the probed `μ_new ∈ [mu_min, mu_max]`. Returns `None` if the
60    /// affine linear solve fails (caller falls back to LOQO).
61    ///
62    /// `tau` is the fraction-to-the-boundary parameter (typically
63    /// `1.0` for the affine probe per upstream's
64    /// `IpProbingMuOracle.cpp:107`).
65    pub fn calculate_mu_with_affine_step(
66        &mut self,
67        data: &IpoptDataHandle,
68        cq: &IpoptCqHandle,
69        nlp: &Rc<RefCell<dyn IpoptNlp>>,
70        pd_search_dir: &mut PdSearchDirCalc,
71        tau: Number,
72    ) -> Option<Number> {
73        if !pd_search_dir.compute_affine_step(data, cq, nlp) {
74            return None;
75        }
76        let delta_aff: IteratesVector = data.borrow().delta_aff.clone()?;
77        let cq_ref = cq.borrow();
78        // Upstream `IpProbingMuOracle.cpp:111` uses `curr_avrg_compl()`
79        // (the actual average s·z at the current iterate), NOT the
80        // stored barrier parameter `data.curr_mu`. Mehrotra's rule
81        // sigma = (mu_aff / mu_curr)^3 only makes sense when mu_curr
82        // is the present complementarity — using the barrier value
83        // produces sigma values orders of magnitude off and causes
84        // huge over-correction steps on iter 1 of LP-shaped problems.
85        let mu_curr = cq_ref.curr_avrg_compl();
86        let alpha_pri = cq_ref.aff_step_alpha_primal_max(&delta_aff, tau);
87        let alpha_du = cq_ref.aff_step_alpha_dual_max(&delta_aff, tau);
88        let mu_aff = cq_ref.aff_step_compl_avrg(&delta_aff, alpha_pri, alpha_du);
89        self.mu_curr = mu_curr;
90        self.mu_aff = mu_aff;
91        let raw = Self::probing_mu(mu_curr, mu_aff, self.sigma_max);
92        Some(raw.clamp(self.mu_min, self.mu_max))
93    }
94}
95
96impl MuOracle for ProbingMuOracle {
97    fn calculate_mu(&mut self) -> Option<Number> {
98        let raw = Self::probing_mu(self.mu_curr, self.mu_aff, self.sigma_max);
99        Some(raw.clamp(self.mu_min, self.mu_max))
100    }
101}
102
103#[cfg(test)]
104mod tests {
105    use super::*;
106
107    #[test]
108    fn probing_aff_equals_curr_keeps_mu() {
109        // ratio=1 → sigma=1 → mu_new=mu_curr.
110        assert_eq!(ProbingMuOracle::probing_mu(1.0, 1.0, 100.0), 1.0);
111    }
112
113    #[test]
114    fn probing_aff_half_curr() {
115        // ratio=0.5 → sigma=0.125 → mu_new=0.125 * 1.0.
116        let m = ProbingMuOracle::probing_mu(1.0, 0.5, 100.0);
117        assert!((m - 0.125).abs() < 1e-15);
118    }
119
120    #[test]
121    fn probing_caps_at_sigma_max() {
122        // ratio=10 → sigma=1000, capped at 100. mu_new=100.
123        let m = ProbingMuOracle::probing_mu(1.0, 10.0, 100.0);
124        assert!((m - 100.0).abs() < 1e-13);
125    }
126
127    #[test]
128    fn calculate_mu_via_trait_clamped() {
129        let mut o = ProbingMuOracle {
130            sigma_max: 100.0,
131            mu_min: 0.5,
132            mu_max: 10.0,
133            mu_curr: 1.0,
134            mu_aff: 0.001, // sigma ≈ 1e-9, mu_new ≈ 1e-9, clamps to 0.5.
135        };
136        assert_eq!(o.calculate_mu(), Some(0.5));
137    }
138}
pounce_algorithm/mu/oracle/probing.rs

pounce_algorithm/mu/oracle/
probing.rs
