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}