pounce_algorithm/ipopt_alg.rs
1//! Main optimization loop — port of
2//! `Algorithm/IpIpoptAlg.{hpp,cpp}`.
3//!
4//! Phase 7 ships the loop scaffold matching `Optimize()` lines
5//! 292-563 in upstream. The body invokes:
6//!
7//! 1. `IterateInitializer::set_initial_iterates`
8//! 2. (loop) `OutputIteration` → `CheckConvergence` →
9//! `UpdateBarrierParameter` → `UpdateHessian` →
10//! `ComputeSearchDirection` → `ComputeAcceptableTrialPoint` →
11//! `AcceptTrialPoint`
12//! 3. `correct_bound_multiplier` (kappa_sigma) per `MAIN_LOOP.md`
13//! §"Bound multiplier reset" lines 1055-1134
14//! 4. exception → `SolverReturn` mapping per the table in
15//! `MAIN_LOOP.md`.
16//!
17//! The NLP handle and search-direction calculator are optional:
18//! when both are present, `iterate()` computes a real Newton step and
19//! drives the line search. Without them, `iterate()` runs the bookkeeping
20//! pieces (mu update, hessian update, conv check, kappa_sigma reset)
21//! and is exercised by structural unit tests. The full path lights up
22//! once `pounce-nlp::OrigIpoptNLP` lands.
23
24use crate::alg_builder::AlgorithmBundle;
25use crate::conv_check::r#trait::ConvergenceStatus;
26use crate::intermediate::{CtxGuard, IntermediateContext};
27use crate::ipopt_cq::IpoptCqHandle;
28use crate::ipopt_data::IpoptDataHandle;
29use crate::ipopt_nlp::IpoptNlp;
30use crate::iter_dump::IterDumper;
31use crate::iterate_dump::emit_record as emit_iterate_record;
32use crate::kkt::pd_search_dir_calc::PdSearchDirCalc;
33use crate::line_search::backtracking::Outcome;
34use crate::restoration::{RestorationOutcome, RestorationPhase};
35use pounce_common::diagnostics::DiagnosticsState;
36use pounce_common::types::{Index, Number};
37use pounce_linalg::Vector;
38use pounce_nlp::alg_types::SolverReturn;
39use pounce_nlp::return_codes::AlgorithmMode;
40use pounce_nlp::tnlp::{IpoptCq as TnlpIpoptCq, IpoptData as TnlpIpoptData, IterStats, TNLP};
41use std::cell::RefCell;
42use std::rc::Rc;
43
44pub struct IpoptAlgorithm {
45 pub data: IpoptDataHandle,
46 pub cq: IpoptCqHandle,
47 pub bundle: AlgorithmBundle,
48 /// Optional NLP handle. Required for any step that evaluates
49 /// problem functions or pulls bound expansion matrices (init,
50 /// search direction, line-search trial-point evaluation). Absent
51 /// in the structural unit tests of Phases 5-6.
52 pub nlp: Option<Rc<RefCell<dyn IpoptNlp>>>,
53 /// Optional TNLP handle — the user-facing problem. When present,
54 /// `iterate()` fires `TNLP::intermediate_callback` once per outer
55 /// iteration so callers can monitor progress or request early
56 /// termination (returning `false` from the callback surfaces as
57 /// `SolverReturn::UserRequestedStop`). Kept separate from `nlp`
58 /// because the algorithm-side NLP is the *compressed* `OrigIpoptNlp`
59 /// view (fixed-variable elimination, c/d split) while the callback
60 /// payload needs to expose the original-coordinate iterate.
61 pub tnlp: Option<Rc<RefCell<dyn TNLP>>>,
62 /// Search-direction calculator (`PdSearchDirCalc`). Lands once a
63 /// concrete `SymLinearSolver` backend (MUMPS / FERAL) is wired
64 /// through `AlgBuilder` in Phase 7's tail.
65 pub search_dir: Option<PdSearchDirCalc>,
66 /// Restoration-phase strategy. Invoked when the line search
67 /// returns [`Outcome::Failed`] (port of upstream
68 /// `IpBacktrackingLineSearch::ActivateLineSearch`'s resto
69 /// fallback). Optional: in its absence, line-search failure maps
70 /// directly to [`SolverReturn::RestorationFailure`] so the main
71 /// loop's exit-code semantics match upstream's "no resto built"
72 /// case.
73 pub restoration: Option<Box<dyn RestorationPhase>>,
74
75 /// `kappa_sigma` for the post-AcceptTrialPoint multiplier reset
76 /// (`IpIpoptAlg.cpp:correct_bound_multiplier`, line 1055-1134).
77 pub kappa_sigma: Number,
78 pub max_iter: Index,
79 /// Initial primal step length offered to the line search at the
80 /// top of each iteration. Mirrors `IpBacktrackingLineSearch`'s
81 /// fraction-to-the-boundary primal step (with τ = `data.curr_tau`).
82 /// In v1.0 the structural value here is 1.0 and the FTB cap is
83 /// applied per-component when the line-search driver computes
84 /// trial slacks; the simplification holds for non-degenerate runs.
85 pub alpha_init: Number,
86 /// Tiny-step relative tolerance — port of upstream
87 /// `IpBacktrackingLineSearch::tiny_step_tol_` (default `10·EPSILON`).
88 /// Step is "tiny" when `max_i |δx_i|/(1+|x_i|) ≤ tiny_step_tol`
89 /// (and same for s, and `c_viol ≤ 1e-4`).
90 pub tiny_step_tol: Number,
91 /// Port of upstream `IpIpoptAlg.cpp` divergence guard: when
92 /// `max_i |x_i|` exceeds this threshold the optimization aborts with
93 /// `SolverReturn::DivergingIterates`. Default `1e20` matches the
94 /// registered `diverging_iterates_tol` option. Catches MESH and
95 /// similar cases where the normal-mode IPM heads off to infinity
96 /// (orig `f` to ±1e33 by iter 90) before line-search failure forces
97 /// a degenerate restoration entry.
98 pub diverging_iterates_tol: Number,
99 /// Companion threshold on the dual step — when both primal and dual
100 /// steps are tiny in two consecutive iterations the algorithm
101 /// declares convergence at the best attainable accuracy. Default
102 /// `1e-2` matches upstream.
103 pub tiny_step_y_tol: Number,
104 /// Set true when the previous iterate was tagged tiny; on the
105 /// second consecutive tiny step the loop sets `data.tiny_step_flag`
106 /// so the mu update can attempt to terminate. Mirrors
107 /// `IpBacktrackingLineSearch::tiny_step_last_iteration_`.
108 pub tiny_step_last_iteration: bool,
109 /// Cycle-detection state for [`Self::invoke_restoration`]: the
110 /// outer `(x, s)` snapshot from the previous restoration entry,
111 /// cleared on any iteration that exits via a normal line-search
112 /// accept. When restoration is invoked twice in a row and the
113 /// outer iterate has not moved between entries (relative
114 /// 2-norm < 1e-10 on both `x` and `s`), the inner resto-IPM is
115 /// returning Recovered points indistinguishable from `curr` — a
116 /// cycle. Surfaces as `ErrorInStepComputation`. Mirrors the
117 /// *intent* of upstream `IpBacktrackingLineSearch.cpp:580-600`'s
118 /// almost-feasible-resto guard while staying robust against the
119 /// `inf_pr` micro-drift seen on ACOPR14 (delta ~3e-12 per entry,
120 /// inf_du essentially constant) where a scalar-`inf_pr` heuristic
121 /// fails. Productive single-restoration sequences (BT8, HIMMELBJ,
122 /// LINSPANH, LSNNODOC, ODFITS, OET3) clear the snapshot via
123 /// `Outcome::Accepted` between entries and are unaffected.
124 last_resto_entry_x: Option<Box<dyn Vector>>,
125 last_resto_entry_s: Option<Box<dyn Vector>>,
126 /// Snapshot of the *recovery* iterate from the previous
127 /// restoration. Compared against the next entry's `(x, s)` to
128 /// detect "outer made no progress between consecutive resto
129 /// invocations". When this distance is below threshold for
130 /// several consecutive entries, terminate — catching
131 /// slow-non-convergence cycles (ACOPR14, TRO3X3, ACOPR30) where
132 /// resto's *inner* moves substantively each call but the *outer*
133 /// makes no progress between calls. Cleared on any LS-accepted
134 /// step.
135 last_resto_recovery_x: Option<Box<dyn Vector>>,
136 last_resto_recovery_s: Option<Box<dyn Vector>>,
137 /// Count of consecutive restoration entries on which the outer
138 /// step (recovery → next-entry) was below the iterate-distance
139 /// threshold. Cleared on any LS-accepted step. Limit chosen to
140 /// let MAKELA3, HAIFAM, HALDMADS, ROBOT, TENBARS2 — which need
141 /// 2-3 consecutive resto entries to recover — pass through.
142 resto_no_outer_progress_count: usize,
143 /// Count of consecutive restoration entries on which the outer
144 /// constraint violation at entry was already below `tol` (the
145 /// outer optimality tolerance). Matches the *intent* of upstream
146 /// `IpBacktrackingLineSearch.cpp:580-600`'s almost-feasible-resto
147 /// guard while using a looser cv threshold (`tol` vs `1e-2·tol`)
148 /// — catches DECONVBNE's resto-thrash where each cycle re-enters
149 /// at cv ≈ 3e-10 < tol with bound multipliers reset to 1, the
150 /// outer's σ-blowup explodes inf_du to 1.9e7, alpha-min triggers
151 /// resto re-entry, and the (inf_pr, inf_du) post-recovery state
152 /// is essentially identical across cycles but `x` drifts enough
153 /// that [`Self::last_resto_recovery_x`]-based detection misses.
154 /// Cumulative (never cleared on LS-accept), since DECONVBNE's
155 /// cycle interleaves R-recoveries with sub-tol accepts that
156 /// accomplish no real outer progress. Fires after 3 near-feasible
157 /// entries — surfaces as `StopAtAcceptablePoint` since the
158 /// recovered point already satisfies constraint feasibility
159 /// within `tol`.
160 resto_near_feasible_count: usize,
161 /// Snapshot of the most recent iterate that the convergence check
162 /// flagged "acceptable" (NLP error ≤ `acceptable_tol`). Mirrors
163 /// upstream `IpBacktrackingLineSearch::acceptable_iterate_`
164 /// (`IpBacktrackingLineSearch.cpp:1286-1310`). Used by
165 /// [`Self::restore_acceptable_point`] to roll back when restoration
166 /// fails — if such an iterate exists, the algorithm exits with
167 /// `SolverReturn::StopAtAcceptablePoint` rather than
168 /// `RestorationFailure`. Cleared/refreshed on every iteration that
169 /// satisfies the acceptable predicate.
170 acceptable_iterate: Option<crate::iterates_vector::IteratesVector>,
171 acceptable_iter_number: Index,
172 /// Shared per-solve diagnostics state. `None` unless the CLI
173 /// requested `--dump <cat>:<spec>`. When set, the outer loop
174 /// advances the state's iter counter and the augmented-system
175 /// solver consults it to gate KKT dumps.
176 diagnostics: Option<Rc<DiagnosticsState>>,
177 /// Optional interactive debugger. Shared (`Rc<RefCell<…>>`) so the
178 /// same debugger instance also drives the restoration inner IPM —
179 /// one debugger sees both levels. Fired at every
180 /// [`crate::debug::Checkpoint`]. See `crate::debug`.
181 debug: Option<Rc<RefCell<dyn crate::debug::DebugHook>>>,
182
183 // ---- Restoration-phase audit counters (pounce#12). ----
184 //
185 // Drained into `SolveStatistics` by `IpoptApplication::optimize_constrained`
186 // after the solve completes. Counts are cumulative across the run.
187 /// Number of `invoke_restoration` entries.
188 pub resto_calls: Index,
189 /// Sum of inner-IPM iter counts across every restoration call.
190 pub resto_inner_iters: Index,
191 /// Number of outer iters that ran in restoration mode (R-line
192 /// equivalents in `print_level=5` output).
193 pub resto_outer_iters: Index,
194 /// Cumulative wall-clock seconds spent inside `perform_restoration`.
195 pub resto_wall_secs: Number,
196
197 // ---- Per-iteration history capture (pounce#8, pounce#71). ----
198 //
199 // The per-iteration trajectory is no longer accumulated on the
200 // algorithm: `iterate()` emits a structured `pounce::iteration`
201 // event each step, and `pounce_observability::IterCollectorLayer`
202 // rebuilds the `IterRecord`s into the active `IterCaptureGuard`
203 // that `IpoptApplication` installs around the solve.
204 /// When `false`, the per-iteration table that `iterate()` writes
205 /// straight to stdout is suppressed. Wired from
206 /// `IpoptApplication`'s `print_level` option: level 0 turns this
207 /// off (matches upstream's "no console output" contract). Default
208 /// `true` so CLI / direct-driver users keep the familiar trace.
209 pub print_iter_output: bool,
210}
211
212impl IpoptAlgorithm {
213 pub fn new(data: IpoptDataHandle, cq: IpoptCqHandle, mut bundle: AlgorithmBundle) -> Self {
214 // The builder may pre-populate `bundle.search_dir` when given a
215 // `LinearBackendFactory`; lift it onto the algorithm so the
216 // iterate body can call into it directly.
217 let search_dir = bundle.search_dir.take();
218 Self {
219 data,
220 cq,
221 bundle,
222 nlp: None,
223 tnlp: None,
224 search_dir,
225 restoration: None,
226 kappa_sigma: 1e10,
227 max_iter: 3000,
228 alpha_init: 1.0,
229 tiny_step_tol: 10.0 * Number::EPSILON,
230 diverging_iterates_tol: 1e20,
231 tiny_step_y_tol: 1e-2,
232 tiny_step_last_iteration: false,
233 last_resto_entry_x: None,
234 last_resto_entry_s: None,
235 last_resto_recovery_x: None,
236 last_resto_recovery_s: None,
237 resto_no_outer_progress_count: 0,
238 resto_near_feasible_count: 0,
239 acceptable_iterate: None,
240 acceptable_iter_number: 0,
241 diagnostics: None,
242 debug: None,
243 resto_calls: 0,
244 resto_inner_iters: 0,
245 resto_outer_iters: 0,
246 resto_wall_secs: 0.0,
247 print_iter_output: true,
248 }
249 }
250
251 /// Stash the current iterate as the "last acceptable" backup —
252 /// port of `IpBacktrackingLineSearch::StoreAcceptablePoint`
253 /// (`IpBacktrackingLineSearch.cpp:1286-1293`).
254 fn store_acceptable_point(&mut self) {
255 let d = self.data.borrow();
256 if let Some(curr) = d.curr.as_ref() {
257 self.acceptable_iterate = Some(curr.clone());
258 self.acceptable_iter_number = d.iter_count;
259 }
260 }
261
262 /// Roll the iterate back to the last acceptable snapshot — port of
263 /// `IpBacktrackingLineSearch::RestoreAcceptablePoint`
264 /// (`IpBacktrackingLineSearch.cpp:1295-1310`). Returns `true` if a
265 /// snapshot was available and applied; `false` otherwise (caller
266 /// then surfaces the original failure status).
267 fn restore_acceptable_point(&mut self) -> bool {
268 let Some(prev) = self.acceptable_iterate.clone() else {
269 return false;
270 };
271 let mut d = self.data.borrow_mut();
272 d.set_trial(prev);
273 // `accept_trial_point` promotes `trial → curr`, mirroring the
274 // upstream sequence `set_trial(...); AcceptTrialPoint();`.
275 d.accept_trial_point();
276 true
277 }
278
279 /// Terminal fallback for a near-feasible numerical breakdown (a
280 /// restoration cycle or a failed step computation). If a finite
281 /// acceptable iterate was recorded earlier in the solve, roll back
282 /// to it and stop at [`SolverReturn::StopAtAcceptablePoint`] (mapped
283 /// by the application layer to `Solved_To_Acceptable_Level`) rather
284 /// than surfacing the hard `fallback` error. This mirrors upstream
285 /// `IpBacktrackingLineSearch`'s `ACCEPTABLE_POINT_REACHED`
286 /// precedence: when the line search exhausts but an acceptable point
287 /// was stored, that point is returned instead of the failure. With
288 /// no snapshot — or if the restored objective is non-finite — the
289 /// original `fallback` status is surfaced unchanged, so genuinely
290 /// failed/infeasible solves keep their honest status. Catches
291 /// degenerate LPs (kleemin8, nsir2) whose μ-endgame reaches the
292 /// optimum, then destabilizes on the ill-conditioned vertex and
293 /// cycles in restoration instead of stopping at the acceptable
294 /// iterate it already passed through.
295 fn terminate_acceptable_or(&mut self, fallback: SolverReturn) -> IterateOutcome {
296 if self.restore_acceptable_point() && self.cq.borrow().curr_f().is_finite() {
297 IterateOutcome::Terminate(SolverReturn::StopAtAcceptablePoint)
298 } else {
299 IterateOutcome::Terminate(fallback)
300 }
301 }
302
303 pub fn with_nlp(mut self, nlp: Rc<RefCell<dyn IpoptNlp>>) -> Self {
304 self.nlp = Some(nlp);
305 self
306 }
307
308 /// Install a user-facing TNLP handle. Enables per-iteration
309 /// `TNLP::intermediate_callback` invocation from `optimize()`.
310 pub fn with_tnlp(mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> Self {
311 self.tnlp = Some(tnlp);
312 self
313 }
314
315 /// Build an [`IterStats`] payload from the current `IpoptData` /
316 /// `IpoptCq` state. Mirrors the field set the upstream Ipopt main
317 /// loop hands to `IntermediateCallback` after each `AcceptTrialPoint`.
318 fn build_iter_stats(&self) -> IterStats {
319 let d = self.data.borrow();
320 let c = self.cq.borrow();
321 let dnrm = match d.delta.as_ref() {
322 Some(delta) => delta.x.amax().max(delta.s.amax()),
323 None => 0.0,
324 };
325 IterStats {
326 // alg_mod tracking (regular vs restoration) is a follow-up;
327 // every callback fire from the outer loop reports RegularMode.
328 mode: AlgorithmMode::RegularMode,
329 iter: d.iter_count,
330 obj_value: c.curr_f(),
331 inf_pr: c.curr_primal_infeasibility_max(),
332 inf_du: c.curr_dual_infeasibility_max(),
333 mu: d.curr_mu,
334 d_norm: dnrm,
335 regularization_size: d.info_regu_x,
336 alpha_du: d.info_alpha_dual,
337 alpha_pr: d.info_alpha_primal,
338 ls_trials: d.info_ls_count,
339 }
340 }
341
342 /// Fire `TNLP::intermediate_callback` if a TNLP handle and NLP
343 /// handle are installed. Wraps the call in an [`IntermediateContext`]
344 /// guard so downstream inspector entry points (the C API's
345 /// `GetIpoptCurrent*`) can read live state for the duration. Returns
346 /// `true` to continue, `false` if the user requested termination.
347 fn fire_intermediate(&self) -> bool {
348 let Some(tnlp) = self.tnlp.as_ref() else {
349 return true;
350 };
351 let Some(nlp) = self.nlp.as_ref() else {
352 return true;
353 };
354 let stats = self.build_iter_stats();
355 let _guard = CtxGuard::install(IntermediateContext {
356 data: Rc::clone(&self.data),
357 cq: Rc::clone(&self.cq),
358 nlp: Rc::clone(nlp),
359 });
360 tnlp.borrow_mut().intermediate_callback(
361 stats,
362 &TnlpIpoptData::default(),
363 &TnlpIpoptCq::default(),
364 )
365 }
366
367 pub fn with_search_dir(mut self, sd: PdSearchDirCalc) -> Self {
368 self.search_dir = Some(sd);
369 self
370 }
371
372 pub fn with_restoration(mut self, resto: Box<dyn RestorationPhase>) -> Self {
373 self.restoration = Some(resto);
374 self
375 }
376
377 /// Install the shared diagnostics state. The state is propagated
378 /// to the augmented-system solver at the top of [`Self::optimize`]
379 /// so dump sites can consult per-iter gating.
380 pub fn with_diagnostics(mut self, diag: Rc<DiagnosticsState>) -> Self {
381 self.diagnostics = Some(diag);
382 self
383 }
384
385 /// Install an interactive debugger hook. Fired at each checkpoint
386 /// in [`Self::optimize`]; returning [`crate::debug::DebugAction::Stop`]
387 /// ends the solve with `SolverReturn::UserRequestedStop`.
388 pub fn with_debug_hook(mut self, hook: Rc<RefCell<dyn crate::debug::DebugHook>>) -> Self {
389 self.debug = Some(hook);
390 self
391 }
392
393 /// Shared handle to the installed debugger, if any — used to forward
394 /// it into the restoration inner IPM.
395 pub fn debug_hook(&self) -> Option<Rc<RefCell<dyn crate::debug::DebugHook>>> {
396 self.debug.as_ref().map(Rc::clone)
397 }
398
399 /// Fire the debugger hook (if installed) at `cp`, building a live
400 /// [`crate::debug::DebugCtx`] over cheap handle clones. Returns the
401 /// requested action, defaulting to `Resume` when no hook is set.
402 fn fire_debug(&mut self, cp: crate::debug::Checkpoint) -> crate::debug::DebugAction {
403 use crate::debug::{DebugAction, DebugCtx};
404 // Clone the Rc so the hook borrow is released before we touch
405 // `self.bundle` to apply any live option changes below.
406 let Some(hook) = self.debug.as_ref().map(Rc::clone) else {
407 return DebugAction::Resume;
408 };
409 let mut ctx = DebugCtx::new(Rc::clone(&self.data), Rc::clone(&self.cq), cp);
410 let action = hook.borrow_mut().at_checkpoint(&mut ctx);
411 // Drain any tolerances the hook asked to hot-swap and write them
412 // into the live convergence-check policy, so the next iteration's
413 // termination test uses the new value (no `resolve` needed).
414 for (name, value) in ctx.take_live_tolerances() {
415 self.bundle.conv_check.set_tolerance(&name, value);
416 }
417 action
418 }
419
420 /// Run the restoration phase, bracketed by the `PreRestoration` /
421 /// `PostRestoration` debug checkpoints so a debugger can inspect the
422 /// iterate just before entry and just after exit. With no debugger
423 /// installed this is exactly `invoke_restoration()`.
424 fn invoke_restoration_debugged(&mut self) -> IterateOutcome {
425 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::PreRestoration) {
426 return o;
427 }
428 let outcome = self.invoke_restoration();
429 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::PostRestoration) {
430 return o;
431 }
432 outcome
433 }
434
435 /// Fire a sub-iteration checkpoint from inside [`Self::iterate`].
436 /// Returns `Some(Terminate(UserRequestedStop))` if the debugger asked
437 /// to stop, so the caller can `return` it; `None` to continue.
438 fn debug_stop(&mut self, cp: crate::debug::Checkpoint) -> Option<IterateOutcome> {
439 if self.debug.is_none() {
440 return None;
441 }
442 if self.fire_debug(cp) == crate::debug::DebugAction::Stop {
443 Some(IterateOutcome::Terminate(SolverReturn::UserRequestedStop))
444 } else {
445 None
446 }
447 }
448
449 /// Fire the terminal post-mortem checkpoint (if a debugger is set),
450 /// carrying the solve outcome so the hook can decide whether to pause
451 /// at the final iterate. The action is advisory — the loop returns
452 /// `result` regardless — so the hook just gets a last look.
453 fn fire_debug_terminal(&mut self, result: SolverReturn) {
454 use crate::debug::{Checkpoint, DebugCtx};
455 let Some(hook) = self.debug.as_ref() else {
456 return;
457 };
458 let mut ctx = DebugCtx::new(
459 Rc::clone(&self.data),
460 Rc::clone(&self.cq),
461 Checkpoint::Terminated,
462 )
463 .with_status(format!("{result:?}"));
464 let _ = hook.borrow_mut().at_checkpoint(&mut ctx);
465 }
466
467 /// One iteration body — port of `Optimize()`'s inner loop.
468 /// Returns either `Continue` to keep iterating or a terminal
469 /// [`SolverReturn`] mirroring upstream's exception → return-code
470 /// translation table (see `MAIN_LOOP.md` §"Exception mapping").
471 fn iterate(&mut self) -> IterateOutcome {
472 // Shared timing accumulator — cheap Rc clone so each phase can
473 // bump its own counter without re-borrowing `data`.
474 let timing = self.data.borrow().timing.clone();
475
476 // Per-iteration span so every event emitted in this body (the
477 // structured iteration record, restoration/linear-solve spans)
478 // is tagged with the iteration index.
479 let _iter_span =
480 tracing::info_span!("iteration", iter = self.data.borrow().iter_count).entered();
481
482 // 1. Output iteration row. Header every 10 iters; the row itself
483 // is built plain by the strategy (so the column widths stay
484 // exact and unit-testable) and wrapped in a tiger/rust style
485 // at the print site (pounce#71). `anstream::stdout()` strips
486 // the escapes automatically when stdout is redirected or
487 // `NO_COLOR` is set, so non-TTY output is plain text.
488 //
489 // Print BEFORE `reset_info` so the row reflects the accepted
490 // step from the previous iteration (alphas, ls count,
491 // alpha_char), matching upstream's `IpIpoptAlgorithm::Optimize`
492 // ordering.
493 timing.output_iteration.start();
494 self.bundle.iter_output.write_output();
495 if self.print_iter_output {
496 use std::io::Write as _;
497 let (iter_count, alpha_pr, alpha_char) = {
498 let d = self.data.borrow();
499 (d.iter_count, d.info_alpha_primal, d.info_alpha_primal_char)
500 };
501 let row = self.bundle.iter_output.format_row(&self.data, &self.cq);
502 // Iteration 0 is the initial point — no step has been taken
503 // yet, so `alpha_primal` is 0; treat it as a full step
504 // (neutral black) rather than a stalling alarm (red).
505 let style_alpha = if iter_count == 0 { 1.0 } else { alpha_pr };
506 let style = pounce_common::style::iteration_row_style(style_alpha, alpha_char);
507 let mut out = anstream::stdout();
508 // Write errors (e.g. a closed pipe / `head` on the output)
509 // are deliberately ignored: a vanished terminal must not
510 // panic the solver, unlike the old `println!`.
511 if iter_count % 10 == 0 {
512 let _ = write!(out, "{}", crate::output::orig::OrigIterationOutput::HEADER);
513 }
514 let _ = writeln!(out, "{}{}{}", style.render(), row, style.render_reset());
515 }
516 timing.output_iteration.end();
517
518 // Structured per-iteration event (pounce#71) — the single source
519 // of truth for the per-iteration trajectory. The JSON log sink
520 // and the solve-report collector
521 // (`pounce_observability::IterCollectorLayer`) both derive from
522 // it. The text console layer filters this target out (its human
523 // form is the colored table above).
524 //
525 // Skipped entirely when nothing consumes it (no iter-history
526 // capture active and JSON logging off) so the default run pays
527 // no per-iteration field-evaluation / allocation cost.
528 if pounce_observability::iteration_event_wanted() {
529 let d = self.data.borrow();
530 let c = self.cq.borrow();
531 let alpha_char = d.info_alpha_primal_char;
532 let alpha_char_s = alpha_char.to_string();
533 let d_norm = match &d.delta {
534 Some(delta) => delta.x.amax().max(delta.s.amax()),
535 None => 0.0,
536 };
537 tracing::info!(
538 target: pounce_observability::ITER_TARGET,
539 iter = d.iter_count,
540 objective = c.unscaled_curr_f(),
541 inf_pr = c.curr_primal_infeasibility_max(),
542 inf_du = c.curr_dual_infeasibility_max(),
543 mu = d.curr_mu,
544 d_norm = d_norm,
545 regularization = d.info_regu_x,
546 alpha_dual = d.info_alpha_dual,
547 alpha_primal = d.info_alpha_primal,
548 ls_trials = d.info_ls_count,
549 alpha_char = alpha_char_s.as_str(),
550 resto_kind = pounce_common::style::resto_kind_str(alpha_char),
551 );
552 }
553
554 // Reset per-iteration info on data (after printing previous
555 // iter's accepted-step info; before the next line search).
556 self.data.borrow_mut().reset_info();
557
558 // 2. Convergence check.
559 timing.check_convergence.start();
560 let nlp_err = self.cq.borrow().curr_nlp_error();
561 let iter_count = self.data.borrow().iter_count;
562 if !nlp_err.is_finite() {
563 timing.check_convergence.end();
564 return IterateOutcome::Terminate(SolverReturn::InvalidNumberDetected);
565 }
566 // Divergence guard — port of upstream `IpIpoptAlg.cpp` post-
567 // AcceptTrialPoint check. When `max_i |x_i|` exceeds the
568 // registered `diverging_iterates_tol` (default `1e20`), exit
569 // cleanly with `DivergingIterates` rather than spiralling into
570 // a degenerate restoration whose inner sub-NLP can't recover
571 // (MESH: orig `f` already at -3.6e33 by iter 90, restoration
572 // entered too late to bound `x`).
573 if let Some(curr) = self.data.borrow().curr.as_ref() {
574 if curr.x.amax() > self.diverging_iterates_tol {
575 timing.check_convergence.end();
576 return IterateOutcome::Terminate(SolverReturn::DivergingIterates);
577 }
578 }
579 let conv_status = self
580 .bundle
581 .conv_check
582 .check_convergence_with_state(nlp_err, iter_count, &self.data, &self.cq);
583 match conv_status {
584 ConvergenceStatus::Continue => {}
585 ConvergenceStatus::Converged => {
586 timing.check_convergence.end();
587 return IterateOutcome::Terminate(SolverReturn::Success);
588 }
589 ConvergenceStatus::ConvergedToAcceptable => {
590 timing.check_convergence.end();
591 return IterateOutcome::Terminate(SolverReturn::StopAtAcceptablePoint);
592 }
593 ConvergenceStatus::MaxIterExceeded => {
594 timing.check_convergence.end();
595 return IterateOutcome::Terminate(SolverReturn::MaxiterExceeded);
596 }
597 ConvergenceStatus::CpuTimeExceeded => {
598 timing.check_convergence.end();
599 return IterateOutcome::Terminate(SolverReturn::CpuTimeExceeded);
600 }
601 ConvergenceStatus::WallTimeExceeded => {
602 timing.check_convergence.end();
603 return IterateOutcome::Terminate(SolverReturn::WallTimeExceeded);
604 }
605 ConvergenceStatus::LocallyInfeasible => {
606 timing.check_convergence.end();
607 return IterateOutcome::Terminate(SolverReturn::LocalInfeasibility);
608 }
609 ConvergenceStatus::Failed => {
610 timing.check_convergence.end();
611 return IterateOutcome::Terminate(SolverReturn::InternalError);
612 }
613 }
614
615 // Stash the iterate if it satisfies the per-component
616 // `acceptable_*_tol` triplet. Mirrors upstream
617 // `IpBacktrackingLineSearch.cpp:282-289` — checked at the top
618 // of every line-search call so the most recent acceptable
619 // iterate is always available as a rollback target if
620 // restoration later fails. The recorder feeds
621 // `acceptable_obj_change_tol`'s stability cross-check on
622 // subsequent iterates.
623 if self
624 .bundle
625 .conv_check
626 .current_is_acceptable_with_state(nlp_err, &self.data, &self.cq)
627 {
628 self.store_acceptable_point();
629 let curr_f = self.cq.borrow().curr_f();
630 self.bundle.conv_check.set_curr_acceptable_obj(curr_f);
631 }
632 timing.check_convergence.end();
633
634 // 3. Hessian update. Must run BEFORE `update_barrier_parameter`
635 // so the adaptive-μ oracles (probing, quality-function) drive
636 // their affine/centering solves against `W(curr_N)`, not the
637 // stale `W(curr_{N-1})` left in `data.w` by the previous iter's
638 // tail-end Hessian update. Upstream calls `UpdateHessian()`
639 // first in every main-loop body (`IpIpoptAlg.cpp:386`); pounce
640 // previously reordered this to the tail, which made iters 1+
641 // pick μ from the prior iterate's Hessian on adaptive-mu +
642 // quality-function — visible on CRESC50 as a catastrophic
643 // early-iter divergence (theta=5.8e5 by iter 61 vs upstream
644 // never entering restoration).
645 timing.update_hessian.start();
646 let _ = self.bundle.hess.update_hessian(&self.data, &self.cq);
647 timing.update_hessian.end();
648
649 // 4. Barrier parameter. Pass nlp + search_dir through so the
650 // adaptive μ oracles (probing, quality-function) can drive
651 // their own affine-step solves; monotone ignores them.
652 // Snapshot the tiny-step flag (set by the previous iteration's
653 // tiny-step branch) and the entry mu — if μ can't reduce while
654 // the flag is on, upstream `IpMonotoneMuUpdate.cpp:158-161`
655 // throws TINY_STEP_DETECTED → STOP_AT_TINY_STEP, which we
656 // realise as a clean termination here. Only the monotone update
657 // throws: `IpAdaptiveMuUpdate.cpp` consumes the tiny-step flag
658 // via its `force_no_progress` path (fix μ, keep iterating), so
659 // the termination is gated on `terminates_on_tiny_step()`.
660 timing.update_barrier_parameter.start();
661 let tiny_at_entry = self.data.borrow().tiny_step_flag;
662 let mu_before = self.data.borrow().curr_mu;
663 let mu_terminates_on_tiny = self.bundle.mu_update.terminates_on_tiny_step();
664 let next_mu = self.bundle.mu_update.update_barrier_parameter(
665 &self.data,
666 &self.cq,
667 self.nlp.as_ref(),
668 self.search_dir.as_mut(),
669 );
670 self.data.borrow_mut().curr_mu = next_mu;
671 timing.update_barrier_parameter.end();
672 if tiny_at_entry && mu_terminates_on_tiny && (next_mu - mu_before).abs() < Number::EPSILON {
673 return IterateOutcome::Terminate(SolverReturn::StopAtTinyStep);
674 }
675
676 // pounce#58 — iterate-quality guard for the probing oracle.
677 // The μ-update layer sets `request_resto` when the input
678 // iterate is too corrupted for the probing rule to produce a
679 // sane μ (see `mu/adaptive.rs` Probing dispatch). Restoration
680 // re-initialises the multipliers and gives the outer loop a
681 // clean iterate to continue from. When no restoration phase
682 // is configured (embedded callers, tests), emit a one-line
683 // notice and continue with the current μ — the guard has
684 // already prevented the destabilising 4-order μ jump.
685 let request_resto = {
686 let mut d = self.data.borrow_mut();
687 let f = d.request_resto;
688 d.request_resto = false;
689 f
690 };
691 if request_resto {
692 if self.restoration.is_some() {
693 return self.invoke_restoration_debugged();
694 } else {
695 tracing::warn!(target: "pounce::algorithm",
696 "[POUNCE] probing-oracle iterate-quality guard fired \
697 at iter {}, but no restoration phase is configured; \
698 continuing with μ={:.3e}.",
699 self.data.borrow().iter_count,
700 next_mu,
701 );
702 }
703 }
704
705 // Mirror upstream `IpAdaptiveMuUpdate.cpp:339, 386, 431` and
706 // `IpMonotoneMuUpdate.cpp:165`: every code path that *changes*
707 // μ calls `linesearch_->Reset()`, which clears the filter via
708 // `FilterLSAcceptor::Reset` (`IpFilterLSAcceptor.cpp:524-532`).
709 // Rationale: filter entries are computed against the current
710 // barrier — when μ changes, prior entries no longer apply and
711 // would over-constrain acceptance. The two upstream paths that
712 // do NOT reset (stay-fixed-no-decrease and fixed→free transition)
713 // both keep μ at curr_mu, so the `mu_changed` check captures
714 // the intended distinction.
715 if next_mu != mu_before {
716 self.bundle.line_search.reset();
717 }
718
719 // Sub-iteration checkpoint: μ has been updated for this iteration.
720 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::AfterBarrierUpdate) {
721 return o;
722 }
723
724 // 5. Search direction. Skipped without an NLP + search_dir.
725 // (Hessian was updated in step 3 above before the barrier-μ
726 // oracle so that adaptive-μ uses W(curr_N), not stale W.)
727 if let (Some(nlp), Some(sd)) = (self.nlp.as_ref(), self.search_dir.as_mut()) {
728 timing.compute_search_direction.start();
729 // Fields are declared `Empty` and filled by the linear
730 // solver (matrix size, factor nnz, inertia, ordering — see
731 // `pounce_feral::record_factor_stats`) and below
732 // (regularization), so the `linear_solve` span carries the
733 // KKT-solve characteristics for the JSON sink (pounce#71).
734 let ls_span = tracing::info_span!(
735 target: "pounce::linsol",
736 "linear_solve",
737 n = tracing::field::Empty,
738 matrix_nnz = tracing::field::Empty,
739 factor_nnz = tracing::field::Empty,
740 inertia_neg = tracing::field::Empty,
741 fill_ratio = tracing::field::Empty,
742 ordering = tracing::field::Empty,
743 regularization = tracing::field::Empty,
744 );
745 let ls_enter = ls_span.enter();
746 let ok = sd.compute_search_direction(&self.data, &self.cq, nlp);
747 ls_span.record("regularization", self.data.borrow().info_regu_x);
748 // Within-span marker so the enriched `linear_solve` fields
749 // (filled by the solver above) surface to the JSON sink at
750 // debug level; off at the default `info` level.
751 tracing::debug!(target: "pounce::linsol", "kkt solve complete");
752 drop(ls_enter);
753 timing.compute_search_direction.end();
754 if !ok {
755 // Mirror upstream `IpIpoptAlg.cpp:417-430`: a failed
756 // step computation puts the algorithm in emergency
757 // mode, which calls `BacktrackingLineSearch::
758 // ActivateFallbackMechanism` (cpp:1312-1328). When a
759 // restoration phase is configured, the next pass of
760 // `ComputeAcceptableTrialPoint` sees `goto_resto` at
761 // cpp:299-306 and hands control to restoration. Only
762 // when neither restoration nor an acceptor-level
763 // fallback is available does upstream throw
764 // `STEP_COMPUTATION_FAILED`.
765 if self.restoration.is_some() {
766 return self.invoke_restoration_debugged();
767 }
768 return IterateOutcome::Terminate(SolverReturn::ErrorInStepComputation);
769 }
770 if std::env::var_os("POUNCE_DBG_DELTA").is_some() {
771 let d = self.data.borrow();
772 let it = d.iter_count;
773 if let Some(delta) = d.delta.as_ref() {
774 use crate::iterates_vector::IteratesVector;
775 use pounce_linalg::{compound_vector::CompoundVector, Vector};
776 let dv: &IteratesVector = delta;
777 tracing::debug!(target: "pounce::algorithm",
778 "[PN_DELTA] iter={} mu={:.6e} dx_amax={:.6e} ds_amax={:.6e} dyc_amax={:.6e} dyd_amax={:.6e} dzL_amax={:.6e} dzU_amax={:.6e} dvL_amax={:.6e} dvU_amax={:.6e}",
779 it, d.curr_mu,
780 dv.x.amax(), dv.s.amax(), dv.y_c.amax(), dv.y_d.amax(),
781 dv.z_l.amax(), dv.z_u.amax(), dv.v_l.amax(), dv.v_u.amax()
782 );
783 if let Some(cdx) = dv.x.as_any().downcast_ref::<CompoundVector>() {
784 tracing::debug!(target: "pounce::algorithm",
785 "[PN_DELTA] iter={} dx_blocks_amax: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
786 it,
787 cdx.comp(0).amax(),
788 cdx.comp(1).amax(),
789 cdx.comp(2).amax(),
790 cdx.comp(3).amax(),
791 cdx.comp(4).amax(),
792 );
793 tracing::debug!(target: "pounce::algorithm",
794 "[PN_DELTA] iter={} dx_blocks_nrm2: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
795 it,
796 cdx.comp(0).nrm2(),
797 cdx.comp(1).nrm2(),
798 cdx.comp(2).nrm2(),
799 cdx.comp(3).nrm2(),
800 cdx.comp(4).nrm2(),
801 );
802 tracing::debug!(target: "pounce::algorithm",
803 "[PN_DELTA] iter={} dx_blocks_asum: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
804 it,
805 cdx.comp(0).asum(),
806 cdx.comp(1).asum(),
807 cdx.comp(2).asum(),
808 cdx.comp(3).asum(),
809 cdx.comp(4).asum(),
810 );
811 // Argmax of orig block via dot with sign — print first few values.
812 if let Some(dv_orig) =
813 cdx.comp(0)
814 .as_any()
815 .downcast_ref::<pounce_linalg::dense_vector::DenseVector>()
816 {
817 let v = dv_orig.values();
818 let mut imax = 0usize;
819 let mut amax = 0.0f64;
820 for (i, &x) in v.iter().enumerate() {
821 if x.abs() > amax {
822 amax = x.abs();
823 imax = i;
824 }
825 }
826 tracing::debug!(target: "pounce::algorithm",
827 "[PN_DELTA] iter={} dx_orig argmax: i={} v={:.17e} (n={})",
828 it,
829 imax,
830 v[imax],
831 v.len()
832 );
833 }
834 }
835 let p = &d.perturbations;
836 tracing::debug!(target: "pounce::algorithm",
837 "[PN_DELTA] iter={} pert: dx={:.6e} ds={:.6e} dc={:.6e} dd={:.6e}",
838 it, p.delta_x, p.delta_s, p.delta_c, p.delta_d
839 );
840 drop(d);
841 let cq = self.cq.borrow();
842 let gf = cq.curr_grad_f();
843 let gl = cq.curr_grad_lag_x();
844 let cc = cq.curr_c();
845 let cd = cq.curr_d_minus_s();
846 let sx = cq.curr_sigma_x();
847 let ss = cq.curr_sigma_s();
848 tracing::debug!(target: "pounce::algorithm",
849 "[PN_DELTA] iter={} cq: gradf_amax={:.6e} gradf_nrm2={:.6e} gradlag_amax={:.6e} gradlag_nrm2={:.6e} c_amax={:.6e} c_nrm2={:.6e} d_amax={:.6e} d_nrm2={:.6e} sigx_amax={:.6e} sigx_nrm2={:.6e} sigs_amax={:.6e} sigs_nrm2={:.6e}",
850 it,
851 gf.amax(), gf.nrm2(),
852 gl.amax(), gl.nrm2(),
853 cc.amax(), cc.nrm2(),
854 cd.amax(), cd.nrm2(),
855 sx.amax(), sx.nrm2(),
856 ss.amax(), ss.nrm2(),
857 );
858 if let Some(cgf) = gf.as_any().downcast_ref::<CompoundVector>() {
859 tracing::debug!(target: "pounce::algorithm",
860 "[PN_DELTA] iter={} gradf_blocks_amax: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
861 it,
862 cgf.comp(0).amax(),
863 cgf.comp(1).amax(),
864 cgf.comp(2).amax(),
865 cgf.comp(3).amax(),
866 cgf.comp(4).amax(),
867 );
868 }
869 if let Some(curr) = self.data.borrow().curr.clone() {
870 tracing::debug!(target: "pounce::algorithm",
871 "[PN_DELTA] iter={} bound_mults: zL_amax={:.6e} zU_amax={:.6e} vL_amax={:.6e} vU_amax={:.6e} s_amax={:.6e} s_nrm2={:.6e} x_amax={:.6e} x_nrm2={:.6e}",
872 it,
873 curr.z_l.amax(), curr.z_u.amax(),
874 curr.v_l.amax(), curr.v_u.amax(),
875 curr.s.amax(), curr.s.nrm2(),
876 curr.x.amax(), curr.x.nrm2(),
877 );
878 if let Some(czl) = curr.z_l.as_any().downcast_ref::<CompoundVector>() {
879 tracing::debug!(target: "pounce::algorithm",
880 "[PN_DELTA] iter={} zL_blocks_amax: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
881 it,
882 czl.comp(0).amax(),
883 czl.comp(1).amax(),
884 czl.comp(2).amax(),
885 czl.comp(3).amax(),
886 czl.comp(4).amax(),
887 );
888 }
889 if let Some(czu) = curr.z_u.as_any().downcast_ref::<CompoundVector>() {
890 tracing::debug!(target: "pounce::algorithm", "[PN_DELTA] iter={} zU_ncomps={}", it, czu.n_comps());
891 for ic in 0..czu.n_comps() {
892 tracing::debug!(target: "pounce::algorithm",
893 "[PN_DELTA] iter={} zU_block[{}]_amax={:.6e} dim={}",
894 it,
895 ic,
896 czu.comp(ic).amax(),
897 czu.comp(ic).dim()
898 );
899 }
900 }
901 }
902 if let Some(csx) = sx.as_any().downcast_ref::<CompoundVector>() {
903 tracing::debug!(target: "pounce::algorithm",
904 "[PN_DELTA] iter={} sigx_blocks_amax: orig={:.6e} nc={:.6e} pc={:.6e} nd={:.6e} pd={:.6e}",
905 it,
906 csx.comp(0).amax(),
907 csx.comp(1).amax(),
908 csx.comp(2).amax(),
909 csx.comp(3).amax(),
910 csx.comp(4).amax(),
911 );
912 }
913 drop(cq);
914 let d = self.data.borrow();
915 // Also dump curr.x_orig argmax
916 if let Some(curr) = d.curr.as_ref() {
917 if let Some(cx) = curr.x.as_any().downcast_ref::<CompoundVector>() {
918 if let Some(xo) =
919 cx.comp(0)
920 .as_any()
921 .downcast_ref::<pounce_linalg::dense_vector::DenseVector>()
922 {
923 let v = xo.values();
924 let mut imax = 0usize;
925 let mut amax = 0.0f64;
926 for (i, &x) in v.iter().enumerate() {
927 if x.abs() > amax {
928 amax = x.abs();
929 imax = i;
930 }
931 }
932 tracing::debug!(target: "pounce::algorithm", "[PN_DELTA] iter={} curr_x_orig argmax: i={} v={:.17e} amax={:.17e} nrm2={:.17e}",
933 it, imax, v[imax], xo.amax(), xo.nrm2());
934 }
935 }
936 }
937 }
938 }
939 }
940
941 // Capture KKT-factorization diagnostics for the debugger before
942 // the line search runs. Only when a debugger is installed. The
943 // inertia/status fields are cheap and always captured; the matrix
944 // triplets and `LDLᵀ` factor are O(nnz) assemblies, so they're
945 // captured only while the debugger is stepping (`wants_kkt_capture`)
946 // — a detached/free-running debugger drops them to keep the run
947 // cheap. `kkt_debug` is overwritten every iteration and never
948 // cleared at `iter_start`, so a stepping session always has the
949 // previous iteration's system to look back at via `viz kkt`/`viz L`.
950 if let Some(hook) = self.debug.as_ref() {
951 let capture_heavy = hook.borrow().wants_kkt_capture();
952 let captured_iter = self.data.borrow().iter_count;
953 let info = self.search_dir.as_ref().map(|sd| {
954 let pd = sd.pd_solver_mut();
955 let aug = pd.aug_solver();
956 let provides = aug.provides_inertia();
957 crate::ipopt_data::KktDebug {
958 iter: captured_iter,
959 dim: aug.system_dim(),
960 n_neg: if provides {
961 aug.number_of_neg_evals()
962 } else {
963 -1
964 },
965 provides_inertia: provides,
966 status: format!("{:?}", aug.last_solve_status()),
967 matrix: if capture_heavy {
968 aug.kkt_triplets()
969 } else {
970 None
971 },
972 l_factor: if capture_heavy {
973 aug.l_factor(true)
974 } else {
975 None
976 },
977 }
978 });
979 self.data.borrow_mut().kkt_debug = info;
980 }
981
982 // Sub-iteration checkpoint: the Newton step `δ` (data.delta) and
983 // the applied regularization are now available, before the line
984 // search consumes them.
985 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::AfterSearchDirection) {
986 return o;
987 }
988
989 // 6. Acceptable trial point — run the line search if we have a
990 // primal/dual step on `data.delta`. Wrap in a guard so all
991 // early-return paths (ErrorInStepComputation, InternalError,
992 // restoration entry) still stop the timer.
993 let _ls_guard = timing.compute_acceptable_trial_point.guard();
994 let have_delta = self.data.borrow().delta.is_some();
995 if have_delta {
996 let delta = match self.data.borrow().delta.as_ref().cloned() {
997 Some(d) => d,
998 None => {
999 return IterateOutcome::Terminate(SolverReturn::ErrorInStepComputation);
1000 }
1001 };
1002 // Cap alpha by the primal fraction-to-the-boundary so the
1003 // first trial cannot push slacks past their bounds, and by
1004 // the dual FTB so bound multipliers stay positive. Mirrors
1005 // upstream `IpBacktrackingLineSearch::FindAcceptableTrialPoint`'s
1006 // calls to `IpCq.primal_frac_to_the_bound` /
1007 // `IpCq.dual_frac_to_the_bound` with τ = `curr_tau`.
1008 let tau = self.data.borrow().curr_tau;
1009 let alpha_p_max = self.cq.borrow().aff_step_alpha_primal_max(&delta, tau);
1010 let alpha_d_max = self.cq.borrow().aff_step_alpha_dual_max(&delta, tau);
1011
1012 // Tiny-step gate — port of `IpBacktrackingLineSearch.cpp:363`
1013 // and the handling block at lines 382-435. When the search
1014 // direction is so small that any nonzero α would just
1015 // bounce inside floating-point noise, take the FTB step
1016 // unchecked and skip the line search; that's the only way
1017 // to hit `STOP_AT_TINY_STEP` cleanly when the iterate is
1018 // already at a converged point but `nlp_error > tol` due to
1019 // scaling or unbounded duals.
1020 if self.detect_tiny_step(&delta) {
1021 let alpha_p = alpha_p_max;
1022 let alpha_d = alpha_d_max;
1023 let curr = match self.data.borrow().curr.clone() {
1024 Some(c) => c,
1025 None => return IterateOutcome::Terminate(SolverReturn::InternalError),
1026 };
1027 let trial_iv = scaled_step_unchecked(&curr, &delta, alpha_p, alpha_d);
1028 {
1029 let mut d = self.data.borrow_mut();
1030 d.set_trial(trial_iv);
1031 d.info_alpha_primal = alpha_p;
1032 d.info_alpha_dual = alpha_d;
1033 d.info_ls_count = 0;
1034 if self.tiny_step_last_iteration {
1035 d.info_alpha_primal_char = 'T';
1036 d.tiny_step_flag = true;
1037 } else {
1038 d.info_alpha_primal_char = 't';
1039 }
1040 }
1041 let dy_amax = delta.y_c.amax().max(delta.y_d.amax());
1042 self.tiny_step_last_iteration = dy_amax < self.tiny_step_y_tol;
1043 } else {
1044 self.tiny_step_last_iteration = false;
1045 let alpha_init = self.alpha_init.min(alpha_p_max);
1046 let alpha_dual = self.alpha_init.min(alpha_d_max);
1047 let outcome = self.bundle.line_search.find_acceptable_trial_point(
1048 &self.data,
1049 &self.cq,
1050 &delta,
1051 alpha_init,
1052 alpha_dual,
1053 self.nlp.as_ref(),
1054 self.search_dir.as_mut(),
1055 );
1056 match outcome {
1057 Outcome::Accepted => {
1058 // A normal LS-accepted step breaks any in-flight
1059 // restoration cycle — clear the cycle detector
1060 // so the next resto entry starts fresh.
1061 self.last_resto_entry_x = None;
1062 self.last_resto_entry_s = None;
1063 self.last_resto_recovery_x = None;
1064 self.last_resto_recovery_s = None;
1065 self.resto_no_outer_progress_count = 0;
1066 // Intentionally *not* clearing
1067 // `resto_near_feasible_count` here: DECONVBNE's
1068 // cycle interleaves R-recoveries with 2-3
1069 // LS-accepted 'f'/'h' steps (which return
1070 // `Outcome::Accepted` but accomplish no real
1071 // outer progress — alpha drops to 1e-6 and
1072 // inf_du remains pinned at 1.9e7), so resetting
1073 // on every accept would zero the counter every
1074 // cycle and never fire. The counter persists
1075 // for the duration of the run and trips after
1076 // 3 cumulative near-feasible entries; legitimate
1077 // solves enter resto at most once at near-
1078 // feasibility (POLAK6, HAIFAM) and stay under
1079 // the limit.
1080 }
1081 Outcome::TinyStep | Outcome::Failed => {
1082 // Debugger stop: the line search rejected the step
1083 // (tiny-step floor or all backtracks failed), before
1084 // we fall into restoration. Lets a "why did the line
1085 // search give up?" inspection happen at the failing
1086 // point distinctly from the restoration entry.
1087 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::StepRejected) {
1088 return o;
1089 }
1090 // Upstream `IpBacktrackingLineSearch.cpp` raises
1091 // `LINE_SEARCH_FAILED` when α drops below
1092 // `alpha_min` or all retries reject, which in
1093 // turn triggers `ActivateLineSearch` →
1094 // restoration.
1095 return self.invoke_restoration_debugged();
1096 }
1097 }
1098 }
1099 }
1100
1101 // End the line-search/trial timer here so the bookkeeping in
1102 // steps 7-8 below is attributed to `accept_trial_point` (which
1103 // mirrors upstream's split: filter update and FTB reset are
1104 // accept-side, not line-search-side).
1105 _ls_guard.stop();
1106
1107 // 7. Accept trial point (promotes `trial` to `curr` if set).
1108 // The acceptor's filter has already been augmented (when
1109 // appropriate) inside `find_acceptable_trial_point` via
1110 // `update_for_next_iteration`, mirroring upstream's call
1111 // chain in `IpBacktrackingLineSearch.cpp:839`.
1112 let _accept_guard = timing.accept_trial_point.guard();
1113
1114 // 7a. Safe-slack bound adjustment. Before promoting `trial`, move
1115 // any `x_L/x_U/d_L/d_U` whose trial slack fell below
1116 // `eps*min(1,mu)` so the slack becomes representable (port of
1117 // the bound-adjustment block in
1118 // `IpoptAlgorithm::AcceptTrialPoint`, `IpIpoptAlg.cpp:664-706`).
1119 self.adjust_variable_bounds_for_small_slacks();
1120
1121 self.data.borrow_mut().accept_trial_point();
1122
1123 // 8. Bound multiplier kappa_sigma reset.
1124 self.correct_bound_multiplier();
1125
1126 // Sub-iteration checkpoint: the trial point was accepted; α and
1127 // the new iterate are in place (before the loop's iter bookkeeping
1128 // and the next `IterStart`).
1129 drop(_accept_guard);
1130 if let Some(o) = self.debug_stop(crate::debug::Checkpoint::AfterStep) {
1131 return o;
1132 }
1133
1134 IterateOutcome::Continue
1135 }
1136
1137 /// Port of `IpBacktrackingLineSearch::DetectTinyStep`
1138 /// (`IpBacktrackingLineSearch.cpp:1219-1278`). Returns true iff
1139 /// `max_i |δx_i|/(1+|x_i|) ≤ tiny_step_tol`,
1140 /// `max_i |δs_i|/(1+|s_i|) ≤ tiny_step_tol`, AND
1141 /// `curr_constraint_violation ≤ 1e-4`. Disabled when
1142 /// `tiny_step_tol == 0`.
1143 fn detect_tiny_step(&self, delta: &crate::iterates_vector::IteratesVector) -> bool {
1144 if self.tiny_step_tol == 0.0 {
1145 return false;
1146 }
1147 let curr = match self.data.borrow().curr.clone() {
1148 Some(c) => c,
1149 None => return false,
1150 };
1151
1152 // |x_i|+1
1153 let mut tmp = curr.x.make_new_copy();
1154 tmp.element_wise_abs();
1155 tmp.add_scalar(1.0);
1156 // |δx_i|/(|x_i|+1) ; checked via Amax of (δx ./ (|x|+1)).
1157 let mut tmp2 = delta.x.make_new_copy();
1158 tmp2.element_wise_divide(&*tmp);
1159 if tmp2.amax() > self.tiny_step_tol {
1160 return false;
1161 }
1162
1163 if curr.s.dim() > 0 {
1164 let mut tmp = curr.s.make_new_copy();
1165 tmp.element_wise_abs();
1166 tmp.add_scalar(1.0);
1167 let mut tmp2 = delta.s.make_new_copy();
1168 tmp2.element_wise_divide(&*tmp);
1169 if tmp2.amax() > self.tiny_step_tol {
1170 return false;
1171 }
1172 }
1173
1174 let cviol = self.cq.borrow().curr_constraint_violation();
1175 if cviol > 1e-4 {
1176 return false;
1177 }
1178 true
1179 }
1180
1181 /// Drive the restoration phase after a line-search failure.
1182 /// Returns `IterateOutcome::Continue` if the restoration driver
1183 /// recovered (the algorithm carries on from the recovered iterate);
1184 /// otherwise terminates with [`SolverReturn::RestorationFailure`].
1185 /// Mirrors upstream's
1186 /// `IpBacktrackingLineSearch::ActivateLineSearch` → `PerformRestoration`
1187 /// chain.
1188 fn invoke_restoration(&mut self) -> IterateOutcome {
1189 // Snapshot the outer reference iterate's `(theta, barr)` and
1190 // build the orig-progress callback the inner IPM will consult
1191 // at every iteration (mirrors upstream
1192 // `IpRestoFilterConvCheck::SetOrigLSAcceptor` plus
1193 // `IpFilterLSAcceptor::Reset`'s `reference_*_` snapshot).
1194 let reference_theta = self.cq.borrow().curr_constraint_violation();
1195 let reference_barr = self.cq.borrow().curr_barrier_obj();
1196
1197 if std::env::var("POUNCE_DBG_RESTO").is_ok() {
1198 let iter = self.data.borrow().iter_count;
1199 tracing::debug!(target: "pounce::algorithm",
1200 "RESTO_ENTRY iter={} theta={:.6e} barr={:.6e} near_feas_ct={}",
1201 iter, reference_theta, reference_barr, self.resto_near_feasible_count,
1202 );
1203 }
1204
1205 // No-progress restoration cycle detector. Two layered checks
1206 // surface as `ErrorInStepComputation` instead of cycling to
1207 // `max_iter` exhaustion (mirrors the *intent* of upstream
1208 // `IpBacktrackingLineSearch.cpp:580-600`'s almost-feasible
1209 // resto guard):
1210 //
1211 // 1. *Static cycle*: entry-to-entry — when the curr `(x, s)`
1212 // at this entry is essentially identical to the snapshot
1213 // from the previous entry, the inner resto-IPM is
1214 // returning recovered iterates indistinguishable from
1215 // entry, AND the outer didn't move either. Fires
1216 // immediately. Catches QCNEW, EQC, MESH, POLAK6, S365,
1217 // S365MOD, SIPOW2M, PFIT4.
1218 //
1219 // 2. *Slow-progress cycle*: recovery-to-entry — when curr at
1220 // this entry is essentially identical to the *recovery*
1221 // iterate from the previous resto, the outer made no
1222 // progress between resto invocations even though resto's
1223 // inner moved substantively. Counted, fires after 5
1224 // consecutive entries. Catches ACOPR14, ACOPR30, TRO3X3
1225 // while letting MAKELA3, HAIFAM, HALDMADS, ROBOT,
1226 // TENBARS2 — which need 2-3 productive resto entries
1227 // before LS accepts — pass through.
1228 //
1229 // A productive single-restoration sequence (BT8, HIMMELBJ,
1230 // LINSPANH, LSNNODOC, ODFITS, OET3) clears both snapshots via
1231 // `Outcome::Accepted` between entries and is unaffected.
1232 let curr = self
1233 .data
1234 .borrow()
1235 .curr
1236 .as_ref()
1237 .expect("curr set before invoke_restoration")
1238 .clone();
1239 // Helper: when the cycle detector fires and the orig cv is
1240 // bounded away from outer tol (e.g. PFIT1's 2.73e-2), the
1241 // outer is stuck at a feasibility-stationary point and the
1242 // honest exit is `LocalInfeasibility`. Below that threshold
1243 // the failure is numerical, not algorithmic, so retain
1244 // `ErrorInStepComputation`.
1245 let outer_tol_for_cycle = self.bundle.conv_check.tol_or_default();
1246 let cycle_exit = if reference_theta > (100.0 * outer_tol_for_cycle).max(1e-4) {
1247 SolverReturn::LocalInfeasibility
1248 } else {
1249 SolverReturn::ErrorInStepComputation
1250 };
1251 let static_cycle = if let (Some(prev_x), Some(prev_s)) = (
1252 self.last_resto_entry_x.as_ref(),
1253 self.last_resto_entry_s.as_ref(),
1254 ) {
1255 let dx_rel = relative_distance(&*curr.x, &**prev_x);
1256 let ds_rel = relative_distance(&*curr.s, &**prev_s);
1257 if std::env::var_os("POUNCE_DBG_RESTO_CYCLE").is_some() {
1258 tracing::debug!(target: "pounce::algorithm",
1259 "[PN_RESTO_CYCLE] entry-vs-entry dx_rel={:.6e} ds_rel={:.6e}",
1260 dx_rel, ds_rel
1261 );
1262 }
1263 dx_rel <= 1e-10 && ds_rel <= 1e-10
1264 } else {
1265 false
1266 };
1267 if static_cycle {
1268 // Prefer the last acceptable point over the cycle error —
1269 // the borrows above are released, so the `&mut self` helper
1270 // is free to roll back.
1271 return self.terminate_acceptable_or(cycle_exit);
1272 }
1273 let recovery_cycle = if let (Some(prev_x), Some(prev_s)) = (
1274 self.last_resto_recovery_x.as_ref(),
1275 self.last_resto_recovery_s.as_ref(),
1276 ) {
1277 let dx_rel = relative_distance(&*curr.x, &**prev_x);
1278 let ds_rel = relative_distance(&*curr.s, &**prev_s);
1279 if std::env::var_os("POUNCE_DBG_RESTO_CYCLE").is_some() {
1280 tracing::debug!(target: "pounce::algorithm",
1281 "[PN_RESTO_CYCLE] entry-vs-recovery dx_rel={:.6e} ds_rel={:.6e} count={}",
1282 dx_rel, ds_rel, self.resto_no_outer_progress_count
1283 );
1284 }
1285 dx_rel <= 1e-10 && ds_rel <= 1e-10
1286 } else {
1287 false
1288 };
1289 if recovery_cycle {
1290 self.resto_no_outer_progress_count =
1291 self.resto_no_outer_progress_count.saturating_add(1);
1292 // 10-strike limit: tuned to give OET7-style traces room
1293 // to break through (inner inf_pr still decreasing across
1294 // strikes) while still bounding DECONVBNE-style cycles
1295 // (which need a guard but tolerate a wider window —
1296 // ~3 outer steps per cycle, so 10 strikes ≈ 30 outer
1297 // iters, well below the 2987-iter pathological run).
1298 if self.resto_no_outer_progress_count >= 10 {
1299 // Prefer the last acceptable point over the cycle error;
1300 // borrows are released, so the `&mut self` helper is free.
1301 return self.terminate_acceptable_or(cycle_exit);
1302 }
1303 } else {
1304 self.resto_no_outer_progress_count = 0;
1305 }
1306 // Near-feasible resto re-entry detector — matches the *intent*
1307 // of upstream `IpBacktrackingLineSearch.cpp:580-600`'s almost-
1308 // feasible-resto guard with a looser cv threshold. When the
1309 // outer enters restoration with the constraint violation
1310 // already at or below `tol`, the resto sub-IPM will produce a
1311 // recovered iterate that's at most marginally more feasible,
1312 // and any post-recovery σ-blowup from the next outer KKT solve
1313 // will re-trigger resto on the next iteration. Counting these
1314 // entries surfaces the cycle as `StopAtAcceptablePoint` —
1315 // primal feasibility is already met, only the dual residual
1316 // remains. Catches DECONVBNE: pounce ran 2987 iters before
1317 // this guard (cycle of ~30-inner-resto + 3 outer per cycle);
1318 // upstream solves in 505 iters via a different x trajectory.
1319 // Single-entry productive restos (BT8, HIMMELBJ, ODFITS) and
1320 // sub-tol-but-recoverable starts pass through under the 3-
1321 // strike limit.
1322 let outer_tol = self.bundle.conv_check.tol_or_default();
1323 if reference_theta <= outer_tol {
1324 self.resto_near_feasible_count = self.resto_near_feasible_count.saturating_add(1);
1325 if self.resto_near_feasible_count >= 3 {
1326 // Constraint feasibility is met, but a near-feasible iterate is
1327 // only "acceptable" if its objective is finite. CUTE `himmelbj`
1328 // reaches a point with cv ≈ 2e-9 where f evaluates to NaN; that
1329 // must surface as Invalid_Number_Detected rather than be
1330 // reported as Solved_To_Acceptable_Level with a `nan` objective.
1331 if !self.cq.borrow().curr_f().is_finite() {
1332 return IterateOutcome::Terminate(SolverReturn::InvalidNumberDetected);
1333 }
1334 return IterateOutcome::Terminate(SolverReturn::StopAtAcceptablePoint);
1335 }
1336 } else {
1337 self.resto_near_feasible_count = 0;
1338 }
1339 self.last_resto_entry_x = Some(curr.x.make_new_copy());
1340 self.last_resto_entry_s = Some(curr.s.make_new_copy());
1341
1342 // Augment the outer's filter with the resto-entry envelope —
1343 // mirrors upstream `IpBacktrackingLineSearch.cpp:566`:
1344 // `acceptor_->PrepareRestoPhaseStart()`. Adds
1345 // `((1-γ_θ)·θ_entry, φ_entry - γ_φ·θ_entry)` to the filter so
1346 // that after restoration recovers, the outer's Newton step is
1347 // forced by the filter to make real progress vs the entry
1348 // point. Without this, the outer accepts null-progress 'h'
1349 // steps and re-enters restoration on the next iteration (root
1350 // cause of DECONVBNE's 323 R-accepts vs ipopt's 21).
1351 self.bundle
1352 .line_search
1353 .acceptor_mut()
1354 .prepare_resto_phase_start(reference_theta, reference_barr);
1355
1356 let orig_progress_cb = self.bundle.line_search.acceptor().make_orig_progress_check(
1357 reference_theta,
1358 reference_barr,
1359 5.0,
1360 );
1361
1362 let (Some(nlp), Some(sd), Some(resto)) = (
1363 self.nlp.as_ref(),
1364 self.search_dir.as_mut(),
1365 self.restoration.as_mut(),
1366 ) else {
1367 return IterateOutcome::Terminate(SolverReturn::RestorationFailure);
1368 };
1369 resto.set_orig_progress_check(orig_progress_cb);
1370 // Forward the shared debugger so it can step the inner solve.
1371 resto.set_debug_hook(self.debug.as_ref().map(Rc::clone));
1372 let mut pd_guard = sd.pd_solver_mut();
1373 let aug = pd_guard.aug_solver_mut();
1374 // Audit counters (pounce#12). Increment call count + outer-iter
1375 // count (one outer iter is consumed per restoration call) and
1376 // wall-time around the inner call. Inner iter count is read
1377 // after via the trait accessor.
1378 self.resto_calls = self.resto_calls.saturating_add(1);
1379 self.resto_outer_iters = self.resto_outer_iters.saturating_add(1);
1380 let resto_t0 = std::time::Instant::now();
1381 let outcome = resto.perform_restoration(&self.data, &self.cq, nlp, aug);
1382 drop(pd_guard);
1383 self.resto_wall_secs += resto_t0.elapsed().as_secs_f64();
1384 self.resto_inner_iters = self
1385 .resto_inner_iters
1386 .saturating_add(resto.last_inner_iter_count());
1387 match outcome {
1388 RestorationOutcome::Recovered => {
1389 // The driver has staged the recovered point on
1390 // `data.trial`; apply the safe-slack bound adjustment
1391 // (as the main accept path does), then promote it and
1392 // continue iterating.
1393 self.adjust_variable_bounds_for_small_slacks();
1394 self.data.borrow_mut().accept_trial_point();
1395 // Snapshot the recovery iterate for the slow-cycle
1396 // detector at the top of the next `invoke_restoration`.
1397 // Compared against next-entry curr, dx_rel ≈ ‖α·d‖ —
1398 // measures purely the outer step. See header comment
1399 // on the cycle detector above.
1400 let recovered = self
1401 .data
1402 .borrow()
1403 .curr
1404 .as_ref()
1405 .expect("accept_trial_point sets curr")
1406 .clone();
1407 self.last_resto_recovery_x = Some(recovered.x.make_new_copy());
1408 self.last_resto_recovery_s = Some(recovered.s.make_new_copy());
1409 // Mirror upstream `IpoptAlgorithm::AcceptTrialPoint`
1410 // (`IpIpoptAlg.cpp:917-963`): kappa_sigma clamp on the
1411 // four bound-multiplier vectors. Upstream applies this
1412 // unconditionally inside AcceptTrialPoint, so the
1413 // post-restoration path inherits it; pounce factored
1414 // the clamp out of the data swap so we must call it
1415 // explicitly here. Without it the all-1 multiplier
1416 // reset (`bound_mult_reset_threshold`) leaves z*s far
1417 // from mu at the recovered iterate, blowing up the
1418 // next KKT solve's σ = z/s diagonal.
1419 self.correct_bound_multiplier();
1420 IterateOutcome::Continue
1421 }
1422 RestorationOutcome::Failed => {
1423 // Mirrors upstream `IpBacktrackingLineSearch.cpp:611-623`:
1424 // when `PerformRestoration` returns false, attempt to
1425 // roll back to the most recent acceptable iterate before
1426 // surfacing failure. If a snapshot is available we exit
1427 // cleanly with `StopAtAcceptablePoint` (mapped by the
1428 // application layer to `Solved_To_Acceptable_Level`),
1429 // matching the upstream `ACCEPTABLE_POINT_REACHED`
1430 // throw. Without a snapshot we surface
1431 // `RestorationFailure` — unless the restoration left the
1432 // iterate diverging (`|x|_∞ > diverging_iterates_tol`), in
1433 // which case we surface `DivergingIterates` to mirror the
1434 // outcome upstream produces on pathological problems like
1435 // MESH (where ipopt reports `Diverging_Iterates` and
1436 // pounce previously reported `Restoration_Failed` with an
1437 // obj of −3.6e+33).
1438 if self.restore_acceptable_point() {
1439 IterateOutcome::Terminate(SolverReturn::StopAtAcceptablePoint)
1440 } else if let Some(curr) = self.data.borrow().curr.as_ref() {
1441 if curr.x.amax() > self.diverging_iterates_tol {
1442 IterateOutcome::Terminate(SolverReturn::DivergingIterates)
1443 } else {
1444 IterateOutcome::Terminate(SolverReturn::RestorationFailure)
1445 }
1446 } else {
1447 IterateOutcome::Terminate(SolverReturn::RestorationFailure)
1448 }
1449 }
1450 RestorationOutcome::LocallyInfeasible => {
1451 // Mirrors upstream's catch of `LOCALLY_INFEASIBLE` thrown
1452 // from `IpRestoConvCheck.cpp:240` — the resto sub-IPM
1453 // settled at a stationary point of `||c(x)||_1` whose
1454 // residual is still well above `tol`. Without this
1455 // detection the outer would re-enter restoration on the
1456 // unchanged iterate forever.
1457 IterateOutcome::Terminate(SolverReturn::LocalInfeasibility)
1458 }
1459 }
1460 }
1461
1462 /// Safe-slack bound adjustment, applied to the staged `trial`
1463 /// iterate before it is promoted to `curr`. When one or more trial
1464 /// slacks fell below `eps*min(1,mu)`, [`IpoptCalculatedQuantities::
1465 /// adjusted_trial_bounds`] returns the moved `x_L/x_U/d_L/d_U`; we
1466 /// install them on the NLP so the slack becomes representable. Port
1467 /// of the bound-adjustment block in `IpoptAlgorithm::AcceptTrialPoint`
1468 /// (`IpIpoptAlg.cpp:664-706`).
1469 fn adjust_variable_bounds_for_small_slacks(&mut self) {
1470 // Compute the moved bounds (releases the CQ/NLP borrows on return).
1471 let adjusted = {
1472 let trial_set = self.data.borrow().trial.is_some();
1473 if !trial_set {
1474 return;
1475 }
1476 self.cq.borrow().adjusted_trial_bounds()
1477 };
1478 let Some(bounds) = adjusted else {
1479 return;
1480 };
1481 tracing::debug!(
1482 target: "pounce::algorithm",
1483 "slack_move: {} slack(s) too small, adjusting variable bound(s) at iter {}",
1484 bounds.adjusted,
1485 self.data.borrow().iter_count,
1486 );
1487 let nlp = Rc::clone(self.cq.borrow().nlp());
1488 nlp.borrow_mut().adjust_variable_bounds(
1489 &*bounds.x_l,
1490 &*bounds.x_u,
1491 &*bounds.d_l,
1492 &*bounds.d_u,
1493 );
1494 }
1495
1496 /// Port of `IpIpoptAlg::correct_bound_multiplier`
1497 /// (`IpIpoptAlg.cpp:1055-1134`). Clamp each bound multiplier
1498 /// component into `[mu/(kappa_sigma * s_i), kappa_sigma * mu / s_i]`
1499 /// for all four bound-multiplier vectors.
1500 fn correct_bound_multiplier(&mut self) {
1501 if self.kappa_sigma < 1.0 {
1502 return;
1503 }
1504 let mu = self.data.borrow().curr_mu;
1505 let curr = match self.data.borrow().curr.clone() {
1506 Some(c) => c,
1507 None => return,
1508 };
1509
1510 let cq = self.cq.borrow();
1511
1512 let z_l_new = clamp_against_slack(&*curr.z_l, &*cq.curr_slack_x_l(), mu, self.kappa_sigma);
1513 let z_u_new = clamp_against_slack(&*curr.z_u, &*cq.curr_slack_x_u(), mu, self.kappa_sigma);
1514 let v_l_new = clamp_against_slack(&*curr.v_l, &*cq.curr_slack_s_l(), mu, self.kappa_sigma);
1515 let v_u_new = clamp_against_slack(&*curr.v_u, &*cq.curr_slack_s_u(), mu, self.kappa_sigma);
1516 drop(cq);
1517
1518 let new_iv = crate::iterates_vector::IteratesVector::new(
1519 curr.x.clone(),
1520 curr.s.clone(),
1521 curr.y_c.clone(),
1522 curr.y_d.clone(),
1523 z_l_new,
1524 z_u_new,
1525 v_l_new,
1526 v_u_new,
1527 );
1528 self.data.borrow_mut().set_curr(new_iv);
1529 }
1530
1531 /// Outer entry point — port of `IpoptAlgorithm::Optimize()`. Calls
1532 /// the iterate-initializer once, then loops `iterate()` until a
1533 /// terminal status. The exception → SolverReturn mapping
1534 /// (TINY_STEP_DETECTED → STEP_BECOMES_TINY,
1535 /// RESTORATION_FAILED → RESTORATION_FAILURE, etc.) lands in
1536 /// Phase 9 alongside the restoration phase.
1537 pub fn optimize(&mut self) -> SolverReturn {
1538 // Top-level span for the whole solve; every iteration / linear
1539 // solve / restoration event nests under it (pounce#71).
1540 let _solve_span = tracing::info_span!("solve").entered();
1541
1542 // Shared timing accumulator — every phase below records into it.
1543 let timing = self.data.borrow().timing.clone();
1544
1545 // Install the shared accumulator on the augmented-system solver
1546 // so its factor / back-solve calls are attributed to
1547 // `linear_system_factorization` / `linear_system_back_solve`.
1548 // Same pattern for the diagnostics state when present, so KKT
1549 // dump sites can consult per-iter gating.
1550 if let Some(sd) = self.search_dir.as_mut() {
1551 sd.pd_solver_mut()
1552 .aug_solver_mut()
1553 .set_timing_stats(std::rc::Rc::clone(&timing));
1554 if let Some(diag) = self.diagnostics.as_ref() {
1555 sd.pd_solver_mut()
1556 .aug_solver_mut()
1557 .set_diagnostics(Rc::clone(diag));
1558 }
1559 }
1560
1561 // 0a. Strategy initialization — port of upstream's
1562 // `IpoptAlgorithm::InitializeImpl` calls. The mu update needs
1563 // `data.curr_mu`/`curr_tau` seeded before the iterate
1564 // initializer runs (`CalculateSafeSlack` reads them).
1565 self.bundle.mu_update.initialize(&self.data);
1566
1567 // 0b. Iterate initializer. Requires NLP; without one the caller
1568 // must have populated `data.curr` themselves.
1569 if let Some(nlp) = self.nlp.as_ref() {
1570 // The initializer needs an aug-system solver for the
1571 // least-square multiplier branch; until that's wired we
1572 // route through whatever the search-direction calculator
1573 // owns when present. For the stub flow we skip the LSM
1574 // path by giving the initializer a dummy solver only if
1575 // the search_dir is present (otherwise the init function
1576 // is responsible for not consulting it).
1577 if let Some(sd) = self.search_dir.as_mut() {
1578 timing.initialize_iterates.start();
1579 let mut pd_guard = sd.pd_solver_mut();
1580 let aug_solver = pd_guard.aug_solver_mut();
1581 let ok = self
1582 .bundle
1583 .init
1584 .set_initial_iterates(&self.data, &self.cq, nlp, aug_solver);
1585 drop(pd_guard);
1586 timing.initialize_iterates.end();
1587 if !ok {
1588 return SolverReturn::InternalError;
1589 }
1590 }
1591 }
1592
1593 // 0c. Seed `IpoptData::w` with the initial-iterate Hessian.
1594 // Redundant with the iter-body `update_hessian` call (which
1595 // now runs BEFORE `update_barrier_parameter`) but kept to
1596 // cover any code path that consults `data.w` between
1597 // `set_initial_iterates` and the first `iterate()` call
1598 // (e.g. the iter-0 trace dump below).
1599 if self.data.borrow().curr.is_some() {
1600 timing.update_hessian.start();
1601 let _ = self.bundle.hess.update_hessian(&self.data, &self.cq);
1602 timing.update_hessian.end();
1603 }
1604
1605 // Track-A iterate-trace dumper. Activated by
1606 // `IPOPT_ITER_DUMP_PATH`; otherwise no-op. See `iter_dump.rs`.
1607 let mut dumper = IterDumper::from_env();
1608 // Iter 0 record — captures the initialised iterate before any
1609 // step. Mirrors upstream's "after InitializeIterates(), before
1610 // the loop" emission point.
1611 if let Some(d) = dumper.as_mut() {
1612 d.write_record(&self.data, &self.cq);
1613 }
1614
1615 // Advance the diagnostics iter counter so the first `iterate()`
1616 // body reports as iter 0 (matches `data.iter_count`). Subsequent
1617 // bumps live at the bottom of the loop alongside the iter_count
1618 // bookkeeping.
1619 if let Some(diag) = self.diagnostics.as_ref() {
1620 diag.bump_iter();
1621 // Iter-0 iterate row (issue #68). Same hook point as
1622 // the binary IterDumper above; emits only when
1623 // `--dump iterates:*` is configured.
1624 emit_iterate_record(diag.as_ref(), &self.data, &self.cq);
1625 }
1626
1627 // Iter 0 intermediate callback — upstream fires once after
1628 // `InitializeIterates` before the loop body starts so users
1629 // observe the initial point.
1630 if !self.fire_intermediate() {
1631 return SolverReturn::UserRequestedStop;
1632 }
1633 if self.fire_debug(crate::debug::Checkpoint::IterStart) == crate::debug::DebugAction::Stop {
1634 return SolverReturn::UserRequestedStop;
1635 }
1636
1637 let result = loop {
1638 match self.iterate() {
1639 IterateOutcome::Terminate(ret) => break ret,
1640 IterateOutcome::Continue => {
1641 // Source the local counter from `data.iter_count`
1642 // each pass so a pre-seeded counter (e.g. the inner
1643 // restoration IPM at `outer.iter + 1`, matching
1644 // upstream `IpRestoMinC_1Nrm.cpp:181`) and any
1645 // restoration step that set
1646 // `data.iter_count = inner.iter_count - 1`
1647 // (mirroring `IpRestoMinC_1Nrm.cpp:Set_iter_count`)
1648 // are honored — without this the local counter
1649 // would advance from its pre-restoration value,
1650 // ignoring the inner-IPM iterations.
1651 let mut iter_count: Index = self.data.borrow().iter_count;
1652 iter_count += 1;
1653 if iter_count >= self.max_iter {
1654 break SolverReturn::MaxiterExceeded;
1655 }
1656 self.data.borrow_mut().iter_count = iter_count;
1657 // Keep the diagnostics counter in lock-step with
1658 // `data.iter_count` so KKT-dump gating reflects the
1659 // about-to-execute iteration.
1660 if let Some(diag) = self.diagnostics.as_ref() {
1661 diag.bump_iter();
1662 // Per-iter iterate row (issue #68). Mirrors
1663 // the binary IterDumper hook below.
1664 emit_iterate_record(diag.as_ref(), &self.data, &self.cq);
1665 }
1666 // Per-iteration record — emitted after the
1667 // iter_count bump so the recorded `iter` field
1668 // matches `IpData().iter_count()` at the moment of
1669 // emission, identical to upstream's writer.
1670 if let Some(d) = dumper.as_mut() {
1671 d.write_record(&self.data, &self.cq);
1672 }
1673 // Per-iteration intermediate callback — fired with
1674 // an `IntermediateContext` guard so downstream
1675 // inspector entry points (the C API
1676 // `GetIpoptCurrent*` family) see live state for the
1677 // duration of the user callback.
1678 if !self.fire_intermediate() {
1679 break SolverReturn::UserRequestedStop;
1680 }
1681 if self.fire_debug(crate::debug::Checkpoint::IterStart)
1682 == crate::debug::DebugAction::Stop
1683 {
1684 break SolverReturn::UserRequestedStop;
1685 }
1686 }
1687 }
1688 };
1689
1690 // Terminal post-mortem checkpoint. Skipped when the user already
1691 // asked to stop (they were just at a prompt); otherwise the
1692 // debugger gets a last look at the final/failing iterate.
1693 if !matches!(result, SolverReturn::UserRequestedStop) {
1694 self.fire_debug_terminal(result);
1695 }
1696 result
1697 }
1698}
1699
1700/// Internal result of one [`IpoptAlgorithm::iterate`] call. Mirrors the
1701/// upstream try/catch around `IpoptAlg::Optimize` — anything that's not
1702/// `Continue` carries the [`SolverReturn`] that the outer loop will
1703/// surface to `IpoptApplication`.
1704enum IterateOutcome {
1705 Continue,
1706 Terminate(SolverReturn),
1707}
1708
1709/// `||a - b||_2 / (1 + ||b||_2)`. Used by the restoration cycle
1710/// detector in [`IpoptAlgorithm::invoke_restoration`] to test whether
1711/// the outer iterate has moved between two consecutive restoration
1712/// entries.
1713fn relative_distance(a: &dyn Vector, b: &dyn Vector) -> Number {
1714 if a.dim() == 0 {
1715 return 0.0;
1716 }
1717 let mut diff = a.make_new_copy();
1718 diff.axpy(-1.0, b);
1719 diff.nrm2() / (1.0 + b.nrm2())
1720}
1721
1722/// `out = curr + α_p · δ` for the primal/equality blocks and
1723/// `out = curr + α_d · δ` for the bound multipliers, returned as a
1724/// fresh frozen `IteratesVector`. Mirrors `scaled_step` in the line
1725/// search; duplicated here for the tiny-step branch which bypasses
1726/// the line-search driver.
1727fn scaled_step_unchecked(
1728 curr: &crate::iterates_vector::IteratesVector,
1729 delta: &crate::iterates_vector::IteratesVector,
1730 alpha_primal: Number,
1731 alpha_dual: Number,
1732) -> crate::iterates_vector::IteratesVector {
1733 let mut out = curr.make_new_zeroed();
1734 out.add_one_vector(1.0, curr, 0.0);
1735 out.x.axpy(alpha_primal, &*delta.x);
1736 out.s.axpy(alpha_primal, &*delta.s);
1737 out.y_c.axpy(alpha_primal, &*delta.y_c);
1738 out.y_d.axpy(alpha_primal, &*delta.y_d);
1739 out.z_l.axpy(alpha_dual, &*delta.z_l);
1740 out.z_u.axpy(alpha_dual, &*delta.z_u);
1741 out.v_l.axpy(alpha_dual, &*delta.v_l);
1742 out.v_u.axpy(alpha_dual, &*delta.v_u);
1743 out.freeze()
1744}
1745
1746/// Allocate a fresh `Rc<dyn Vector>` with `kappa_sigma_clamp`
1747/// applied component-wise against the supplied `slack`. Inputs are
1748/// borrowed; the original `z` is never mutated. Ports the per-vector
1749/// piece of `IpIpoptAlg.cpp:1080-1133`.
1750fn clamp_against_slack(
1751 z: &dyn Vector,
1752 slack: &dyn Vector,
1753 mu: Number,
1754 kappa_sigma: Number,
1755) -> Rc<dyn Vector> {
1756 debug_assert_eq!(z.dim(), slack.dim());
1757 let n = z.dim() as usize;
1758 // Flatten both z and slack into contiguous slices so the
1759 // elementwise clamp doesn't care whether the inputs are
1760 // [`DenseVector`] (regular IPM path) or [`CompoundVector`]
1761 // (resto IPM path). The result is reconstructed into a
1762 // same-shape Vector via `Vector::make_new` + a flat-write
1763 // helper so the caller sees a vector with the same blocking as
1764 // its input.
1765 let mut buf = vec![0.0_f64; n];
1766 flat_read_into(z, &mut buf);
1767 let s_vals = flat_read_owned(slack);
1768 let _ = kappa_sigma_clamp(&mut buf, &s_vals, mu, kappa_sigma);
1769 let mut out: Box<dyn Vector> = z.make_new();
1770 flat_write_into(&mut *out, &buf);
1771 Rc::from(out)
1772}
1773
1774pub(crate) fn flat_read_into(v: &dyn Vector, dst: &mut [Number]) {
1775 if let Some(dv) = v
1776 .as_any()
1777 .downcast_ref::<pounce_linalg::dense_vector::DenseVector>()
1778 {
1779 let vs = dv.expanded_values();
1780 dst.copy_from_slice(&vs);
1781 return;
1782 }
1783 if let Some(cv) = v.as_any().downcast_ref::<pounce_linalg::CompoundVector>() {
1784 let mut off = 0usize;
1785 for k in 0..cv.n_comps() {
1786 let blk = cv.comp(k);
1787 let dim = blk.dim() as usize;
1788 let dblk = blk
1789 .as_any()
1790 .downcast_ref::<pounce_linalg::dense_vector::DenseVector>()
1791 .expect("clamp_against_slack: CompoundVector blocks must be DenseVectors");
1792 let vs = dblk.expanded_values();
1793 dst[off..off + dim].copy_from_slice(&vs);
1794 off += dim;
1795 }
1796 return;
1797 }
1798 panic!("clamp_against_slack: unsupported Vector kind");
1799}
1800
1801pub(crate) fn flat_read_owned(v: &dyn Vector) -> Vec<Number> {
1802 let mut out = vec![0.0; v.dim() as usize];
1803 flat_read_into(v, &mut out);
1804 out
1805}
1806
1807pub(crate) fn flat_write_into(v: &mut dyn Vector, src: &[Number]) {
1808 if let Some(dv) = v
1809 .as_any_mut()
1810 .downcast_mut::<pounce_linalg::dense_vector::DenseVector>()
1811 {
1812 dv.set_values(src);
1813 return;
1814 }
1815 if let Some(cv) = v
1816 .as_any_mut()
1817 .downcast_mut::<pounce_linalg::CompoundVector>()
1818 {
1819 let mut off = 0usize;
1820 for k in 0..cv.n_comps() {
1821 let blk = cv.comp_mut(k);
1822 let dim = blk.dim() as usize;
1823 let dblk = blk
1824 .as_any_mut()
1825 .downcast_mut::<pounce_linalg::dense_vector::DenseVector>()
1826 .expect("clamp_against_slack: CompoundVector blocks must be DenseVectors");
1827 dblk.set_values(&src[off..off + dim]);
1828 off += dim;
1829 }
1830 return;
1831 }
1832 panic!("clamp_against_slack: unsupported Vector kind");
1833}
1834
1835/// Per-element kappa-sigma clamp — the elementwise arithmetic at the
1836/// heart of `IpIpoptAlg.cpp:correct_bound_multiplier` (lines
1837/// 1090-1133). For each index `i`:
1838///
1839/// ```text
1840/// slack_i = max(slack_i, tiny_double) // avoid /0
1841/// z_lo_i = mu / (kappa_sigma * slack_i)
1842/// z_hi_i = kappa_sigma * mu / slack_i
1843/// z_i ← clamp(z_i, z_lo_i, z_hi_i)
1844/// ```
1845///
1846/// Returns the maximum elementwise correction magnitude (matching
1847/// upstream's `Max(max_correction_up, max_correction_low)`).
1848///
1849/// `kappa_sigma < 1` short-circuits to the identity per upstream's
1850/// guard at line 1065.
1851pub fn kappa_sigma_clamp(
1852 z: &mut [Number],
1853 slack: &[Number],
1854 mu: Number,
1855 kappa_sigma: Number,
1856) -> Number {
1857 debug_assert_eq!(z.len(), slack.len());
1858 if kappa_sigma < 1.0 {
1859 return 0.0;
1860 }
1861 let mut max_correction = 0.0_f64;
1862 for (zi, &si) in z.iter_mut().zip(slack.iter()) {
1863 let s_safe = si.max(Number::MIN_POSITIVE);
1864 let lo = mu / (kappa_sigma * s_safe);
1865 let hi = kappa_sigma * mu / s_safe;
1866 let clamped = zi.clamp(lo, hi);
1867 let delta = (clamped - *zi).abs();
1868 if delta > max_correction {
1869 max_correction = delta;
1870 }
1871 *zi = clamped;
1872 }
1873 max_correction
1874}
1875
1876#[cfg(test)]
1877mod tests {
1878 use super::*;
1879
1880 #[test]
1881 fn kappa_sigma_below_one_is_identity() {
1882 let mut z = vec![1.0, 2.0, 3.0];
1883 let slack = [1.0, 1.0, 1.0];
1884 let m = kappa_sigma_clamp(&mut z, &slack, 1.0, 0.5);
1885 assert_eq!(m, 0.0);
1886 assert_eq!(z, [1.0, 2.0, 3.0]);
1887 }
1888
1889 #[test]
1890 fn within_band_is_unchanged() {
1891 // mu=1, kappa=10, slack=1 → band [0.1, 10]. z=1 → unchanged.
1892 let mut z = vec![1.0];
1893 let slack = [1.0];
1894 let m = kappa_sigma_clamp(&mut z, &slack, 1.0, 10.0);
1895 assert_eq!(m, 0.0);
1896 assert_eq!(z, [1.0]);
1897 }
1898
1899 #[test]
1900 fn above_upper_clamped_down() {
1901 // mu=1, kappa=10, slack=1 → upper = 10. z=100 → 10.
1902 let mut z = vec![100.0];
1903 let slack = [1.0];
1904 let m = kappa_sigma_clamp(&mut z, &slack, 1.0, 10.0);
1905 assert!((m - 90.0).abs() < 1e-13);
1906 assert_eq!(z, [10.0]);
1907 }
1908
1909 #[test]
1910 fn below_lower_clamped_up() {
1911 // mu=1, kappa=10, slack=1 → lower = 0.1. z=0.001 → 0.1.
1912 let mut z = vec![0.001];
1913 let slack = [1.0];
1914 let m = kappa_sigma_clamp(&mut z, &slack, 1.0, 10.0);
1915 assert!((m - 0.099).abs() < 1e-13);
1916 assert!((z[0] - 0.1).abs() < 1e-15);
1917 }
1918
1919 #[test]
1920 fn returns_max_over_components() {
1921 let mut z = vec![100.0, 0.001];
1922 let slack = [1.0, 1.0];
1923 let m = kappa_sigma_clamp(&mut z, &slack, 1.0, 10.0);
1924 assert!((m - 90.0).abs() < 1e-13);
1925 assert_eq!(z[0], 10.0);
1926 assert!((z[1] - 0.1).abs() < 1e-15);
1927 }
1928
1929 #[test]
1930 fn slack_clamped_to_min_positive_avoids_division_by_zero() {
1931 let mut z = vec![1e100];
1932 let slack = [0.0];
1933 let _ = kappa_sigma_clamp(&mut z, &slack, 1.0, 10.0);
1934 assert!(z[0].is_finite() || z[0] == 1e100);
1935 }
1936
1937 /// The restoration slot is exercised structurally:
1938 /// `IpoptAlgorithm::with_restoration` accepts a
1939 /// `Box<dyn RestorationPhase>` and the trait's default
1940 /// `perform_restoration` returns `Failed`. End-to-end coverage
1941 /// (iterate() → line-search-Failed → restoration → recovered)
1942 /// lands in the Phase 9 integration suite alongside the nested
1943 /// IPM driver.
1944 struct _DummyResto;
1945 impl RestorationPhase for _DummyResto {}
1946}