Skip to main content

pounce_algorithm/mu/
adaptive.rs

1//! Adaptive mu update — port of `IpAdaptiveMuUpdate.{hpp,cpp}`.
2//!
3//! Phase 10. The full update reaches into `IpoptCq` for residuals and
4//! into a `MuOracle` for the candidate σ; this file ships:
5//!
6//! * the option struct with upstream defaults from `RegisterOptions`,
7//! * the `lower_mu_safeguard` scalar core (lines 753-786),
8//! * the globalization-mode enum and the FreeMuMode/FixedMuMode state
9//!   machine (`UpdateBarrierParameter` lines 252-444),
10//! * the `mu_oracle` selector ([`MuOracleKind`]) — `Loqo` runs the
11//!   closed form; `Probing` / `QualityFunction` drive an affine /
12//!   centring solve when [`MuUpdate`] is given the search-dir + nlp
13//!   handles, otherwise fall through to LOQO (mirrors upstream's
14//!   "oracle returned no candidate" branch at lines 402-408).
15
16use crate::ipopt_cq::IpoptCqHandle;
17use crate::ipopt_data::IpoptDataHandle;
18use crate::ipopt_nlp::IpoptNlp;
19use crate::iterates_vector::IteratesVector;
20use crate::kkt::pd_search_dir_calc::PdSearchDirCalc;
21use crate::line_search::filter::Filter;
22use crate::mu::oracle::loqo::LoqoMuOracle;
23use crate::mu::oracle::probing::ProbingMuOracle;
24use crate::mu::oracle::quality_function::QualityFunctionMuOracle;
25use crate::mu::oracle::r#trait::MuOracle;
26use crate::mu::r#trait::MuUpdate;
27use pounce_common::types::Number;
28use std::cell::RefCell;
29use std::collections::VecDeque;
30use std::rc::Rc;
31
32/// `mu_oracle` option from `IpAdaptiveMuUpdate.cpp:RegisterOptions`.
33/// Default `QualityFunction` matches upstream (`"quality-function"`).
34#[derive(Debug, Clone, Copy, PartialEq, Eq)]
35pub enum MuOracleKind {
36    /// Closed-form LOQO rule. No predictor solve required.
37    Loqo,
38    /// Mehrotra probing oracle. Needs an affine-step solve.
39    Probing,
40    /// Golden-section minimisation of the q(σ) quality function.
41    /// Needs an affine-step solve plus a centring evaluator.
42    QualityFunction,
43}
44
45#[derive(Debug, Clone, Copy, PartialEq, Eq)]
46pub enum AdaptiveMuGlobalization {
47    KktError,
48    ObjConstrFilter,
49    NeverMonotoneMode,
50}
51
52#[derive(Debug, Clone, Copy, PartialEq, Eq)]
53pub enum AdaptiveMuKktNorm {
54    OneNorm,
55    TwoNormSquared,
56    MaxNorm,
57    TwoNorm,
58}
59
60pub struct AdaptiveMuUpdate {
61    pub mu_oracle: MuOracleKind,
62    pub adaptive_mu_globalization: AdaptiveMuGlobalization,
63    pub adaptive_mu_kkt_norm: AdaptiveMuKktNorm,
64    pub adaptive_mu_safeguard_factor: Number,
65    pub adaptive_mu_kkterror_red_iters: usize,
66    pub adaptive_mu_kkterror_red_fact: Number,
67    pub filter_max_margin: Number,
68    pub filter_margin_fact: Number,
69    pub mu_min: Number,
70    /// Upper bound on μ. Sentinel `-1.0` means "not yet computed; init
71    /// lazily on the first `update_barrier_parameter` call to
72    /// `mu_max_fact * curr_avrg_compl()`". Mirrors
73    /// `IpAdaptiveMuUpdate.cpp:160-165` (load step) and
74    /// `IpAdaptiveMuUpdate.cpp:267-274` (lazy init).
75    pub mu_max: Number,
76    /// `mu_max_fact` (default 1e3) — factor for lazy init of `mu_max`.
77    /// Upstream `IpAdaptiveMuUpdate.cpp:RegisterOptions` line 42.
78    /// Ignored if the user explicitly sets `mu_max` to a non-sentinel
79    /// value.
80    pub mu_max_fact: Number,
81    /// `tau_min` from `IpAdaptiveMuUpdate.cpp:RegisterOptions`. Used to
82    /// derive `curr_tau = max(tau_min, 1 - mu)` after each update,
83    /// mirroring upstream's `IpAdaptiveMuUpdate.cpp:UpdateBarrierParameter`
84    /// at the post-oracle update.
85    pub tau_min: Number,
86    /// Initial mu seed — `mu_init` from `IpoptAlgorithm` registered
87    /// options. Used to seed `curr_mu` in `initialize`.
88    pub mu_init: Number,
89    /// `barrier_tol_factor` (default 10) from upstream
90    /// `IpMonotoneMuUpdate::RegisterOptions`. Threshold for fixed-mode
91    /// barrier subproblem completion: reduce μ when
92    /// `curr_barrier_error ≤ barrier_tol_factor · μ`.
93    pub barrier_tol_factor: Number,
94    /// `mu_linear_decrease_factor` (default 0.2) — fixed-mode update
95    /// uses `min(linear · μ, μ^superlinear_power)`.
96    pub mu_linear_decrease_factor: Number,
97    /// `mu_superlinear_decrease_power` (default 1.5).
98    pub mu_superlinear_decrease_power: Number,
99    /// `adaptive_mu_monotone_init_factor` (default 0.8). Used by
100    /// `new_fixed_mu` when no `fix_mu_oracle_` is configured.
101    pub adaptive_mu_monotone_init_factor: Number,
102    /// `adaptive_mu_restore_previous_iterate` (default false).
103    pub restore_accepted_iterate: bool,
104    /// `sigma_max` / `sigma_min` forwarded to
105    /// `QualityFunctionMuOracle` on every free-mode call. Defaults
106    /// from `IpQualityFunctionMuOracle.cpp:RegisterOptions`.
107    pub sigma_max: Number,
108    pub sigma_min: Number,
109    /// `quality_function_norm_type` (default `2-norm-squared`) —
110    /// norm used to aggregate the three KKT components inside the
111    /// quality function. Forwarded to `QualityFunctionMuOracle` on
112    /// every free-mode call. Mirrors
113    /// `IpQualityFunctionMuOracle.cpp:RegisterOptions`.
114    pub qf_norm_type: crate::mu::oracle::quality_function::NormType,
115    /// `quality_function_centrality` (default `none`) — penalty term
116    /// added to the quality function for centrality deviation.
117    pub qf_centrality_type: crate::mu::oracle::quality_function::CentralityType,
118    /// `quality_function_balancing_term` (default `none`) — penalty
119    /// term added to the quality function when the complementarity
120    /// is far smaller than the infeasibilities.
121    pub qf_balancing_term: crate::mu::oracle::quality_function::BalancingTermType,
122    /// `quality_function_max_section_steps` (default 8) — cap on
123    /// golden-section iterations when picking σ.
124    pub qf_max_section_steps: i32,
125    /// `quality_function_section_sigma_tol` (default 1e-2) — width
126    /// tolerance in σ-space for the golden-section search.
127    pub qf_section_sigma_tol: Number,
128    /// `quality_function_section_qf_tol` (default 0.0) — relative
129    /// flatness tolerance for the golden-section search.
130    pub qf_section_qf_tol: Number,
131
132    /// `probing_iterate_quality_factor` (default 1e4, pounce-specific;
133    /// see pounce#58). When the probing (Mehrotra) μ-oracle is about
134    /// to read `curr_avrg_compl()` for its `mu_curr` input, a single
135    /// imbalanced `(s_i, z_i)` pair can inflate the average 5+ orders
136    /// above the stored `data.curr_mu`. Probing then mathematically
137    /// correctly returns `σ·mu_curr` ≫ previous μ, which throws the
138    /// iterate out of the convergence neighborhood. This guard
139    /// short-circuits that case: when `curr_avrg_compl / curr_mu >
140    /// probing_iterate_quality_factor`, we signal restoration via
141    /// [`IpoptData::request_resto`] and keep μ unchanged. Set to 0 or
142    /// any non-positive value to disable.
143    pub probing_iterate_quality_factor: Number,
144
145    /// Upstream tracks `init_*_inf` lazily — sentinel −1 means
146    /// "not yet captured".
147    init_dual_inf: Number,
148    init_primal_inf: Number,
149
150    /// FreeMuMode/FixedMuMode flag — port of
151    /// `IpoptData::FreeMuMode()`. `true` means "let the oracle drive
152    /// μ"; `false` means "monotone decrease until sufficient progress
153    /// is made". Initialised to `true` in [`MuUpdate::initialize`]
154    /// (matches upstream `InitializeImpl` line 239).
155    free_mu_mode: bool,
156    /// KKT-error history for `KKT_ERROR` globalization. Bounded length
157    /// = `adaptive_mu_kkterror_red_iters`. Mirrors `refs_vals_`.
158    refs_vals: VecDeque<Number>,
159    /// 2-D `(theta, phi)` filter for `OBJ_CONSTR_FILTER` globalization.
160    /// Mirrors `filter_` (constructed with `Filter(2)`).
161    filter: Filter,
162    /// Snapshot of `curr` at the most recent successful free-mode
163    /// iterate; restored when switching to fixed mode if
164    /// `restore_accepted_iterate` is on. Mirrors `accepted_point_`.
165    accepted_point: Option<IteratesVector>,
166    /// `no_bounds_` flag — port of `IpAdaptiveMuUpdate.cpp:282-287`.
167    /// Set to `true` on the first `update_barrier_parameter` call when
168    /// the iterate has zero bound multipliers (z_l, z_u, v_l, v_u all
169    /// have dim 0 — e.g. BT3, GENHS28, HS50, equality-only TNLPs).
170    /// Subsequent calls return `mu_min` immediately. Without this,
171    /// `mu_max = mu_max_fact * curr_avrg_compl()` evaluates to 0 (no
172    /// slacks → zero complementarity) and the later `clamp(mu_min,
173    /// mu_max)` panics with `min > max`.
174    no_bounds: bool,
175}
176
177impl Default for AdaptiveMuUpdate {
178    fn default() -> Self {
179        // Defaults from `IpAdaptiveMuUpdate.cpp:RegisterOptions`.
180        Self {
181            mu_oracle: MuOracleKind::QualityFunction,
182            adaptive_mu_globalization: AdaptiveMuGlobalization::ObjConstrFilter,
183            adaptive_mu_kkt_norm: AdaptiveMuKktNorm::TwoNormSquared,
184            adaptive_mu_safeguard_factor: 0.0,
185            adaptive_mu_kkterror_red_iters: 4,
186            adaptive_mu_kkterror_red_fact: 0.9999,
187            filter_max_margin: 1.0,
188            filter_margin_fact: 1e-5,
189            mu_min: 1e-11,
190            // Sentinel; lazy-initialised to `mu_max_fact * avrg_compl`
191            // on the first `update_barrier_parameter` call. Upstream
192            // `IpAdaptiveMuUpdate.cpp:164` sets `mu_max_ = -1.` when
193            // the option is not user-specified.
194            mu_max: -1.0,
195            mu_max_fact: 1e3,
196            tau_min: 0.99,
197            mu_init: 0.1,
198            barrier_tol_factor: 10.0,
199            mu_linear_decrease_factor: 0.2,
200            mu_superlinear_decrease_power: 1.5,
201            adaptive_mu_monotone_init_factor: 0.8,
202            restore_accepted_iterate: false,
203            sigma_max: 1e2,
204            sigma_min: 1e-6,
205            qf_norm_type: crate::mu::oracle::quality_function::NormType::TwoNormSquared,
206            qf_centrality_type: crate::mu::oracle::quality_function::CentralityType::None,
207            qf_balancing_term: crate::mu::oracle::quality_function::BalancingTermType::None,
208            qf_max_section_steps: 8,
209            qf_section_sigma_tol: 1e-2,
210            qf_section_qf_tol: 0.0,
211            probing_iterate_quality_factor: 1e4,
212            init_dual_inf: -1.0,
213            init_primal_inf: -1.0,
214            free_mu_mode: true,
215            refs_vals: VecDeque::new(),
216            filter: Filter::new(),
217            accepted_point: None,
218            no_bounds: false,
219        }
220    }
221}
222
223impl AdaptiveMuUpdate {
224    pub fn new() -> Self {
225        Self::default()
226    }
227
228    /// Pure-arithmetic predicate behind the probing-oracle iterate-
229    /// quality guard (pounce#58). Returns `true` when the ratio
230    /// `avrg_compl / curr_mu` exceeds `factor`. The two non-strict
231    /// gates (`factor > 0`, `curr_mu > 0`) keep the predicate
232    /// well-defined when the guard is disabled or when an unusual
233    /// μ-strategy zeroes `curr_mu`.
234    pub fn probing_iterate_guard_fires(
235        factor: Number,
236        curr_mu: Number,
237        avrg_compl: Number,
238    ) -> bool {
239        factor > 0.0 && curr_mu > 0.0 && avrg_compl > factor * curr_mu
240    }
241
242    /// Scalar core of `AdaptiveMuUpdate::lower_mu_safeguard`
243    /// (`IpAdaptiveMuUpdate.cpp:753-786`):
244    /// ```text
245    ///   init_dual_inf   ← max(1, dual_inf)   if not yet set
246    ///   init_primal_inf ← max(1, primal_inf) if not yet set
247    ///   lower = max(safeguard_factor * dual_inf / init_dual_inf,
248    ///               safeguard_factor * primal_inf / init_primal_inf)
249    ///   if globalization == KKT_ERROR: lower = min(lower, min_ref_val)
250    /// ```
251    pub fn lower_mu_safeguard(
252        &mut self,
253        dual_inf: Number,
254        primal_inf: Number,
255        min_ref_val: Number,
256    ) -> Number {
257        if self.init_dual_inf < 0.0 {
258            self.init_dual_inf = dual_inf.max(1.0);
259        }
260        if self.init_primal_inf < 0.0 {
261            self.init_primal_inf = primal_inf.max(1.0);
262        }
263        let dual_term = self.adaptive_mu_safeguard_factor * (dual_inf / self.init_dual_inf);
264        let prim_term = self.adaptive_mu_safeguard_factor * (primal_inf / self.init_primal_inf);
265        let mut lower = dual_term.max(prim_term);
266        if self.adaptive_mu_globalization == AdaptiveMuGlobalization::KktError {
267            lower = lower.min(min_ref_val);
268        }
269        lower
270    }
271
272    pub fn reset_init_inf(&mut self) {
273        self.init_dual_inf = -1.0;
274        self.init_primal_inf = -1.0;
275    }
276
277    /// Globalization KKT-error proxy — port of
278    /// `AdaptiveMuUpdate::quality_function_pd_system`
279    /// (`IpAdaptiveMuUpdate.cpp:629-744`). v1.0 hardwires the
280    /// max-norm variant (`adaptive_mu_kkt_norm_type=max-norm`,
281    /// upstream "NM_NORM_MAX") because the existing CQ surface
282    /// exposes max-norm primal/dual infeasibility cheaply; the
283    /// other three norm variants follow once `curr_*_infeasibility`
284    /// learns to dispatch on `NormEnum`. The score sums primal +
285    /// dual + complementarity (+ optional centrality / balancing
286    /// — both default off; left as `0`).
287    fn quality_function_pd_system(&self, cq: &IpoptCqHandle) -> Number {
288        let cq_ref = cq.borrow();
289        let primal_inf = cq_ref.curr_primal_infeasibility_max();
290        let dual_inf = cq_ref.curr_dual_infeasibility_max();
291        // Max-norm complementarity ≈ avrg_compl is a cheap proxy.
292        // Upstream's `curr_complementarity(0., NORM_MAX)` would use
293        // `||s ⊙ z||_∞`; absent that accessor, fall through to the
294        // average. For the monotonicity test inside
295        // `check_sufficient_progress` only ratios matter, so the
296        // proxy preserves the convergence criterion.
297        let complty = cq_ref.curr_avrg_compl();
298        primal_inf + dual_inf + complty
299    }
300
301    /// Port of `AdaptiveMuUpdate::CheckSufficientProgress`
302    /// (`IpAdaptiveMuUpdate.cpp:446-490`). Returns `true` if the
303    /// current iterate makes acceptable progress under the active
304    /// globalization rule.
305    fn check_sufficient_progress(&self, cq: &IpoptCqHandle) -> bool {
306        match self.adaptive_mu_globalization {
307            AdaptiveMuGlobalization::KktError => {
308                if self.refs_vals.len() < self.adaptive_mu_kkterror_red_iters.max(1) {
309                    // Not enough history yet — accept (matches
310                    // upstream's `num_refs >= num_refs_max_` guard).
311                    return true;
312                }
313                let curr_error = self.quality_function_pd_system(cq);
314                self.refs_vals
315                    .iter()
316                    .any(|&r| curr_error <= self.adaptive_mu_kkterror_red_fact * r)
317            }
318            AdaptiveMuGlobalization::ObjConstrFilter => {
319                let cq_ref = cq.borrow();
320                let curr_f = cq_ref.curr_f();
321                let curr_theta = cq_ref.curr_constraint_violation();
322                // `curr_nlp_error` is our analogue of upstream's
323                // global error margin driver.
324                let curr_err = cq_ref.curr_nlp_error();
325                drop(cq_ref);
326                let margin = self.filter_margin_fact * self.filter_max_margin.min(curr_err);
327                !self
328                    .filter
329                    .dominated_by_any(curr_theta + margin, curr_f + margin)
330            }
331            AdaptiveMuGlobalization::NeverMonotoneMode => true,
332        }
333    }
334
335    /// Port of `AdaptiveMuUpdate::RememberCurrentPointAsAccepted`
336    /// (`IpAdaptiveMuUpdate.cpp:492-546`). Records the iterate state
337    /// for the next sufficient-progress check.
338    fn remember_current_point_as_accepted(&mut self, data: &IpoptDataHandle, cq: &IpoptCqHandle) {
339        match self.adaptive_mu_globalization {
340            AdaptiveMuGlobalization::KktError => {
341                let curr_error = self.quality_function_pd_system(cq);
342                if self.refs_vals.len() >= self.adaptive_mu_kkterror_red_iters.max(1) {
343                    self.refs_vals.pop_front();
344                }
345                self.refs_vals.push_back(curr_error);
346            }
347            AdaptiveMuGlobalization::ObjConstrFilter => {
348                let cq_ref = cq.borrow();
349                let f = cq_ref.curr_f();
350                let theta = cq_ref.curr_constraint_violation();
351                let it = data.borrow().iter_count;
352                drop(cq_ref);
353                self.filter.add(theta, f, it);
354            }
355            AdaptiveMuGlobalization::NeverMonotoneMode => {}
356        }
357        if self.restore_accepted_iterate {
358            self.accepted_point = data.borrow().curr.clone();
359        }
360    }
361
362    /// Port of `AdaptiveMuUpdate::NewFixedMu`
363    /// (`IpAdaptiveMuUpdate.cpp:583-627`). Selects μ when the state
364    /// machine drops out of free mode. v1.0 always uses the
365    /// "average complementarity" branch (no `fix_mu_oracle_` is
366    /// wired; matches `fixed_mu_oracle = average_compl`).
367    fn new_fixed_mu(&self, cq: &IpoptCqHandle) -> Number {
368        let avrg = cq.borrow().curr_avrg_compl();
369        let new_mu = self.adaptive_mu_monotone_init_factor * avrg;
370        new_mu.clamp(self.mu_min, self.mu_max)
371    }
372}
373
374impl MuUpdate for AdaptiveMuUpdate {
375    /// Port of `IpAdaptiveMuUpdate.cpp:InitializeImpl`. Seeds
376    /// `curr_mu = mu_init`, `curr_tau = max(tau_min, 1 - mu_init)`,
377    /// resets the globalization state, and starts in free-μ mode
378    /// (`SetFreeMuMode(true)` at line 239).
379    fn initialize(&mut self, data: &IpoptDataHandle) {
380        // Mirror upstream `IpAdaptiveMuUpdate.cpp:246-247`:
381        //   IpData().Set_mu(1.);
382        //   IpData().Set_tau(0.);
383        // These are placeholder values so `CalculateSafeSlack` and the
384        // first output line have something to work with; the actual μ
385        // is computed by the oracle at iter 0's `update_barrier_parameter`.
386        // Setting curr_mu = mu_init here (as we used to) skipped the
387        // oracle's iter-0 call and locked μ at mu_init for the first
388        // Newton step — diverging from upstream's iter-0 behaviour
389        // (PFIT3: upstream iter 0 oracle picked μ=1.6e-6, pounce was
390        // stuck at μ=0.1, producing different iter-1 trial point).
391        let mut d = data.borrow_mut();
392        d.curr_mu = 1.0;
393        d.curr_tau = 0.0;
394        drop(d);
395        self.free_mu_mode = true;
396        self.refs_vals.clear();
397        self.filter.clear();
398        self.accepted_point = None;
399        self.init_dual_inf = -1.0;
400        self.init_primal_inf = -1.0;
401        // Reset mu_max sentinel so a re-solve re-runs the lazy init
402        // against the fresh starting iterate's curr_avrg_compl.
403        // Upstream re-enters InitializeImpl on each solve which
404        // (lines 160-165) resets `mu_max_ = -1.` when not user-set.
405        self.mu_max = -1.0;
406        // Reset no-bounds detection on re-solve.
407        self.no_bounds = false;
408    }
409
410    /// Adaptive μ update — port of `UpdateBarrierParameter`
411    /// (`IpAdaptiveMuUpdate.cpp:252-444`). Runs the FreeMuMode /
412    /// FixedMuMode state machine:
413    ///
414    /// * **FreeMuMode**: ask the configured oracle for a candidate
415    ///   (LOQO closed-form, Probing predictor solve, or
416    ///   QualityFunction golden-section). If progress is sufficient,
417    ///   stay in free mode and remember the iterate; otherwise switch
418    ///   to fixed mode at `new_fixed_mu`.
419    /// * **FixedMuMode**: monotone Fiacco-McCormick reduction
420    ///   (`min(linear · μ, μ^superlinear_power)`). Switch back to
421    ///   free mode once the globalization criterion is satisfied
422    ///   again.
423    ///
424    /// Probing / QualityFunction silently fall back to LOQO when
425    /// `nlp` / `pd_search_dir` are unavailable (mirrors upstream
426    /// lines 402-408).
427    ///
428    /// Note: line-search reset (upstream's `linesearch_->Reset()` at
429    /// lines 339, 386, 431) is not yet wired here — that handle is
430    /// not part of the [`MuUpdate`] trait surface. This is a
431    /// deliberate v1.0 deviation; it primarily affects the watchdog
432    /// counter, not convergence.
433    fn update_barrier_parameter(
434        &mut self,
435        data: &IpoptDataHandle,
436        cq: &IpoptCqHandle,
437        nlp: Option<&Rc<RefCell<dyn IpoptNlp>>>,
438        pd_search_dir: Option<&mut PdSearchDirCalc>,
439    ) -> Number {
440        // Lazy `mu_max` init — port of `IpAdaptiveMuUpdate.cpp:267-274`.
441        // Upstream computes `mu_max = mu_max_fact * curr_avrg_compl()`
442        // on the first call when the user did not set `mu_max`
443        // explicitly. Pounce previously hard-coded `mu_max = 1e5`,
444        // which let `new_fixed_mu = 0.8 * curr_avrg_compl` cap at 1e5
445        // — on DECONVBNE that allowed μ to jump from 2.5e-3 to ~2000
446        // at iter 198, destabilising the rest of the run.
447        if self.mu_max < 0.0 {
448            let avrg = cq.borrow().curr_avrg_compl();
449            self.mu_max = self.mu_max_fact * avrg;
450        }
451
452        // No-bounds short-circuit — port of `IpAdaptiveMuUpdate.cpp:282-296`.
453        // Detect once on the first call whether the iterate has any
454        // bound multipliers (z_l, z_u, v_l, v_u). When all four are
455        // dim-zero (equality-only TNLPs: BT3, GENHS28, HS50, METHANL8,
456        // ...), `curr_avrg_compl()` is 0, hence `mu_max = 0`, and the
457        // later `clamp(mu_min, mu_max)` panics with `min > max`.
458        // Upstream sets `mu = mu_min`, `tau = tau_min`, and short-
459        // circuits all subsequent oracle work; we mirror that.
460        if !self.no_bounds {
461            let n_bounds = {
462                let d = data.borrow();
463                let c = d.curr.as_ref().expect("curr set");
464                c.z_l.dim() + c.z_u.dim() + c.v_l.dim() + c.v_u.dim()
465            };
466            if n_bounds == 0 {
467                self.no_bounds = true;
468                let mut d = data.borrow_mut();
469                d.curr_mu = self.mu_min;
470                d.curr_tau = self.tau_min;
471                return self.mu_min;
472            }
473        }
474        if self.no_bounds {
475            let mut d = data.borrow_mut();
476            d.curr_mu = self.mu_min;
477            d.curr_tau = self.tau_min;
478            return self.mu_min;
479        }
480
481        // Read-and-clear `tiny_step_flag` — mirrors upstream
482        // `IpAdaptiveMuUpdate.cpp:297-298`. The flag is consumed by
483        // this call: without the clear, a single tiny-step detection
484        // would persist forever and suppress `sufficient_progress` on
485        // every later outer iter.
486        let (curr_mu, iter_count, tiny_step_flag) = {
487            let mut d = data.borrow_mut();
488            let out = (d.curr_mu, d.iter_count, d.tiny_step_flag);
489            d.tiny_step_flag = false;
490            out
491        };
492
493        // NB: do NOT short-circuit at iter_count==0. Upstream's
494        // `UpdateBarrierParameter` runs the oracle at iter 0 (the
495        // initialize() above set μ=1.0 as a placeholder only). Skipping
496        // the oracle here locked μ at the placeholder for the first
497        // Newton step. Letting the iter-0 path flow through the
498        // free-μ branch picks up the oracle's choice — the empty
499        // `refs_vals_` makes `check_sufficient_progress` return true,
500        // we remember the iterate, then call the oracle below.
501        // `tiny_step_flag` (and upstream's `CheckSkippedLineSearch()`,
502        // which is only set in non-rigorous resto mode) forces
503        // `sufficient_progress = false` when not in `NEVER_MONOTONE_MODE`
504        // — see `IpAdaptiveMuUpdate.cpp:347-351`. This is what lets a
505        // stalled outer iter drop into fixed-μ and re-seed μ via
506        // `new_fixed_mu` instead of the oracle re-driving μ further down.
507        let force_no_progress = tiny_step_flag
508            && self.adaptive_mu_globalization != AdaptiveMuGlobalization::NeverMonotoneMode;
509
510        if !self.free_mu_mode {
511            // Fixed-mu branch — `cpp:299-342`.
512            let sufficient_progress = !force_no_progress && self.check_sufficient_progress(cq);
513            if sufficient_progress {
514                // Switch back to free mode and record the iterate —
515                // upstream `cpp:303-311`. Upstream does NOT return
516                // here: after flipping `FreeMuMode` to true the first
517                // if/else ends and control reaches the `if
518                // FreeMuMode()` block at `cpp:391`, which runs the
519                // oracle and picks a fresh μ in the SAME iteration.
520                // Returning `curr_mu` here froze μ on the transition
521                // iter — PALMER4's iter-15 fixed→free transition kept
522                // μ at 2.4e-7 instead of letting the oracle drop it to
523                // mu_min, stalling to Maximum_Iterations_Exceeded.
524                // Fall through to the oracle call below.
525                self.free_mu_mode = true;
526                self.remember_current_point_as_accepted(data, cq);
527            } else {
528                // Keep reducing μ Fiacco-McCormick style if the
529                // barrier subproblem is solved to within
530                // `barrier_tol_factor · μ`, OR if a tiny step was
531                // just detected (`cpp:320` `|| tiny_step_flag`).
532                let sub_problem_error = cq.borrow().curr_barrier_error();
533                if sub_problem_error <= self.barrier_tol_factor * curr_mu || tiny_step_flag {
534                    let lin = self.mu_linear_decrease_factor * curr_mu;
535                    let sup = curr_mu.powf(self.mu_superlinear_decrease_power);
536                    let new_mu = lin.min(sup).max(self.mu_min).min(self.mu_max);
537                    let new_tau = self.tau_min.max(1.0 - new_mu);
538                    data.borrow_mut().curr_tau = new_tau;
539                    return new_mu;
540                }
541                // Subproblem not yet solved — keep μ.
542                let new_tau = self.tau_min.max(1.0 - curr_mu);
543                data.borrow_mut().curr_tau = new_tau;
544                return curr_mu;
545            }
546        } else {
547            // Free-mu branch — `cpp:343-389`.
548            let sufficient_progress = !force_no_progress && self.check_sufficient_progress(cq);
549            if sufficient_progress {
550                self.remember_current_point_as_accepted(data, cq);
551                // Fall through to the oracle call below.
552            } else {
553                if std::env::var("POUNCE_DBG_AMU").is_ok() {
554                    let cqr = cq.borrow();
555                    let theta = cqr.curr_constraint_violation();
556                    let f = cqr.curr_f();
557                    let nlp_err = cqr.curr_nlp_error();
558                    let avrg = cqr.curr_avrg_compl();
559                    drop(cqr);
560                    let margin = self.filter_margin_fact * self.filter_max_margin.min(nlp_err);
561                    let entries: Vec<(Number, Number, i32)> = self
562                        .filter
563                        .entries()
564                        .iter()
565                        .map(|e| (e.theta, e.phi, e.iter))
566                        .collect();
567                    tracing::debug!(target: "pounce::mu",
568                        "[AMU] iter={} free->fixed: curr_mu={:.3e} theta={:.3e} f={:.3e} nlp_err={:.3e} margin={:.3e} avrg_compl={:.3e} new_mu={:.3e} | filter={:?} | force_no_progress={} tiny={}",
569                        iter_count,
570                        curr_mu,
571                        theta,
572                        f,
573                        nlp_err,
574                        margin,
575                        avrg,
576                        self.adaptive_mu_monotone_init_factor * avrg,
577                        entries,
578                        force_no_progress,
579                        tiny_step_flag,
580                    );
581                }
582                // Switch into fixed mode.
583                self.free_mu_mode = false;
584                if self.restore_accepted_iterate {
585                    if let Some(prev) = self.accepted_point.clone() {
586                        let mut d = data.borrow_mut();
587                        d.set_trial(prev);
588                        d.accept_trial_point();
589                    }
590                }
591                let new_mu = self.new_fixed_mu(cq);
592                let new_tau = self.tau_min.max(1.0 - new_mu);
593                data.borrow_mut().curr_tau = new_tau;
594                return new_mu;
595            }
596        }
597
598        // ----- Free-mu oracle call (cpp:391-436) -----
599        let cq_ref = cq.borrow();
600        let dual_inf = cq_ref.curr_dual_infeasibility_max();
601        let primal_inf = cq_ref.curr_primal_infeasibility_max();
602        let avrg_compl = cq_ref.curr_avrg_compl();
603        let centrality_xi = cq_ref.curr_centrality_measure();
604        let nlp_error = cq_ref.curr_nlp_error();
605        drop(cq_ref);
606
607        // τ = max(tau_min, 1 - curr_nlp_error) — upstream cpp:397.
608        let tau = self.tau_min.max(1.0 - nlp_error);
609        data.borrow_mut().curr_tau = tau;
610
611        let loqo_candidate = || {
612            let mut oracle = LoqoMuOracle {
613                mu_min: self.mu_min,
614                mu_max: self.mu_max,
615                avrg_compl,
616                centrality_xi,
617            };
618            oracle.calculate_mu().unwrap_or(curr_mu)
619        };
620
621        let candidate = match self.mu_oracle {
622            MuOracleKind::Loqo => loqo_candidate(),
623            MuOracleKind::Probing => {
624                // Iterate-quality guard (pounce#58). The probing
625                // oracle uses `curr_avrg_compl()` for its `mu_curr`
626                // input (see `mu/oracle/probing.rs:85`). When a single
627                // imbalanced `(s_i, z_i)` pair inflates the average
628                // many orders above the stored `data.curr_mu`,
629                // probing's `σ·mu_curr` correctly returns the inflated
630                // value and the resulting search direction throws the
631                // iterate out of the convergence neighborhood. On
632                // arki0012 this manifests as μ jumping 5 orders at
633                // iter 155 followed by divergence to "Local
634                // Infeasibility" at iter 284. We short-circuit by
635                // signalling restoration and keeping μ unchanged; the
636                // main loop in `ipopt_alg.rs` consumes the flag
637                // before the search-direction step.
638                if Self::probing_iterate_guard_fires(
639                    self.probing_iterate_quality_factor,
640                    curr_mu,
641                    avrg_compl,
642                ) {
643                    if std::env::var("POUNCE_DBG_ORACLE").is_ok() {
644                        tracing::debug!(target: "pounce::mu",
645                            "[PN_PROBE_GUARD] iter={} curr_mu={:.3e} avrg_compl={:.3e} ratio={:.3e} > factor={:.3e} → request_resto",
646                            iter_count,
647                            curr_mu,
648                            avrg_compl,
649                            avrg_compl / curr_mu,
650                            self.probing_iterate_quality_factor,
651                        );
652                    }
653                    data.borrow_mut().request_resto = true;
654                    return curr_mu;
655                }
656                match (nlp, pd_search_dir) {
657                    (Some(nlp), Some(sd)) => {
658                        let mut oracle = ProbingMuOracle {
659                            sigma_max: 100.0,
660                            mu_min: self.mu_min,
661                            mu_max: self.mu_max,
662                            mu_curr: curr_mu,
663                            mu_aff: curr_mu,
664                        };
665                        oracle
666                            .calculate_mu_with_affine_step(data, cq, nlp, sd, 1.0)
667                            .unwrap_or_else(loqo_candidate)
668                    }
669                    _ => loqo_candidate(),
670                }
671            }
672            MuOracleKind::QualityFunction => match (nlp, pd_search_dir) {
673                (Some(nlp), Some(sd)) => {
674                    let mut oracle = QualityFunctionMuOracle::new();
675                    oracle.mu_min = self.mu_min;
676                    oracle.mu_max = self.mu_max;
677                    oracle.sigma_min = self.sigma_min;
678                    oracle.sigma_max = self.sigma_max;
679                    oracle.norm_type = self.qf_norm_type;
680                    oracle.centrality_type = self.qf_centrality_type;
681                    oracle.balancing_term = self.qf_balancing_term;
682                    oracle.max_section_steps = self.qf_max_section_steps;
683                    oracle.section_sigma_tol = self.qf_section_sigma_tol;
684                    oracle.section_qf_tol = self.qf_section_qf_tol;
685                    // Mirrors upstream's `quality_function_search` timer
686                    // around `CalculateMu` in `IpQualityFunctionMuOracle.cpp`.
687                    let timing = data.borrow().timing.clone();
688                    let _qf_guard = timing.quality_function_search.guard();
689                    oracle
690                        .calculate_mu_with_predictor_centering(data, cq, nlp, sd)
691                        .unwrap_or_else(loqo_candidate)
692                }
693                _ => loqo_candidate(),
694            },
695        };
696
697        // Safeguard floor + global band clamp (cpp:410-426).
698        let lower = self.lower_mu_safeguard(dual_inf, primal_inf, candidate);
699        let mu = candidate.max(self.mu_min).max(lower).min(self.mu_max);
700
701        // NB: upstream `IpAdaptiveMuUpdate.cpp:410-426` does NOT require
702        // `mu ≤ curr_mu` in free mode — the oracle is allowed to bump
703        // μ back up. A prior attempt to cap growth here ("HAIFAM
704        // stability hack") let DECONVBNE's μ plunge from 0.1 to 5e-10
705        // in ~20 iters and never recover (upstream oscillates μ in
706        // [-8,-1] for the same range), trapping `inf_du` at 1e13.
707        // Tiny-step skips are already handled by the
708        // `tiny_step_flag → force_no_progress → new_fixed_mu` path
709        // above, which can raise μ via the fixed-mode branch.
710        mu
711    }
712}
713
714#[cfg(test)]
715mod tests {
716    use super::*;
717
718    #[test]
719    fn lower_mu_safeguard_initializes_from_first_call() {
720        let mut a = AdaptiveMuUpdate::new();
721        a.adaptive_mu_safeguard_factor = 1e-2;
722        // First call captures init values.
723        let _ = a.lower_mu_safeguard(0.5, 2.0, 1.0);
724        assert_eq!(a.init_dual_inf, 1.0); // max(1, 0.5)
725        assert_eq!(a.init_primal_inf, 2.0); // max(1, 2.0)
726    }
727
728    #[test]
729    fn lower_mu_safeguard_takes_max_of_dual_and_primal_terms() {
730        let mut a = AdaptiveMuUpdate::new();
731        a.adaptive_mu_safeguard_factor = 1.0;
732        // Primal term dominates.
733        let r = a.lower_mu_safeguard(0.1, 5.0, 1e9);
734        // init_dual = 1, init_primal = 5 → terms: 0.1, 1.0 → max = 1.0.
735        assert!((r - 1.0).abs() < 1e-15);
736    }
737
738    #[test]
739    fn kkt_error_globalization_clips_to_min_ref_val() {
740        let mut a = AdaptiveMuUpdate::new();
741        a.adaptive_mu_globalization = AdaptiveMuGlobalization::KktError;
742        a.adaptive_mu_safeguard_factor = 1.0;
743        // Without clip, safeguard would be 5.0; min_ref_val = 0.1 wins.
744        let r = a.lower_mu_safeguard(0.1, 5.0, 0.1);
745        assert!((r - 0.1).abs() < 1e-15);
746    }
747
748    #[test]
749    fn reset_clears_init_inf() {
750        let mut a = AdaptiveMuUpdate::new();
751        a.adaptive_mu_safeguard_factor = 1.0;
752        let _ = a.lower_mu_safeguard(0.5, 2.0, 1.0);
753        a.reset_init_inf();
754        assert_eq!(a.init_dual_inf, -1.0);
755        assert_eq!(a.init_primal_inf, -1.0);
756    }
757
758    // The trait `update_barrier_parameter` now takes
759    // `(&IpoptDataHandle, &IpoptCqHandle)`. End-to-end coverage of the
760    // adaptive path lands alongside the integration test that drives
761    // `IpoptAlgorithm::optimize` with `mu_strategy=adaptive`; in
762    // isolation the unit tests above exercise the safeguard
763    // arithmetic and option defaults.
764
765    #[test]
766    fn default_mu_oracle_is_quality_function() {
767        let a = AdaptiveMuUpdate::new();
768        assert_eq!(a.mu_oracle, MuOracleKind::QualityFunction);
769    }
770
771    #[test]
772    fn mu_oracle_kind_is_distinct() {
773        assert_ne!(MuOracleKind::Loqo, MuOracleKind::Probing);
774        assert_ne!(MuOracleKind::Probing, MuOracleKind::QualityFunction);
775        assert_ne!(MuOracleKind::Loqo, MuOracleKind::QualityFunction);
776    }
777
778    // pounce#58 guard predicate. Numbers below come from the issue
779    // body's iter 154-155 trace on arki0012.
780    #[test]
781    fn probing_iterate_guard_fires_on_arki0012_iter155() {
782        let curr_mu = 1.98e-11;
783        let avrg_compl = 8.90e-6;
784        assert!(AdaptiveMuUpdate::probing_iterate_guard_fires(
785            1e4, curr_mu, avrg_compl
786        ));
787    }
788
789    #[test]
790    fn probing_iterate_guard_quiet_on_healthy_iter() {
791        // iter 154 in the same trace — ratio ≈ 2.2; ought not fire.
792        let curr_mu = 1.02e-11;
793        let avrg_compl = 2.24e-11;
794        assert!(!AdaptiveMuUpdate::probing_iterate_guard_fires(
795            1e4, curr_mu, avrg_compl
796        ));
797    }
798
799    #[test]
800    fn probing_iterate_guard_disabled_at_zero_factor() {
801        // factor=0 ⇒ guard off, even with extreme ratio.
802        assert!(!AdaptiveMuUpdate::probing_iterate_guard_fires(
803            0.0, 1e-11, 1.0
804        ));
805    }
806
807    #[test]
808    fn probing_iterate_guard_disabled_at_negative_factor() {
809        assert!(!AdaptiveMuUpdate::probing_iterate_guard_fires(
810            -1.0, 1e-11, 1.0
811        ));
812    }
813
814    #[test]
815    fn probing_iterate_guard_quiet_when_curr_mu_zero() {
816        // Pathological `curr_mu = 0` (no-bounds branch zeroes it out).
817        // Predicate must stay quiet rather than division-by-zero.
818        assert!(!AdaptiveMuUpdate::probing_iterate_guard_fires(
819            1e4, 0.0, 1e-6
820        ));
821    }
822
823    #[test]
824    fn probing_iterate_guard_threshold_at_factor_times_mu() {
825        // Boundary: equality does NOT fire (strict >).
826        let curr_mu = 1.0e-10;
827        let factor = 1e4;
828        assert!(!AdaptiveMuUpdate::probing_iterate_guard_fires(
829            factor,
830            curr_mu,
831            factor * curr_mu
832        ));
833        // Just above the boundary fires.
834        assert!(AdaptiveMuUpdate::probing_iterate_guard_fires(
835            factor,
836            curr_mu,
837            factor * curr_mu * (1.0 + 1e-12)
838        ));
839    }
840}