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