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}