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