Skip to main content

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}