Skip to main content

pounce_algorithm/
application.rs

1//! User-facing application object — port of `Interfaces/IpIpoptApplication.{hpp,cpp}`.
2//!
3//! # Crate placement
4//!
5//! `IpoptApplication` lives in `pounce-algorithm` (rather than
6//! alongside the other Interfaces-side ports in `pounce-nlp`) because
7//! `optimize_tnlp` needs to drive the full IPM: it constructs a
8//! `TNLPAdapter` + `OrigIpoptNlp` (from `pounce-nlp`) and hands the
9//! NLP off to an [`IpoptAlgorithm`] (this crate). `pounce-nlp` cannot
10//! depend on `pounce-algorithm` (the reverse already exists), so
11//! orchestration must live on the algorithm side. Public callers
12//! continue to import via `pounce_algorithm::IpoptApplication`.
13//!
14//! `optimize_tnlp` routes every problem — constrained or not —
15//! through the same primal-dual IPM, exactly as upstream Ipopt does:
16//! it builds the algorithm via [`crate::alg_builder::AlgorithmBuilder`]
17//! (default backend MA57 from `pounce-hsl`) and runs
18//! [`IpoptAlgorithm::optimize`].
19
20use crate::alg_builder::{
21    AlgorithmBuilder, HessianApproxChoice, LineSearchChoice, LinearBackendFactory,
22    LinearSolverChoice, MuStrategyChoice,
23};
24use crate::hess::lim_mem_quasi_newton::UpdateType;
25use crate::ipopt_alg::IpoptAlgorithm;
26use crate::ipopt_cq::IpoptCalculatedQuantities;
27use crate::ipopt_data::IpoptData as AlgIpoptData;
28use crate::ipopt_nlp::IpoptNlp;
29use crate::iterates_vector::IteratesVector;
30use crate::restoration::RestorationPhase;
31use crate::upstream_options::register_all_upstream_options;
32
33/// Factory that constructs a fresh restoration-phase strategy on
34/// demand. The outer algorithm owns at most one restoration object,
35/// so the factory is invoked once per `optimize_tnlp` call. The
36/// factory is `FnMut` to allow callers to capture a builder that
37/// internally reuses caches across builds.
38pub type RestorationFactory = Box<dyn FnMut() -> Box<dyn RestorationPhase>>;
39
40/// Provider that mints fresh [`RestorationFactory`] instances on
41/// demand. Used by drivers that need to run the inner IPM more than
42/// once per `optimize_tnlp` call — notably the Phase-3 ℓ₁-exact
43/// penalty-barrier outer loop (pounce#10), which the existing
44/// `RestorationFactory` cannot support because pounce's default
45/// `make_default_restoration_factory` is a one-shot. Callers wire
46/// this via [`IpoptApplication::set_restoration_factory_provider`].
47pub type RestorationFactoryProvider = Box<dyn FnMut() -> RestorationFactory>;
48
49/// Callback fired by [`IpoptApplication::optimize_constrained`] once
50/// the IPM has converged (status `SolveSucceeded` or
51/// `SolvedToAcceptableLevel`) and before the user TNLP's
52/// `finalize_solution` runs. Receives borrowed handles into the
53/// algorithm's converged state.
54///
55/// **Use case**: post-optimal sensitivity analysis (pounce#7 /
56/// `pounce-sensitivity`). The callback receives a shared handle to
57/// the PD solver so a `SensBacksolver` adapter can run backsolves
58/// against the converged KKT factor — and so that handle may outlive
59/// the call frame (e.g. the public `Solver` session API retains the
60/// factor for repeated `parametric_step` / `kkt_solve` calls);
61/// receives the data / cq / nlp handles so the adapter can reproduce
62/// the augmented-system coefficient layout the IPM converged at.
63///
64/// **Not** the same as `set_intermediate_callback` (per-iteration
65/// progress notification) — this fires exactly once per `optimize_*`
66/// call, only on success.
67pub type ConvergedCallback = Box<
68    dyn FnMut(
69        &crate::ipopt_data::IpoptDataHandle,
70        &crate::ipopt_cq::IpoptCqHandle,
71        &Rc<RefCell<dyn pounce_nlp::ipopt_nlp::IpoptNlp>>,
72        Rc<RefCell<crate::kkt::pd_full_space_solver::PdFullSpaceSolver>>,
73    ),
74>;
75use pounce_common::diagnostics::DiagnosticsState;
76use pounce_common::exception::{ExceptionKind, SolverException};
77use pounce_common::journalist::{JournalLevel, Journalist};
78use pounce_common::options_list::OptionsList;
79use pounce_common::reg_options::{PrintOptionsMode, RegisteredOptions};
80use pounce_common::timing::TimingStatistics;
81use pounce_common::types::{Index, Number};
82use pounce_linalg::dense_vector::DenseVectorSpace;
83use pounce_linsol::summary::LinearSolverSummary;
84use pounce_linsol::SparseSymLinearSolverInterface;
85use pounce_nlp::alg_types::SolverReturn;
86use pounce_nlp::orig_ipopt_nlp::{ConstObjScaling, OrigIpoptNlp, ScalingMethod};
87use pounce_nlp::return_codes::ApplicationReturnStatus;
88use pounce_nlp::solve_statistics::SolveStatistics;
89use pounce_nlp::tnlp::{
90    IpoptCq as TnlpIpoptCq, IpoptData as TnlpIpoptData, NlpInfo, Solution, TNLP,
91};
92use pounce_nlp::tnlp_adapter::{
93    FixedVarTreatment, TNLPAdapter, DEFAULT_NLP_LOWER_BOUND_INF, DEFAULT_NLP_UPPER_BOUND_INF,
94};
95use std::cell::RefCell;
96use std::fmt;
97use std::path::Path;
98use std::rc::Rc;
99use std::sync::{Arc, Mutex};
100use std::time::Instant;
101
102pub struct IpoptApplication {
103    options: OptionsList,
104    reg_options: Rc<RegisteredOptions>,
105    journalist: Rc<Journalist>,
106    statistics: RefCell<SolveStatistics>,
107    /// Shared per-subsystem timing accumulator. Re-created at the top of
108    /// every solve (so back-to-back `optimize_tnlp` calls don't bleed
109    /// timings across invocations) and handed to the data, the NLP, and
110    /// any other consumer via `Rc`. Reported by [`Self::timing_stats`]
111    /// after the solve completes.
112    timing: RefCell<Rc<TimingStatistics>>,
113    /// Optional override factory for the symmetric linear-solver
114    /// backend. When `None`, we ship the workspace default (MA57 via
115    /// `pounce-hsl`). Tests can plug a stub via [`Self::set_linear_backend_factory`].
116    linear_backend_factory: Option<LinearBackendFactory>,
117    /// Optional factory for the restoration phase. Lives outside this
118    /// crate because `pounce-algorithm` cannot depend on
119    /// `pounce-restoration` (the dep edge is the other way). Callers
120    /// that need restoration plug a factory via
121    /// [`Self::set_restoration_factory`]; when unset, the outer
122    /// algorithm runs without a restoration fallback and surfaces
123    /// `RestorationFailure` as soon as the line-search would otherwise
124    /// jump into restoration.
125    restoration_factory: Option<RestorationFactory>,
126    /// Shared diagnostic-dump state, installed by the CLI when the
127    /// user passes `--dump <cat>:<spec>`. When set, the application
128    /// propagates an `Rc<DiagnosticsState>` into [`IpoptAlgorithm`]
129    /// via [`IpoptAlgorithm::with_diagnostics`] so the KKT solver and
130    /// other dump sites can consult per-iter gating.
131    diagnostics: Option<Rc<DiagnosticsState>>,
132    /// Optional interactive debugger hook. When set, it is moved into
133    /// the main [`IpoptAlgorithm`] for the next `optimize_*` call via
134    /// [`IpoptAlgorithm::with_debug_hook`], so a REPL or agent can pause
135    /// at each iteration to inspect / mutate live state. Consumed on use
136    /// (one solve per installed hook).
137    debug_hook: Option<std::rc::Rc<std::cell::RefCell<dyn crate::debug::DebugHook>>>,
138    /// Provider for the BNW outer loop (pounce#10 Phase 3). When set,
139    /// `optimize_constrained` consults the provider before each inner
140    /// solve, replacing `restoration_factory` with a fresh one so
141    /// multi-pass drivers can run the inner IPM repeatedly without
142    /// tripping the default factory's one-shot guard.
143    restoration_factory_provider: Option<RestorationFactoryProvider>,
144    /// Optional hook fired once per `optimize_*` call on convergence,
145    /// before the user TNLP's `finalize_solution`. See
146    /// [`ConvergedCallback`].
147    on_converged: Option<ConvergedCallback>,
148    /// When `true`, the per-iteration `IterRecord` trajectory is
149    /// captured into [`SolveStatistics::iterations`] for downstream
150    /// consumers (the JSON solve report in pounce-cli, pounce#8). Off
151    /// by default so library callers that never read the iterations
152    /// vector don't pay the per-iter alloc.
153    record_iter_history: bool,
154    /// Shared sink that the linear-solver backend writes a rolling
155    /// [`LinearSolverSummary`] into after every factor. Reset at the
156    /// top of every solve (so back-to-back `optimize_tnlp` calls don't
157    /// bleed stats across invocations) and read out via
158    /// [`Self::linear_solver_summary`] once the solve returns. Only
159    /// the workspace-default FERAL backend (via
160    /// [`default_backend_factory_with_sink`]) wires the sink today;
161    /// custom factories plugged through [`Self::set_linear_backend_factory`]
162    /// and the HSL MA57 backend leave the sink empty.
163    linsol_summary_sink: Arc<Mutex<LinearSolverSummary>>,
164    /// Phase 5c (§6) SQP warm-start input. When `Some`, the next
165    /// `optimize_tnlp` call on the SQP path consumes the iterate
166    /// instead of cold-starting; consumed once per solve, then
167    /// auto-cleared. The IPM path ignores this field. Wire-set
168    /// via [`Self::set_sqp_warm_start`].
169    sqp_warm_start: Option<crate::sqp::SqpIterates>,
170    /// Phase 5c (§6) SQP warm-start output. Populated by every
171    /// `optimize_sqp_tnlp` call with the final QP working set.
172    /// Stays valid until the next solve (which overwrites it).
173    /// Accessed via [`Self::last_sqp_working_set`].
174    sqp_last_working_set: Option<pounce_qp::WorkingSet>,
175    /// Full primal-dual warm-start iterate for the IPM path, captured by
176    /// the interactive debugger's `resolve` command. When `Some`, the
177    /// next `optimize_tnlp` installs this 8-vector (algorithm space)
178    /// directly onto `data.curr` before the iterate initializer runs, so
179    /// a warm `resolve` continues from the paused interior point rather
180    /// than cold-restarting the duals. Consumed once per solve, then
181    /// auto-cleared. Requires `warm_start_init_point=yes` so the
182    /// re-optimize branch of `WarmStartIterateInitializer` keeps the
183    /// installed iterate. Wire-set via [`Self::set_warm_start_iterate`].
184    warm_start_iterate: Option<crate::debug::IterateSnapshot>,
185}
186
187impl fmt::Debug for IpoptApplication {
188    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
189        f.debug_struct("IpoptApplication")
190            .field("options", &self.options)
191            .field("statistics", &self.statistics)
192            .finish_non_exhaustive()
193    }
194}
195
196impl Default for IpoptApplication {
197    fn default() -> Self {
198        Self::new()
199    }
200}
201
202impl IpoptApplication {
203    /// New application with empty options and a default journalist.
204    /// Equivalent to `IpoptApplication::IpoptApplication(true,true)`.
205    pub fn new() -> Self {
206        let reg = RegisteredOptions::default();
207        // Registration of a fresh registry can only fail on a duplicate
208        // name, which would be a programming error in `reg_op`.
209        register_all_upstream_options(&reg)
210            .unwrap_or_else(|e| panic!("Upstream options registration failed: {e}"));
211        pounce_presolve::register_options(&reg)
212            .unwrap_or_else(|e| panic!("Presolve options registration failed: {e}"));
213        let reg = Rc::new(reg);
214        Self {
215            options: OptionsList::with_registered(Rc::clone(&reg)),
216            reg_options: reg,
217            journalist: Rc::new(Journalist::new()),
218            statistics: RefCell::new(SolveStatistics::new()),
219            timing: RefCell::new(Rc::new(TimingStatistics::new())),
220            linear_backend_factory: None,
221            restoration_factory: None,
222            diagnostics: None,
223            debug_hook: None,
224            restoration_factory_provider: None,
225            on_converged: None,
226            record_iter_history: false,
227            linsol_summary_sink: Arc::new(Mutex::new(LinearSolverSummary::default())),
228            sqp_warm_start: None,
229            sqp_last_working_set: None,
230            warm_start_iterate: None,
231        }
232    }
233
234    pub fn options(&self) -> &OptionsList {
235        &self.options
236    }
237
238    pub fn options_mut(&mut self) -> &mut OptionsList {
239        &mut self.options
240    }
241
242    pub fn registered_options(&self) -> &Rc<RegisteredOptions> {
243        &self.reg_options
244    }
245
246    pub fn journalist(&self) -> &Rc<Journalist> {
247        &self.journalist
248    }
249
250    /// Plug a custom symmetric-linear-solver factory. Useful for tests
251    /// that want to swap MA57 for a stub. Production callers should
252    /// leave this unset — the default ([`default_backend_factory`])
253    /// returns the workspace's MA57 binding.
254    pub fn set_linear_backend_factory(&mut self, factory: LinearBackendFactory) {
255        self.linear_backend_factory = Some(factory);
256    }
257
258    /// Plug a restoration-phase factory. Called once per
259    /// `optimize_tnlp` invocation to mint a fresh
260    /// `Box<dyn RestorationPhase>` that the outer algorithm uses as
261    /// its line-search restoration fallback. Lives behind a setter
262    /// (rather than at construction) because the concrete restoration
263    /// strategies live in `pounce-restoration`, which depends on this
264    /// crate; consumers in `pounce-cli` / integration tests wire the
265    /// factory at the application boundary.
266    pub fn set_restoration_factory(&mut self, factory: RestorationFactory) {
267        self.restoration_factory = Some(factory);
268    }
269
270    /// Install the shared diagnostics state. Once set, every
271    /// subsequent `optimize_tnlp` call forwards the state into the
272    /// algorithm via [`IpoptAlgorithm::with_diagnostics`] so the KKT
273    /// solver can emit `--dump kkt:...` artifacts.
274    pub fn set_diagnostics(&mut self, diag: Rc<DiagnosticsState>) {
275        self.diagnostics = Some(diag);
276    }
277
278    /// Install an interactive debugger hook for the next `optimize_*`
279    /// call. The hook is moved into the main [`IpoptAlgorithm`] and
280    /// consumed by that solve; reinstall it to debug a subsequent solve.
281    pub fn set_debug_hook(
282        &mut self,
283        hook: std::rc::Rc<std::cell::RefCell<dyn crate::debug::DebugHook>>,
284    ) {
285        self.debug_hook = Some(hook);
286    }
287
288    /// Read-side accessor for the installed diagnostics state, if any.
289    /// Lets the CLI write the top-level manifest/timing files after
290    /// the solve completes.
291    pub fn diagnostics(&self) -> Option<Rc<DiagnosticsState>> {
292        self.diagnostics.as_ref().map(Rc::clone)
293    }
294
295    /// Plug a restoration-phase **factory provider** for drivers that
296    /// need to run the inner IPM more than once per `optimize_tnlp`
297    /// call (notably the Phase-3 ℓ₁-exact penalty-barrier outer loop,
298    /// pounce#10). On each inner solve, the application consults the
299    /// provider to mint a fresh [`RestorationFactory`], replacing any
300    /// stale one, so the default one-shot restoration factory does
301    /// not panic on its second invocation. If both `set_restoration_factory`
302    /// and this are configured, the provider wins.
303    pub fn set_restoration_factory_provider(&mut self, provider: RestorationFactoryProvider) {
304        self.restoration_factory_provider = Some(provider);
305    }
306
307    /// Register a callback to run once the IPM has converged (status
308    /// [`ApplicationReturnStatus::SolveSucceeded`] or
309    /// [`ApplicationReturnStatus::SolvedToAcceptableLevel`]) but before
310    /// `finalize_solution` flows back to the TNLP. See
311    /// [`ConvergedCallback`] for the use case (post-optimal sensitivity).
312    pub fn set_on_converged(&mut self, cb: ConvergedCallback) {
313        self.on_converged = Some(cb);
314    }
315
316    /// Enable per-iteration trajectory capture. After the solve
317    /// returns, [`Self::statistics()`] exposes
318    /// [`pounce_nlp::solve_statistics::SolveStatistics::iterations`]
319    /// populated with one [`pounce_nlp::solve_statistics::IterRecord`]
320    /// per accepted iterate. Off by default — the `pounce_sens` and
321    /// `pounce` binaries opt in when `--json-output` is passed.
322    pub fn enable_iter_history(&mut self) {
323        self.record_iter_history = true;
324    }
325
326    /// Read an `ipopt.opt`-format options file. Equivalent to
327    /// `IpoptApplication::Initialize(const std::string& options_file)`.
328    pub fn initialize_with_options_file(&mut self, path: &Path) -> Result<(), SolverException> {
329        let txt = std::fs::read_to_string(path).map_err(|e| {
330            SolverException::new(
331                ExceptionKind::IPOPT_APPLICATION_ERROR,
332                format!("could not read options file {}: {}", path.display(), e),
333                file!(),
334                line!() as Index,
335            )
336        })?;
337        self.options.read_from_str(&txt, true)?;
338        self.open_output_file_journal();
339        Ok(())
340    }
341
342    /// Read options from a string in `ipopt.opt` format. Useful for
343    /// tests and embedded callers.
344    pub fn initialize_with_options_str(&mut self, s: &str) -> Result<(), SolverException> {
345        self.options.read_from_str(s, true)?;
346        self.open_output_file_journal();
347        Ok(())
348    }
349
350    /// Honor `output_file` / `file_print_level` / `file_append`: when
351    /// `output_file` is non-empty, attach a `FileJournal` named
352    /// `"OutputFile:<fname>"` at the requested level. Mirrors
353    /// `IpoptApplication::OpenOutputFile` (called from `Initialize`).
354    /// No-op if `output_file` is unset, empty, or could not be opened.
355    ///
356    /// NOTE: pounce's iteration output currently bypasses the
357    /// journalist and writes directly to stdout. The file journal is
358    /// attached and the timing report (gated by `print_timing_statistics`)
359    /// is mirrored to it; per-iter rows will start landing in the file
360    /// once the iter-output path is routed through the journalist.
361    fn open_output_file_journal(&self) {
362        let fname = match self.options.get_string_value("output_file", "") {
363            Ok((v, true)) if !v.is_empty() => v,
364            _ => return,
365        };
366        let level_int = self
367            .options
368            .get_integer_value("file_print_level", "")
369            .ok()
370            .and_then(|(v, f)| f.then_some(v))
371            .unwrap_or(5);
372        let level = journal_level_from_int(level_int);
373        let append = self
374            .options
375            .get_bool_value("file_append", "")
376            .ok()
377            .and_then(|(v, f)| f.then_some(v))
378            .unwrap_or(false);
379        let jname = format!("OutputFile:{}", fname);
380        let _ = self
381            .journalist
382            .add_file_journal(&jname, &fname, level, append);
383    }
384
385    /// No-op initialize (just succeeds). Mirrors
386    /// `IpoptApplication::Initialize(bool allow_clobber)` with no
387    /// options file.
388    pub fn initialize(&mut self) -> Result<(), SolverException> {
389        Ok(())
390    }
391
392    /// Mirror `IpoptApplication::OpenOutputFile`. Sets the `output_file`
393    /// / `file_print_level` options and attaches a matching
394    /// `FileJournal` named `OutputFile:<fname>` to the journalist.
395    /// Returns `false` if the file could not be opened or the option
396    /// store rejected the request (e.g. clamped print level).
397    pub fn open_output_file(&mut self, fname: &str, print_level: i32) -> bool {
398        if self
399            .options
400            .set_string_value("output_file", fname, true, false)
401            .is_err()
402        {
403            return false;
404        }
405        if self
406            .options
407            .set_integer_value("file_print_level", print_level as Index, true, false)
408            .is_err()
409        {
410            return false;
411        }
412        let level = journal_level_from_int(print_level);
413        let jname = format!("OutputFile:{}", fname);
414        // Drop any previous file journal so a second call switches files
415        // cleanly. `add_file_journal` would otherwise refuse to attach
416        // a duplicate by name; remove-by-name isn't in the journalist
417        // API, so we settle for the name-collision case here.
418        self.journalist
419            .add_file_journal(&jname, fname, level, false)
420            .is_some()
421    }
422
423    /// Wrap a TNLP and report problem dimensions. Used in tests until
424    /// the full IPM path covers every entry shape.
425    pub fn problem_dimensions(&self, tnlp: &mut dyn TNLP) -> Option<NlpInfo> {
426        tnlp.get_nlp_info()
427    }
428
429    pub fn statistics(&self) -> SolveStatistics {
430        self.statistics.borrow().clone()
431    }
432
433    /// Shared timing accumulator from the most recent `optimize_tnlp`
434    /// call. Each subsystem (algorithm, NLP, KKT solver) bumped its own
435    /// fields during the solve; consumers read totals out of the
436    /// returned `Rc`. The instance is replaced at the top of every
437    /// subsequent solve, so cloning the `Rc` and holding it past a
438    /// re-solve will give you the previous solve's timings — by design.
439    pub fn timing_stats(&self) -> Rc<TimingStatistics> {
440        Rc::clone(&self.timing.borrow())
441    }
442
443    /// Aggregate linear-solver post-mortem from the most recent
444    /// `optimize_tnlp` call. `Some` when the workspace-default FERAL
445    /// backend ran at least one factor; `None` when no factors were
446    /// recorded (custom factory plugged via
447    /// [`Self::set_linear_backend_factory`], or solve aborted before
448    /// the first KKT factor). Reset at the top of every solve.
449    pub fn linear_solver_summary(&self) -> Option<LinearSolverSummary> {
450        let guard = self.linsol_summary_sink.lock().ok()?;
451        if guard.is_empty() {
452            None
453        } else {
454            Some(guard.clone())
455        }
456    }
457
458    /// Drive a solve.
459    ///
460    /// * Constrained problems (`m > 0`) take the primal-dual IPM path:
461    ///   build a `TNLPAdapter` → `OrigIpoptNlp`, run the
462    ///   [`AlgorithmBuilder`] with the workspace MA57 backend, and
463    ///   call [`IpoptAlgorithm::optimize`]. The `SolverReturn` →
464    ///   `ApplicationReturnStatus` mapping mirrors the table in
465    ///   `ref/Ipopt/AGENT_REFERENCE/MAIN_LOOP.md` ("exception →
466    ///   SolverReturn map").
467    /// * Unconstrained problems (`m == 0`) keep going through the
468    ///   in-`pounce-nlp` Newton driver so the trivial path is
469    ///   independent of the linear-solver backend.
470    pub fn optimize_tnlp(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
471        // Top-level algorithm dispatch (Phase 5b §7.1). When the
472        // `algorithm` option resolves to "active-set-sqp", route
473        // to the Phase 5b SQP path; otherwise fall through to the
474        // existing IPM flow unchanged.
475        if self.is_sqp_algorithm_selected() {
476            return self.optimize_sqp_tnlp(tnlp);
477        }
478        let info = match tnlp.borrow_mut().get_nlp_info() {
479            Some(info) => info,
480            None => return ApplicationReturnStatus::InvalidProblemDefinition,
481        };
482        // ℓ₁-exact penalty-barrier opt-in (pounce#10).
483        // Phase 3 wraps the user TNLP and runs an outer Byrd-Nocedal-
484        // Waltz ρ-escalation loop around the constrained IPM, with a
485        // honest-infeasibility status upgrade when the slacks fail to
486        // collapse at saturated ρ. Phase-1/2 one-shot use is preserved
487        // when `l1_penalty_max_outer_iter == 1`. The wrapper is a
488        // no-op for problems with no equality rows, so the
489        // unconstrained dispatch below is unaffected when there is
490        // nothing to wrap.
491        if info.m > 0 && self.is_l1_penalty_enabled() {
492            if let Some(status) = self.run_l1_penalty_outer_loop(Rc::clone(&tnlp)) {
493                return status;
494            }
495            // Falls through: wrapper construction failed (inner refused
496            // get_nlp_info / get_bounds_info) or no equality rows to
497            // slack. Standard dispatch runs unmodified.
498        }
499        // Phase 3.5 auto-fallback (pounce#10): if the standard solve
500        // ends in a trigger-class status, retry transparently with
501        // the wrapper. Promote the retry's status only if it returns
502        // SolveSucceeded — otherwise return the original. Skipped if
503        // the user already opted into the wrapper above (this avoids
504        // a double pass and keeps semantics predictable).
505        if info.m > 0 && self.is_l1_fallback_enabled() && !self.is_l1_penalty_enabled() {
506            return self.run_with_l1_fallback(tnlp);
507        }
508        // μ-strategy auto-fallback (pounce#138): if the standard solve
509        // only reaches Solved_To_Acceptable_Level, retry once with the
510        // opposite mu_strategy and promote only on Solve_Succeeded.
511        // Applies to constrained and unconstrained alike (both run the
512        // same IPM). Independent of, and lower priority than, the ℓ₁
513        // fallback above.
514        if self.is_mu_strategy_fallback_enabled() {
515            return self.run_with_mu_strategy_fallback(tnlp);
516        }
517        // Every problem — constrained or not — goes through the same
518        // primal-dual IPM, exactly as upstream Ipopt does. There is no
519        // separate "unconstrained Newton" path: the linear-solver
520        // backend (FERAL/MA57) handles the augmented system, so the
521        // sparse IPM covers `m == 0` at any `n` without a dense-Hessian
522        // blowup.
523        self.optimize_constrained(tnlp)
524    }
525
526    /// Read the ℓ₁ wrapper master switch from the OptionsList.
527    /// Default `false` when the option is not set.
528    fn is_l1_penalty_enabled(&self) -> bool {
529        self.options
530            .get_bool_value("l1_exact_penalty_barrier", "")
531            .ok()
532            .and_then(|(v, found)| found.then_some(v))
533            .unwrap_or(false)
534    }
535
536    fn l1_penalty_init(&self) -> Number {
537        self.options
538            .get_numeric_value("l1_penalty_init", "")
539            .ok()
540            .and_then(|(v, found)| found.then_some(v))
541            .unwrap_or(1.0)
542    }
543    fn l1_penalty_max(&self) -> Number {
544        self.options
545            .get_numeric_value("l1_penalty_max", "")
546            .ok()
547            .and_then(|(v, found)| found.then_some(v))
548            .unwrap_or(1.0e6)
549    }
550    fn l1_penalty_increase_factor(&self) -> Number {
551        self.options
552            .get_numeric_value("l1_penalty_increase_factor", "")
553            .ok()
554            .and_then(|(v, found)| found.then_some(v))
555            .unwrap_or(8.0)
556    }
557    fn l1_penalty_max_outer_iter(&self) -> usize {
558        self.options
559            .get_integer_value("l1_penalty_max_outer_iter", "")
560            .ok()
561            .and_then(|(v, found)| found.then_some(v))
562            .unwrap_or(8) as usize
563    }
564    fn l1_slack_tol(&self) -> Number {
565        self.options
566            .get_numeric_value("l1_slack_tol", "")
567            .ok()
568            .and_then(|(v, found)| found.then_some(v))
569            .unwrap_or(1.0e-6)
570    }
571    fn l1_steering_factor(&self) -> Number {
572        self.options
573            .get_numeric_value("l1_steering_factor", "")
574            .ok()
575            .and_then(|(v, found)| found.then_some(v))
576            .unwrap_or(10.0)
577    }
578    fn is_l1_fallback_enabled(&self) -> bool {
579        self.options
580            .get_bool_value("l1_fallback_on_restoration_failure", "")
581            .ok()
582            .and_then(|(v, found)| found.then_some(v))
583            .unwrap_or(false)
584    }
585
586    /// Read the μ-strategy auto-fallback switch (pounce#138).
587    /// Default `false` when the option is not set.
588    fn is_mu_strategy_fallback_enabled(&self) -> bool {
589        self.options
590            .get_bool_value("mu_strategy_fallback", "")
591            .ok()
592            .and_then(|(v, found)| found.then_some(v))
593            .unwrap_or(false)
594    }
595
596    /// Has the user set `algorithm = active-set-sqp`? Reads the
597    /// string option and matches case-insensitively against the
598    /// design-note §7.1 spelling. Any value other than
599    /// "active-set-sqp" (including absence) routes to the
600    /// default IPM path.
601    /// Stash a warm-start iterate for the SQP path. Consumed by
602    /// the next `optimize_tnlp` call when the `algorithm` option
603    /// resolves to `active-set-sqp`; the IPM path ignores it.
604    /// Phase 5c (§6) — the parametric / MPC warm-start hand-off.
605    ///
606    /// The iterate is auto-cleared after use, so a follow-up
607    /// solve without an intervening `set_sqp_warm_start` call
608    /// cold-starts.
609    pub fn set_sqp_warm_start(&mut self, warm: crate::sqp::SqpIterates) {
610        self.sqp_warm_start = Some(warm);
611    }
612
613    /// Drop any pending warm-start iterate without solving.
614    pub fn clear_sqp_warm_start(&mut self) {
615        self.sqp_warm_start = None;
616    }
617
618    /// Install a full primal-dual warm-start iterate for the next IPM
619    /// `optimize_tnlp`. Captured by the debugger's `resolve` so the
620    /// re-solve continues from the paused interior point. The caller is
621    /// responsible for also enabling `warm_start_init_point=yes` (and
622    /// usually `warm_start_target_mu=<μ>`) so the re-optimize branch of
623    /// `WarmStartIterateInitializer` preserves the installed iterate.
624    /// Consumed once per solve, then auto-cleared.
625    pub fn set_warm_start_iterate(&mut self, snap: crate::debug::IterateSnapshot) {
626        self.warm_start_iterate = Some(snap);
627    }
628
629    /// Return the final QP working set from the most recent SQP
630    /// solve, or `None` if the last solve wasn't SQP, didn't
631    /// produce a working set (cold-start declared the iterate
632    /// optimal before solving any QP), or no SQP solve has run.
633    pub fn last_sqp_working_set(&self) -> Option<&pounce_qp::WorkingSet> {
634        self.sqp_last_working_set.as_ref()
635    }
636
637    fn is_sqp_algorithm_selected(&self) -> bool {
638        match self.options.get_string_value("algorithm", "") {
639            Ok((v, true)) => v.eq_ignore_ascii_case("active-set-sqp"),
640            _ => false,
641        }
642    }
643
644    /// Phase 5b SQP entry point. Builds the same NLP chain
645    /// (`TNLPAdapter` → `OrigIpoptNlp` → `IpoptNlpAdapter`) the
646    /// IPM uses, then runs `SqpAlgorithm::optimize`. Maps the
647    /// `SqpResult.status` back to `ApplicationReturnStatus` and
648    /// hands the final iterate to the user TNLP's
649    /// `finalize_solution` callback via `finalize_via_sqp`.
650    fn optimize_sqp_tnlp(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
651        use pounce_nlp::orig_ipopt_nlp::OrigIpoptNlp;
652        use pounce_nlp::tnlp_adapter::TNLPAdapter;
653        use pounce_nlp::ConstObjScaling;
654
655        let adapter = match TNLPAdapter::new(Rc::clone(&tnlp)) {
656            Ok(a) => Rc::new(RefCell::new(a)),
657            Err(_) => return ApplicationReturnStatus::InvalidProblemDefinition,
658        };
659        // The SQP path never runs gradient-based scaling, but the
660        // constant `obj_scaling_factor` (negative ⇒ maximize) still
661        // applies via the OrigIpoptNlp constructor.
662        let obj_scaling_factor = self
663            .options
664            .get_numeric_value("obj_scaling_factor", "")
665            .ok()
666            .and_then(|(v, f)| f.then_some(v))
667            .unwrap_or(1.0);
668        let orig_nlp = match OrigIpoptNlp::new(
669            Rc::clone(&adapter),
670            Rc::new(ConstObjScaling(obj_scaling_factor)),
671        ) {
672            Ok(n) => n,
673            Err(_) => return ApplicationReturnStatus::InternalError,
674        };
675        let nlp_rc: Rc<RefCell<dyn IpoptNlp>> = Rc::new(RefCell::new(orig_nlp));
676
677        let mut sqp_adapter = crate::sqp::IpoptNlpAdapter::new(Rc::clone(&nlp_rc));
678
679        let mut builder = self.algorithm_builder_snapshot();
680        builder.algorithm = crate::alg_builder::AlgorithmChoice::ActiveSetSqp;
681        let factory = self.make_backend_factory();
682        let mut alg = match builder.build_sqp_with_backend(factory) {
683            Some(a) => a,
684            None => return ApplicationReturnStatus::InternalError,
685        };
686
687        // Phase 5c (§6): consume any stashed warm-start iterate.
688        // `optimize_with_warm_start(warm=None)` is equivalent to
689        // `optimize`, so cold callers see no change.
690        let warm = self.sqp_warm_start.take();
691        let res = match alg.optimize_with_warm_start(&mut sqp_adapter, warm) {
692            Ok(r) => r,
693            Err(e) => {
694                if std::env::var_os("POUNCE_DBG_SQP").is_some() {
695                    tracing::warn!(target: "pounce::sqp", "[SQP] optimize_with_warm_start error: {e:?}");
696                }
697                return ApplicationReturnStatus::InternalError;
698            }
699        };
700        // Stash the result's working set so the next solve in a
701        // sequence can fetch it via `last_sqp_working_set`.
702        self.sqp_last_working_set = res.working_set.clone();
703        // Populate the shared `SolveStatistics` so the Python /
704        // C-API post-solve accessors (`GetIpoptIterCount`,
705        // `info["iter_count"]`, etc.) report the SQP outer-iter
706        // count rather than zero. Constraint-violation /
707        // dual-infeasibility residuals get the SQP-side values
708        // too. The IPM path overwrites this dict on its own
709        // solves, so SQP-vs-IPM mixing across solves stays
710        // honest.
711        {
712            let mut stats = self.statistics.borrow_mut();
713            stats.iteration_count = res.n_iter as Index;
714            stats.final_objective = res.obj;
715            stats.final_dual_inf = res.final_stationarity;
716            stats.final_constr_viol = res.final_constr_viol;
717            stats.final_compl = 0.0; // SQP has no barrier — no compl term.
718                                     // Unscaled residuals (pounce#173). The SQP path does not thread
719                                     // the nlp_scaling factors through to its residuals yet, so these
720                                     // mirror the SQP-side values: correct when no scaling is active
721                                     // (the common case) and a conservative proxy otherwise. Populated
722                                     // here so the info dict's `final_unscaled_*` keys are honest
723                                     // rather than left at the 0.0 default.
724            stats.final_unscaled_dual_inf = res.final_stationarity;
725            stats.final_unscaled_constr_viol = res.final_constr_viol;
726            stats.final_unscaled_compl = 0.0;
727            stats.final_unscaled_kkt_error = res.final_stationarity.max(res.final_constr_viol);
728        }
729        let (app_status, solver_status) = match res.status {
730            crate::sqp::SqpStatus::Optimal => (
731                ApplicationReturnStatus::SolveSucceeded,
732                pounce_nlp::SolverReturn::Success,
733            ),
734            crate::sqp::SqpStatus::MaxIter => (
735                ApplicationReturnStatus::MaximumIterationsExceeded,
736                pounce_nlp::SolverReturn::MaxiterExceeded,
737            ),
738            crate::sqp::SqpStatus::InfeasibleSubproblem => (
739                ApplicationReturnStatus::InfeasibleProblemDetected,
740                pounce_nlp::SolverReturn::LocalInfeasibility,
741            ),
742            crate::sqp::SqpStatus::LineSearchFailed => (
743                ApplicationReturnStatus::SearchDirectionBecomesTooSmall,
744                pounce_nlp::SolverReturn::ErrorInStepComputation,
745            ),
746        };
747
748        // Forward to the user TNLP's finalize_solution. We pass
749        // the SQP iterate and recovered multipliers via the
750        // OrigIpoptNlp's lifting hooks. Failure here is silent
751        // (we still return the algorithm's status) — the user
752        // sees the right ApplicationReturnStatus regardless.
753        let _ = finalize_via_sqp(&nlp_rc, &res, solver_status, &tnlp);
754
755        // Honor the opt-in status-fidelity gate on the SQP path too
756        // (pounce#173).
757        self.apply_kkt_fidelity_gate(app_status)
758    }
759
760    /// Opt-in status-fidelity gate (pounce#173), shared by the IPM and
761    /// SQP solve paths. When the user sets a positive `kkt_fidelity_tol`,
762    /// a reported `Solve_Succeeded` whose max-norm UNSCALED KKT error
763    /// (`SolveStatistics::final_unscaled_kkt_error`) exceeds it is
764    /// downgraded to `Solved_To_Acceptable_Level` — the honest "this is a
765    /// point, but not converged to the requested fidelity" status. This
766    /// catches the ill-conditioned / nlp_scaling-deflated case where the
767    /// scaled convergence test passes but the user-space duals have
768    /// drifted. It is a pure relabel at termination (no extra iterations);
769    /// unset or non-positive (the default) is a strict no-op, so every
770    /// existing caller keeps the Ipopt-faithful status.
771    fn apply_kkt_fidelity_gate(
772        &self,
773        app_status: ApplicationReturnStatus,
774    ) -> ApplicationReturnStatus {
775        if !matches!(app_status, ApplicationReturnStatus::SolveSucceeded) {
776            return app_status;
777        }
778        if let Ok((ftol, true)) = self.options.get_numeric_value("kkt_fidelity_tol", "") {
779            if ftol > 0.0 {
780                let unscaled_kkt = self.statistics.borrow().final_unscaled_kkt_error;
781                if unscaled_kkt > ftol {
782                    tracing::info!(target: "pounce::diagnostics",
783                        "kkt_fidelity_tol={ftol:.3e}: unscaled KKT error {unscaled_kkt:.3e} \
784                         exceeds it — downgrading Solve_Succeeded → \
785                         Solved_To_Acceptable_Level (pounce#173)");
786                    return ApplicationReturnStatus::SolvedToAcceptableLevel;
787                }
788            }
789        }
790        app_status
791    }
792
793    /// Build a *copy* of the algorithm builder configured per the
794    /// current options. The SQP path uses this so it gets a
795    /// fresh builder without mutating the application's state.
796    fn algorithm_builder_snapshot(&self) -> AlgorithmBuilder {
797        let mut builder = AlgorithmBuilder::default();
798        apply_sqp_options(&self.options, &mut builder.sqp);
799        apply_qp_subproblem_options(&self.options, &mut builder.sqp_qp);
800        builder
801    }
802
803    /// Construct a LinearBackendFactory honoring the
804    /// `linear_solver` option. Default FERAL; HSL MA57 when
805    /// built with the `ma57` feature.
806    fn make_backend_factory(&self) -> LinearBackendFactory {
807        Box::new(
808            |_choice| -> Box<dyn pounce_linsol::SparseSymLinearSolverInterface> {
809                Box::new(pounce_feral::FeralSolverInterface::new())
810            },
811        )
812    }
813
814    /// Phase 3.5 auto-fallback driver.
815    ///
816    /// Runs the standard solve (no wrapper) first. If it ends in a
817    /// trigger-class status (`Restoration_Failed`, `Infeasible_Problem_Detected`,
818    /// `Solved_To_Acceptable_Level`, `Maximum_Iterations_Exceeded`, or
819    /// `Not_Enough_Degrees_Of_Freedom`), retries transparently with
820    /// the ℓ₁ wrapper enabled. Promotes the retry's status only if
821    /// it returns `Solve_Succeeded`; otherwise returns the original
822    /// status.
823    ///
824    /// Caveat: the user TNLP's `finalize_solution` runs once per
825    /// attempt. When the retry doesn't promote, the user's captured
826    /// fields hold the retry's iterate (the ℓ₁-best least-infeasible
827    /// point) even though the returned status is the original's.
828    /// Documented on the option's help text; tightening this is a
829    /// Phase-4 follow-up.
830    fn run_with_l1_fallback(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
831        // First attempt: the standard IPM solve, no ℓ₁ wrapper. Only
832        // reached for `m > 0`, so `optimize_constrained` is exact.
833        let first_status = self.optimize_constrained(Rc::clone(&tnlp));
834        if !is_l1_fallback_trigger(first_status) {
835            return first_status;
836        }
837        // Trigger fired. Flip the wrapper option for the retry and
838        // restore it after — keeps the user's option-table view of the
839        // session exactly as they left it.
840        let prev = self
841            .options
842            .get_string_value("l1_exact_penalty_barrier", "")
843            .ok();
844        let _ = self
845            .options
846            .set_string_value("l1_exact_penalty_barrier", "yes", true, false);
847        let retry_status = self
848            .run_l1_penalty_outer_loop(Rc::clone(&tnlp))
849            .unwrap_or(ApplicationReturnStatus::InternalError);
850        let _ = self.options.set_string_value(
851            "l1_exact_penalty_barrier",
852            prev.as_ref().map(|(v, _)| v.as_str()).unwrap_or("no"),
853            true,
854            false,
855        );
856        if matches!(retry_status, ApplicationReturnStatus::SolveSucceeded) {
857            retry_status
858        } else {
859            first_status
860        }
861    }
862
863    /// μ-strategy auto-fallback driver (pounce#138).
864    ///
865    /// Runs the standard solve first. If it stalls short of optimal in a
866    /// way a μ-strategy flip can plausibly fix — `Solved_To_Acceptable_Level`
867    /// or `Maximum_Iterations_Exceeded`, the two signatures seen on the
868    /// princetonlib instances where the dual infeasibility parks above
869    /// `tol` while constraint violation and complementarity are already
870    /// deeply converged — it flips `mu_strategy` (adaptive↔monotone) and
871    /// solves once more. The retry's status is promoted only if it returns
872    /// `Solve_Succeeded`; otherwise the original status is returned.
873    ///
874    /// (maxcut/price stall at acceptable-level under adaptive; fermat2_vareps
875    /// stalls at `max_iter` — hence both triggers. flosp2tm is μ-independent
876    /// and correctly does not promote.)
877    ///
878    /// The flip direction is taken from the option's current value:
879    /// `adaptive` → `monotone`, anything else (including absent, which
880    /// the builder treats as monotone) → `adaptive`. The option table is
881    /// restored to the user's original view afterward.
882    ///
883    /// Caveat (shared with the ℓ₁ fallback): the user TNLP's
884    /// `finalize_solution` runs once per attempt, so when the retry
885    /// doesn't promote the captured fields hold the retry's iterate.
886    fn run_with_mu_strategy_fallback(
887        &mut self,
888        tnlp: Rc<RefCell<dyn TNLP>>,
889    ) -> ApplicationReturnStatus {
890        let first_status = self.optimize_constrained(Rc::clone(&tnlp));
891        if !matches!(
892            first_status,
893            ApplicationReturnStatus::SolvedToAcceptableLevel
894                | ApplicationReturnStatus::MaximumIterationsExceeded
895        ) {
896            return first_status;
897        }
898        // Flip the strategy for one retry. The parser maps "adaptive" →
899        // Adaptive and every other value (incl. unset) → Monotone, so the
900        // opposite of an explicit "adaptive" is "monotone" and the
901        // opposite of anything else is "adaptive".
902        let prev = self.options.get_string_value("mu_strategy", "").ok();
903        let was_adaptive = prev
904            .as_ref()
905            .map(|(v, found)| *found && v == "adaptive")
906            .unwrap_or(false);
907        let flipped = if was_adaptive { "monotone" } else { "adaptive" };
908        let _ = self
909            .options
910            .set_string_value("mu_strategy", flipped, true, false);
911        let retry_status = self.optimize_constrained(Rc::clone(&tnlp));
912        // Restore the user's original option-table view.
913        let _ = self.options.set_string_value(
914            "mu_strategy",
915            prev.as_ref()
916                .filter(|(_, found)| *found)
917                .map(|(v, _)| v.as_str())
918                .unwrap_or("monotone"),
919            true,
920            false,
921        );
922        if matches!(retry_status, ApplicationReturnStatus::SolveSucceeded) {
923            retry_status
924        } else {
925            first_status
926        }
927    }
928
929    /// Phase-3 ℓ₁-exact penalty-barrier outer loop.
930    ///
931    /// Builds an [`L1PenaltyBarrierTnlp`] wrapper around the user
932    /// TNLP, runs the constrained IPM at the current ρ, escalates ρ
933    /// per Byrd-Nocedal-Waltz steering, and terminates on any of:
934    ///   - slack sum collapses (`Σ(p+n) ≤ l1_slack_tol`)
935    ///   - inner solve returns non-Optimal (escalation won't fix
936    ///     numerical / restoration failure at this ρ)
937    ///   - ρ already at `l1_penalty_max`
938    ///   - `l1_penalty_max_outer_iter` reached
939    ///
940    /// After the loop, if the inner status is `SolveSucceeded` or
941    /// `SolvedToAcceptableLevel` but slacks didn't collapse, override
942    /// to `Infeasible_Problem_Detected` — the returned point is the
943    /// ℓ₁-best least-infeasible iterate, which is informative even
944    /// though the original constraints are not satisfied.
945    ///
946    /// Returns `Some(status)` if the wrapper ran the solve, `None` if
947    /// wrapper construction failed (caller should fall through to the
948    /// standard dispatch path).
949    fn run_l1_penalty_outer_loop(
950        &mut self,
951        tnlp: Rc<RefCell<dyn TNLP>>,
952    ) -> Option<ApplicationReturnStatus> {
953        let rho_init = self.l1_penalty_init();
954        let rho_max = self.l1_penalty_max().max(rho_init);
955        let factor = self.l1_penalty_increase_factor().max(1.0);
956        let tau = self.l1_steering_factor();
957        let slack_tol = self.l1_slack_tol();
958        let max_outer = self.l1_penalty_max_outer_iter().max(1);
959
960        let mut wrapper = pounce_l1penalty::L1PenaltyBarrierTnlp::new(Rc::clone(&tnlp), rho_init)?;
961        if wrapper.m_eq() == 0 {
962            // Nothing to slack — let the standard dispatch path handle
963            // this TNLP unmodified.
964            return None;
965        }
966        wrapper.set_defer_inner_finalize(true);
967        let wrapper_rc = Rc::new(RefCell::new(wrapper));
968
969        let mut rho = rho_init;
970        let mut last_status = ApplicationReturnStatus::InternalError;
971        for _outer in 0..max_outer {
972            wrapper_rc.borrow_mut().set_rho(rho);
973            let dyn_tnlp: Rc<RefCell<dyn TNLP>> = wrapper_rc.clone();
974            last_status = self.optimize_constrained(dyn_tnlp);
975
976            let w = wrapper_rc.borrow();
977            if !w.has_solution() {
978                // Inner solve aborted before producing an iterate.
979                drop(w);
980                break;
981            }
982            let slack_sum = w.last_slack_sum();
983            let y_eq_inf = w.last_y_eq_inf_norm();
984            drop(w);
985
986            // Termination decisions.
987            let inner_ok = matches!(
988                last_status,
989                ApplicationReturnStatus::SolveSucceeded
990                    | ApplicationReturnStatus::SolvedToAcceptableLevel
991            );
992            if !inner_ok {
993                break;
994            }
995            if slack_sum.is_finite() && slack_sum <= slack_tol {
996                break;
997            }
998            if rho >= rho_max {
999                break;
1000            }
1001            // BNW steering: ρ_new = max(ρ·factor, τ·‖y_eq‖∞ + ε)
1002            let geom = rho * factor;
1003            let steer = tau * y_eq_inf + 1.0e-12;
1004            rho = geom.max(steer).min(rho_max);
1005        }
1006
1007        // Forward to the user's inner.finalize_solution exactly once.
1008        let w = wrapper_rc.borrow();
1009        if w.has_solution() {
1010            let x_trunc: Vec<Number> = w.last_x_trunc().to_vec();
1011            let lambda: Vec<Number> = w.last_lambda().to_vec();
1012            let z_l: Vec<Number> = w.last_z_l_trunc().to_vec();
1013            let z_u: Vec<Number> = w.last_z_u_trunc().to_vec();
1014            let solver_status = w.last_status().unwrap_or(SolverReturn::InternalError);
1015            let slack_sum = w.last_slack_sum();
1016            drop(w);
1017
1018            // Honest-infeasibility upgrade (Phase 3): if the inner
1019            // solve says SolveSucceeded / SolvedToAcceptableLevel but
1020            // the slacks did not collapse, the original problem is
1021            // locally infeasible at the returned point. Override the
1022            // application status; the user-visible Solution.status is
1023            // updated below to the matching SolverReturn so the inner
1024            // TNLP sees a consistent picture.
1025            let infeasible_certificate = matches!(
1026                last_status,
1027                ApplicationReturnStatus::SolveSucceeded
1028                    | ApplicationReturnStatus::SolvedToAcceptableLevel
1029            ) && slack_sum.is_finite()
1030                && slack_sum > slack_tol;
1031            let final_app_status = if infeasible_certificate {
1032                ApplicationReturnStatus::InfeasibleProblemDetected
1033            } else {
1034                last_status
1035            };
1036            let final_solver_status = if infeasible_certificate {
1037                SolverReturn::LocalInfeasibility
1038            } else {
1039                solver_status
1040            };
1041
1042            // Recompute f(x*) and c(x*) on the inner.
1043            let f_inner = tnlp
1044                .borrow_mut()
1045                .eval_f(&x_trunc, true)
1046                .unwrap_or(Number::NAN);
1047            let m = tnlp
1048                .borrow_mut()
1049                .get_nlp_info()
1050                .map(|i| i.m as usize)
1051                .unwrap_or(0);
1052            let mut g_inner = vec![0.0; m];
1053            if m > 0 {
1054                let _ = tnlp.borrow_mut().eval_g(&x_trunc, false, &mut g_inner);
1055            }
1056            tnlp.borrow_mut().finalize_solution(
1057                Solution {
1058                    status: final_solver_status,
1059                    x: &x_trunc,
1060                    z_l: &z_l,
1061                    z_u: &z_u,
1062                    g: &g_inner,
1063                    lambda: &lambda,
1064                    obj_value: f_inner,
1065                },
1066                &TnlpIpoptData::default(),
1067                &TnlpIpoptCq::default(),
1068            );
1069            return Some(final_app_status);
1070        }
1071        // No solution captured at all — pass the inner status through.
1072        Some(last_status)
1073    }
1074
1075    /// Constrained-NLP path: build adapter → OrigIpoptNlp → algorithm
1076    /// bundle, run `optimize`, populate statistics, and call
1077    /// `finalize_solution` on the user's TNLP.
1078    fn optimize_constrained(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
1079        let t_start = Instant::now();
1080
1081        // `print_user_options yes` — dump the OptionsList before the
1082        // solve. Mirrors `IpoptApplication::call_optimize` (upstream
1083        // calls `Jnlst().Printf(.., "%s", options_->PrintUserOptions())`).
1084        let print_opts = self
1085            .options
1086            .get_bool_value("print_user_options", "")
1087            .ok()
1088            .and_then(|(v, f)| f.then_some(v))
1089            .unwrap_or(false);
1090        if print_opts {
1091            print!(
1092                "\nList of user-set options:\n\n{}",
1093                self.options.print_user_options()
1094            );
1095        }
1096
1097        // `print_options_documentation yes` — dump the full registry
1098        // (every option with type, default, valid range/strings, and
1099        // long description) before the solve. Honors
1100        // `print_options_mode` (`text` / `latex` / `doxygen`; only
1101        // `text` is implemented today, the others fall through with a
1102        // one-line note) and `print_advanced_options`. Mirrors
1103        // upstream `IpoptApplication::call_optimize`'s
1104        // `print_options_documentation` branch and `Common/IpRegOptions.cpp`
1105        // `OutputOptionDocumentation`.
1106        let print_doc = self
1107            .options
1108            .get_bool_value("print_options_documentation", "")
1109            .ok()
1110            .and_then(|(v, f)| f.then_some(v))
1111            .unwrap_or(false);
1112        if print_doc {
1113            let mode = self
1114                .options
1115                .get_string_value("print_options_mode", "")
1116                .ok()
1117                .map(|(v, _)| PrintOptionsMode::from_tag(&v))
1118                .unwrap_or(PrintOptionsMode::Text);
1119            let advanced = self
1120                .options
1121                .get_bool_value("print_advanced_options", "")
1122                .ok()
1123                .map(|(v, _)| v)
1124                .unwrap_or(false);
1125            print!(
1126                "\n# Pounce options registry\n\n{}",
1127                self.reg_options.print_options_documentation(mode, advanced)
1128            );
1129        }
1130
1131        // Mint a fresh `TimingStatistics` for this solve — shared (via
1132        // `Rc`) with the data and the NLP below so every `eval_*` and
1133        // every iterate-phase records into the same accumulator. The
1134        // application keeps its own `Rc` so callers can read totals out
1135        // via [`Self::timing_stats`].
1136        let timing = Rc::new(TimingStatistics::new());
1137        *self.timing.borrow_mut() = Rc::clone(&timing);
1138        timing.overall_alg.start();
1139
1140        // Reset the linear-solver summary sink so back-to-back solves
1141        // don't bleed factor counters / extremal pivots into each
1142        // other. Surviving the lock failure with a debug-assert keeps
1143        // a poisoned mutex from sinking a release build that doesn't
1144        // even consume the summary.
1145        if let Ok(mut guard) = self.linsol_summary_sink.lock() {
1146            *guard = LinearSolverSummary::default();
1147        } else {
1148            debug_assert!(false, "linsol summary sink mutex poisoned");
1149        }
1150
1151        // Build adapter + Nlp. Honor `fixed_variable_treatment` (default
1152        // `make_parameter`; pounce additionally implements `relax_bounds`,
1153        // which the adapter also auto-selects as a fallback when
1154        // `make_parameter` would leave `n_x_var < n_c` — mirrors upstream
1155        // `IpTNLPAdapter.cpp:623-633`).
1156        let lo_inf = self
1157            .options
1158            .get_numeric_value("nlp_lower_bound_inf", "")
1159            .ok()
1160            .and_then(|(v, f)| f.then_some(v))
1161            .unwrap_or(DEFAULT_NLP_LOWER_BOUND_INF);
1162        let up_inf = self
1163            .options
1164            .get_numeric_value("nlp_upper_bound_inf", "")
1165            .ok()
1166            .and_then(|(v, f)| f.then_some(v))
1167            .unwrap_or(DEFAULT_NLP_UPPER_BOUND_INF);
1168        let fixed_treatment = match self
1169            .options
1170            .get_string_value("fixed_variable_treatment", "")
1171            .ok()
1172            .and_then(|(v, f)| f.then_some(v))
1173            .as_deref()
1174        {
1175            Some("relax_bounds") => FixedVarTreatment::RelaxBounds,
1176            // `make_constraint` / `make_parameter_nodual` not yet
1177            // implemented; fall back to `make_parameter` (auto-retry to
1178            // `relax_bounds` will still kick in if DOF runs short).
1179            _ => FixedVarTreatment::MakeParameter,
1180        };
1181        let adapter = match TNLPAdapter::new_with_options(
1182            Rc::clone(&tnlp),
1183            lo_inf,
1184            up_inf,
1185            fixed_treatment,
1186        ) {
1187            Ok(a) => Rc::new(RefCell::new(a)),
1188            Err(_) => {
1189                timing.overall_alg.end();
1190                return ApplicationReturnStatus::InvalidProblemDefinition;
1191            }
1192        };
1193        // Carry the user's constant `obj_scaling_factor` (default 1.0;
1194        // negative ⇒ maximize) into the NLP. Until pounce#128's
1195        // follow-up this option was registered but never read, so it
1196        // was silently a no-op — maximization diverged because the
1197        // algorithm minimized the unscaled objective.
1198        let obj_scaling_factor = self
1199            .options
1200            .get_numeric_value("obj_scaling_factor", "")
1201            .ok()
1202            .and_then(|(v, f)| f.then_some(v))
1203            .unwrap_or(1.0);
1204        let mut orig_nlp = match OrigIpoptNlp::new(
1205            Rc::clone(&adapter),
1206            Rc::new(ConstObjScaling(obj_scaling_factor)),
1207        ) {
1208            Ok(n) => n,
1209            Err(_) => {
1210                timing.overall_alg.end();
1211                return ApplicationReturnStatus::InternalError;
1212            }
1213        };
1214        orig_nlp.set_timing_stats(Rc::clone(&timing));
1215
1216        // Mirror upstream `OrigIpoptNLP::InitializeStructures` (IpOrigIpoptNLP.cpp:299):
1217        // bail out with NotEnoughDegreesOfFreedom when there are fewer free
1218        // variables than equality constraints. Without this gate, square /
1219        // over-determined systems push the algorithm into restoration on
1220        // iter 0 and exit Restoration_Failed instead of the cleaner DOF code.
1221        let n_x_var = orig_nlp.x_space().dim();
1222        let n_c = orig_nlp.c_space().dim();
1223        if n_x_var > 0 && n_x_var < n_c {
1224            timing.overall_alg.end();
1225            return ApplicationReturnStatus::NotEnoughDegreesOfFreedom;
1226        }
1227
1228        // Relax `x_L / x_U / d_L / d_U` by `bound_relax_factor` (default
1229        // 1e-8), capped by `constr_viol_tol` (default 1e-4). Matches
1230        // `OrigIpoptNLP::InitializeStructures` lines 343-358.
1231        let bound_relax_factor = self
1232            .options
1233            .get_numeric_value("bound_relax_factor", "")
1234            .ok()
1235            .and_then(|(v, f)| f.then_some(v))
1236            .unwrap_or(1e-8);
1237        let constr_viol_tol = self
1238            .options
1239            .get_numeric_value("constr_viol_tol", "")
1240            .ok()
1241            .and_then(|(v, f)| f.then_some(v))
1242            .unwrap_or(1e-4);
1243        orig_nlp.relax_bounds(bound_relax_factor, constr_viol_tol);
1244
1245        // Apply automatic NLP scaling per `nlp_scaling_method` option
1246        // (port of `OrigIpoptNLP::InitializeStructures` →
1247        // `NLPScalingObject::DetermineScaling`). Default is
1248        // `gradient-based` to match upstream Ipopt 3.14.
1249        let scaling_method = self
1250            .options
1251            .get_string_value("nlp_scaling_method", "")
1252            .ok()
1253            .and_then(|(v, f)| f.then_some(v))
1254            .unwrap_or_else(|| "gradient-based".to_string());
1255        let scaling_method = match scaling_method.as_str() {
1256            "none" => ScalingMethod::None,
1257            "gradient-based" => ScalingMethod::GradientBased,
1258            "user-scaling" => ScalingMethod::UserScaling,
1259            // `equilibration-based` is registered upstream but not yet
1260            // implemented in pounce; fall back to gradient-based (the
1261            // upstream default) to keep behavior predictable.
1262            _ => ScalingMethod::GradientBased,
1263        };
1264        let max_gradient = self
1265            .options
1266            .get_numeric_value("nlp_scaling_max_gradient", "")
1267            .ok()
1268            .and_then(|(v, f)| f.then_some(v))
1269            .unwrap_or(100.0);
1270        let min_value = self
1271            .options
1272            .get_numeric_value("nlp_scaling_min_value", "")
1273            .ok()
1274            .and_then(|(v, f)| f.then_some(v))
1275            .unwrap_or(1e-8);
1276        let obj_target_gradient = self
1277            .options
1278            .get_numeric_value("nlp_scaling_obj_target_gradient", "")
1279            .ok()
1280            .and_then(|(v, f)| f.then_some(v))
1281            .unwrap_or(0.0);
1282        let constr_target_gradient = self
1283            .options
1284            .get_numeric_value("nlp_scaling_constr_target_gradient", "")
1285            .ok()
1286            .and_then(|(v, f)| f.then_some(v))
1287            .unwrap_or(0.0);
1288        orig_nlp.determine_scaling_from_starting_point(
1289            scaling_method,
1290            max_gradient,
1291            min_value,
1292            obj_target_gradient,
1293            constr_target_gradient,
1294        );
1295
1296        let nlp_handle: Rc<RefCell<dyn IpoptNlp>> = Rc::new(RefCell::new(orig_nlp));
1297
1298        // Build the algorithm strategy bundle. Read coarse knobs from
1299        // the OptionsList where we have them; fall through to defaults
1300        // otherwise. The full upstream parsing surface (mu_strategy,
1301        // hessian_approximation, line_search_method, ...) is wired by
1302        // `AlgBuilder::RegisterOptions` in upstream — that registry
1303        // hookup lands as a follow-up; default builder is correct for
1304        // HS71-class problems.
1305        let builder = self.algorithm_builder_from_options();
1306
1307        // Linear-solver backend. The default factory is option-aware
1308        // — it reads the `feral_*` extension options off the same
1309        // `OptionsList` that drove the IPM-level builder above so
1310        // per-problem `.opt` files can flip backend knobs without
1311        // rebuilding pounce.
1312        let feral_cfg = feral_config_from_options(&self.options);
1313        let factory = self.linear_backend_factory.take().unwrap_or_else(|| {
1314            default_backend_factory_with_sink(feral_cfg, Arc::clone(&self.linsol_summary_sink))
1315        });
1316        let bundle = builder.build_with_backend(factory);
1317
1318        // Wire the data / cq pair around the NLP. Install the shared
1319        // `TimingStatistics` so the algorithm's iterate phases
1320        // (output, convergence, hessian, μ, search-direction,
1321        // line-search, accept) all record into the same accumulator
1322        // the application exposes via `timing_stats()`.
1323        let data: crate::ipopt_data::IpoptDataHandle = Rc::new(RefCell::new(AlgIpoptData::new()));
1324        data.borrow_mut().timing = Rc::clone(&timing);
1325        let cq: crate::ipopt_cq::IpoptCqHandle = Rc::new(RefCell::new(
1326            IpoptCalculatedQuantities::new(Rc::clone(&data), Rc::clone(&nlp_handle)),
1327        ));
1328        // Correction size for very small slacks (default mach_eps^{3/4});
1329        // drives the safe-slack bound-adjustment mechanism.
1330        if let Ok((v, true)) = self.options.get_numeric_value("slack_move", "") {
1331            cq.borrow_mut().slack_move = v;
1332        }
1333
1334        // Seed `data.curr` with a zero-valued iterate of the correct
1335        // dimensions. The `IterateInitializer` consumes these as its
1336        // template (it overwrites `x`, `s`, multipliers in place); we
1337        // just need the dim metadata.
1338        {
1339            let nlp_borrow = nlp_handle.borrow();
1340            let n_x = nlp_borrow.n();
1341            let n_s = nlp_borrow.m_ineq();
1342            let n_yc = nlp_borrow.m_eq();
1343            let n_yd = nlp_borrow.m_ineq();
1344            let n_zl = nlp_borrow.x_l().dim();
1345            let n_zu = nlp_borrow.x_u().dim();
1346            let n_vl = nlp_borrow.d_l().dim();
1347            let n_vu = nlp_borrow.d_u().dim();
1348            drop(nlp_borrow);
1349            let iv = IteratesVector::new(
1350                Rc::new(DenseVectorSpace::new(n_x).make_new_dense()),
1351                Rc::new(DenseVectorSpace::new(n_s).make_new_dense()),
1352                Rc::new(DenseVectorSpace::new(n_yc).make_new_dense()),
1353                Rc::new(DenseVectorSpace::new(n_yd).make_new_dense()),
1354                Rc::new(DenseVectorSpace::new(n_zl).make_new_dense()),
1355                Rc::new(DenseVectorSpace::new(n_zu).make_new_dense()),
1356                Rc::new(DenseVectorSpace::new(n_vl).make_new_dense()),
1357                Rc::new(DenseVectorSpace::new(n_vu).make_new_dense()),
1358            );
1359            data.borrow_mut().set_curr(iv);
1360        }
1361
1362        // Full primal-dual warm restart (debugger `resolve`): if a
1363        // captured iterate is queued, install it onto `data.curr` over
1364        // the placeholder so the `WarmStartIterateInitializer`'s
1365        // re-optimize branch (x already initialized) keeps it and only
1366        // clamps multipliers / sets target_mu — no cold re-seed from the
1367        // NLP. Skipped (with a warning) if the dimensions don't line up,
1368        // e.g. an option changed the problem structure between solves.
1369        if let Some(snap) = self.warm_start_iterate.take() {
1370            let dims_match = {
1371                let borrow = data.borrow();
1372                borrow
1373                    .curr
1374                    .as_ref()
1375                    .map(|c| iterates_dims(c) == iterates_dims(snap.iterates()))
1376                    .unwrap_or(false)
1377            };
1378            if dims_match {
1379                data.borrow_mut().set_curr(snap.iterates().clone());
1380                data.borrow_mut().curr_mu = snap.mu();
1381            } else {
1382                tracing::warn!(
1383                    target: "pounce::warm_start",
1384                    "debugger warm-restart iterate dimensions differ from the fresh \
1385                     solve; ignoring the captured iterate and seeding normally"
1386                );
1387            }
1388        }
1389
1390        let max_iter = self
1391            .options
1392            .get_integer_value("max_iter", "")
1393            .ok()
1394            .and_then(|(v, f)| f.then_some(v))
1395            .unwrap_or(3000);
1396        let tol = self
1397            .options
1398            .get_numeric_value("tol", "")
1399            .ok()
1400            .and_then(|(v, f)| f.then_some(v))
1401            .unwrap_or(1e-8);
1402        data.borrow_mut().tol = tol;
1403
1404        let mut alg = IpoptAlgorithm::new(data, cq, bundle)
1405            .with_nlp(Rc::clone(&nlp_handle))
1406            .with_tnlp(Rc::clone(&tnlp));
1407        // Mint a fresh restoration factory per inner solve if a
1408        // provider is configured (pounce#10 Phase 3). Falls back to
1409        // the legacy one-shot `restoration_factory` slot when no
1410        // provider is set, preserving single-shot caller behavior.
1411        if let Some(provider) = self.restoration_factory_provider.as_mut() {
1412            self.restoration_factory = Some(provider());
1413        }
1414        if let Some(factory) = self.restoration_factory.as_mut() {
1415            alg = alg.with_restoration(factory());
1416        }
1417        if let Some(diag) = self.diagnostics.as_ref() {
1418            alg = alg.with_diagnostics(Rc::clone(diag));
1419        }
1420        // Move the interactive debugger hook (if any) into the main
1421        // algorithm. Taken — not cloned — so it drives exactly this
1422        // solve; a subsequent solve must reinstall it.
1423        if let Some(hook) = self.debug_hook.take() {
1424            alg = alg.with_debug_hook(hook);
1425        }
1426        alg.max_iter = max_iter;
1427        // Honor `print_level == 0`: suppress the per-iteration table
1428        // that the algorithm writes straight to stdout. (The Phase-7
1429        // journalist surface respects `print_level` already; this is
1430        // the legacy direct-print site that needs the same gate.)
1431        if let Ok((v, found)) = self.options.get_integer_value("print_level", "") {
1432            if found && v <= 0 {
1433                alg.print_iter_output = false;
1434                // The nested restoration IPM is built inside the
1435                // restoration driver, not by `IpoptAlgorithm::new`, so
1436                // it never sees this gate unless we forward it.
1437                if let Some(resto) = alg.restoration.as_mut() {
1438                    resto.set_print_iter_output(false);
1439                }
1440            }
1441        }
1442
1443        // Per-iteration history (pounce#71): when requested, capture the
1444        // `pounce::iteration` events emitted during the solve into an
1445        // `IterRecord` trajectory via the observability collector layer.
1446        // This replaces the old in-loop `iter_history` accumulation; it
1447        // requires the collector to be installed in the active
1448        // subscriber (the CLI / Python / C frontends install it via
1449        // `pounce_observability::init_subscriber`; tests call
1450        // `init_for_tests`). The collector scopes out restoration
1451        // sub-solve iterations via the `restoration` span, so the
1452        // trajectory matches the previous behavior (outer iters only).
1453        let iter_capture = self
1454            .record_iter_history
1455            .then(pounce_observability::IterCaptureGuard::start);
1456
1457        let solver_status = alg.optimize();
1458
1459        let captured_iters = iter_capture.map(|g| g.finish()).unwrap_or_default();
1460        // Close the overall-algorithm timer on the success path. The
1461        // early-return arms above end it themselves before bailing out;
1462        // this one matches upstream `IpoptApplication::call_optimize`
1463        // (which calls `EndCpuTime()` on overall_alg right after
1464        // `Optimize` returns, regardless of solver_status).
1465        timing.overall_alg.end();
1466
1467        // Drain counters / iter count off the algorithm.
1468        {
1469            let mut stats = self.statistics.borrow_mut();
1470            {
1471                let d = alg.data.borrow();
1472                stats.iteration_count = d.iter_count;
1473                // Converged barrier parameter μ — threaded forward into a
1474                // warm-started corrector's `mu_init` / `warm_start_target_mu`
1475                // for predictor–corrector path following (pounce#86).
1476                stats.final_mu = d.curr_mu;
1477            }
1478            stats.total_wallclock_time_secs = t_start.elapsed().as_secs_f64();
1479            // Restoration-phase audit counters (pounce#12). Zero on
1480            // problems where restoration never fires; populated by
1481            // `IpoptAlgorithm::invoke_restoration`.
1482            stats.restoration_calls = alg.resto_calls;
1483            stats.restoration_inner_iters = alg.resto_inner_iters;
1484            stats.restoration_outer_iters = alg.resto_outer_iters;
1485            stats.restoration_wall_secs = alg.resto_wall_secs;
1486            stats.iterations = captured_iters;
1487            // Capture the final *scaled* objective at the algorithm's
1488            // (compressed `x_var`-space) iterate via the NLP: the
1489            // algorithm-side `eval_f` returns `f * obj_scale_factor`.
1490            // `final_objective` is seeded with it only as a best-effort
1491            // fallback; the success path below overwrites it with the
1492            // true unscaled objective from `finalize_via_orig_nlp`
1493            // (which evaluates the user TNLP directly).
1494            let curr_x = alg.data.borrow().curr.as_ref().map(|c| c.x.clone());
1495            if let Some(x) = curr_x {
1496                if let Ok(f) = try_eval_curr_f(&nlp_handle, &x) {
1497                    stats.final_objective = f;
1498                    stats.final_scaled_objective = f;
1499                }
1500            }
1501            // Final residuals straight off the cq cache. These mirror
1502            // the values upstream prints in its end-of-run summary
1503            // ("Dual infeasibility / Constraint violation /
1504            // Complementarity / Overall NLP error").
1505            let cq = alg.cq.borrow();
1506            stats.final_dual_inf = cq.curr_dual_infeasibility_max();
1507            stats.final_constr_viol = cq.curr_primal_infeasibility_max();
1508            // Infinity-norm complementarity, max over all four bound
1509            // blocks (s_xl·z_l, s_xu·z_u, s_sl·v_l, s_su·v_u). The
1510            // empty-bound blocks return `0` from amax(), so the max is
1511            // safe even when only one side has bounds.
1512            let compl = cq
1513                .curr_compl_x_l()
1514                .amax()
1515                .max(cq.curr_compl_x_u().amax())
1516                .max(cq.curr_compl_s_l().amax())
1517                .max(cq.curr_compl_s_u().amax());
1518            stats.final_compl = compl;
1519            stats.final_kkt_error = cq.curr_nlp_error();
1520            // Unscaled (user-space) counterparts — divide the nlp_scaling
1521            // back out so a consumer can verify the certificate in its own
1522            // units (pounce#173). Identical to the scaled fields when no
1523            // scaling is active.
1524            stats.final_unscaled_dual_inf = cq.curr_unscaled_dual_infeasibility_max();
1525            stats.final_unscaled_constr_viol = cq.curr_unscaled_primal_infeasibility_max();
1526            stats.final_unscaled_compl = cq.curr_unscaled_complementarity_max();
1527            stats.final_unscaled_kkt_error = cq.curr_unscaled_nlp_error();
1528        }
1529
1530        // Map SolverReturn → ApplicationReturnStatus per
1531        // MAIN_LOOP.md's exception table, then apply the opt-in
1532        // status-fidelity gate (pounce#173).
1533        let app_status = self.apply_kkt_fidelity_gate(solver_return_to_app_status(solver_status));
1534
1535        // On convergence, fire the user-supplied callback (post-optimal
1536        // sensitivity hook, pounce#16) before flowing back through
1537        // `finalize_via_orig_nlp`. Borrowed handles into the converged
1538        // KKT state stay alive for the duration of the closure.
1539        if matches!(
1540            app_status,
1541            ApplicationReturnStatus::SolveSucceeded
1542                | ApplicationReturnStatus::SolvedToAcceptableLevel
1543        ) {
1544            if let Some(cb) = self.on_converged.as_mut() {
1545                if let Some(sd) = alg.search_dir.as_mut() {
1546                    let pd = sd.pd_solver_rc();
1547                    cb(&alg.data, &alg.cq, &nlp_handle, pd);
1548                }
1549            }
1550        }
1551
1552        // Finalize: forward the final iterate to the user's TNLP. The
1553        // returned objective is evaluated on the *user* TNLP at the
1554        // unscaled iterate, so it overrides the scaled best-effort
1555        // value stashed in `final_objective` above (the algorithm-side
1556        // `eval_f` returns `f * obj_scale_factor`).
1557        match finalize_via_orig_nlp(&nlp_handle, &alg, solver_status, app_status, &tnlp) {
1558            Ok(f_unscaled) => {
1559                self.statistics.borrow_mut().final_objective = f_unscaled;
1560            }
1561            Err(()) => {
1562                // Couldn't finalize; keep the scaled fallback and
1563                // surface the original status.
1564            }
1565        }
1566
1567        // End-of-solve timing report. Gated on `print_timing_statistics`
1568        // (default "no"); mirrors upstream's
1569        // `IpoptApplication::call_optimize` →
1570        // `IpTimingStatistics::PrintAllValues` call site. The report
1571        // goes to stdout (for parity with the banner / iter-row output
1572        // path) and is also fanned out to the journalist so an
1573        // `output_file` attached via `Initialize` picks it up.
1574        let print_timing = self
1575            .options
1576            .get_bool_value("print_timing_statistics", "")
1577            .ok()
1578            .and_then(|(v, f)| f.then_some(v))
1579            .unwrap_or(false);
1580        if print_timing {
1581            let report = timing.report();
1582            print!("{}", report);
1583            use pounce_common::journalist::{JournalCategory, JournalLevel};
1584            self.journalist.print(
1585                JournalLevel::J_SUMMARY,
1586                JournalCategory::J_TIMING_STATISTICS,
1587                &report,
1588            );
1589        }
1590
1591        app_status
1592    }
1593
1594    /// Build an [`AlgorithmBuilder`] populated from the app's
1595    /// [`OptionsList`]. Public so callers wiring the restoration
1596    /// factory can hand the *inner* IPM a builder that mirrors the
1597    /// outer's `mu_strategy`/`mu_oracle`/line-search choices —
1598    /// matching upstream `IpAlgBuilder::BuildRestoIpoptAlgorithm`,
1599    /// which reads the same `mu_strategy` option with prefix `"resto."
1600    /// + prefix` and falls back to the outer setting.
1601    pub fn algorithm_builder_from_options(&self) -> AlgorithmBuilder {
1602        let mut builder = AlgorithmBuilder::new();
1603
1604        // `mehrotra_algorithm` is parsed first so its cascading
1605        // defaults (mu_strategy=adaptive, mu_oracle=probing) can be
1606        // overridden by an explicit user setting of those keys
1607        // below. Mirrors `IpAlgBuilder.cpp:Mehrotra`.
1608        let mut mehrotra_on = false;
1609        if let Ok((v, found)) = self.options.get_string_value("mehrotra_algorithm", "") {
1610            if found && v == "yes" {
1611                mehrotra_on = true;
1612                builder.mehrotra_algorithm = true;
1613                builder.mu_strategy = MuStrategyChoice::Adaptive;
1614                builder.mu_oracle = crate::mu::adaptive::MuOracleKind::Probing;
1615                // `accept_every_trial_step` short-circuits the alpha
1616                // loop / filter — Mehrotra steps would otherwise be
1617                // rejected by the filter on LP-shaped problems because
1618                // the barrier objective is non-monotone along the
1619                // corrector. Mirrors upstream `IpAlgBuilder.cpp:Mehrotra`.
1620                builder.line_search.accept_every_trial_step = true;
1621                // Aggressive iterate-push defaults (`SetNumericValueIfUnset`
1622                // in upstream). The explicit user parses below will
1623                // overwrite these if the user set them explicitly.
1624                builder.init.bound_push = 10.0;
1625                builder.init.bound_frac = 0.2;
1626                builder.init.slack_bound_push = 10.0;
1627                builder.init.slack_bound_frac = 0.2;
1628                builder.init.bound_mult_init_val = 10.0;
1629                builder.init.constr_mult_init_max = 0.0;
1630                // `alpha_for_y=bound_mult` — Mehrotra wants the
1631                // equality multipliers to advance with the dual
1632                // alpha so they stay in step with z/v. Mirrors
1633                // upstream `IpIpoptAlg.cpp:InitializeImpl`.
1634                builder.line_search.alpha_for_y =
1635                    crate::line_search::backtracking::AlphaForY::BoundMult;
1636                // `adaptive_mu_globalization=never-monotone-mode` —
1637                // upstream `IpIpoptAlg.cpp:148-154` enforces this:
1638                // Mehrotra disables the globalization switch entirely
1639                // (no fallback to monotone mode when convergence
1640                // stalls). Required for the unsafeguarded Mehrotra
1641                // path to function.
1642                builder.mu.adaptive_mu_globalization =
1643                    crate::mu::adaptive::AdaptiveMuGlobalization::NeverMonotoneMode;
1644                // `least_square_init_primal=yes` — upstream
1645                // `IpIpoptAlg.cpp:182` enables this for the Mehrotra
1646                // cascade. Replaces the user's starting `x` with the
1647                // min-norm primal that satisfies the linearized
1648                // equality+inequality constraints. Critical on
1649                // LP-shaped problems where the user's starting point
1650                // can be wildly infeasible (e.g. nuffield2_trap).
1651                builder.init.least_square_init_primal = true;
1652            }
1653        }
1654
1655        if let Ok((v, found)) = self.options.get_string_value("mu_strategy", "") {
1656            if found {
1657                let parsed = match v.as_str() {
1658                    "adaptive" => MuStrategyChoice::Adaptive,
1659                    _ => MuStrategyChoice::Monotone,
1660                };
1661                if mehrotra_on && matches!(parsed, MuStrategyChoice::Monotone) {
1662                    // Upstream Ipopt refuses this combination: Mehrotra
1663                    // needs an affine step every iter, which only the
1664                    // adaptive path computes. Keep adaptive and warn.
1665                    tracing::warn!(target: "pounce::algorithm",
1666                        "pounce: mehrotra_algorithm=yes requires \
1667                         mu_strategy=adaptive; ignoring \
1668                         mu_strategy=monotone."
1669                    );
1670                } else {
1671                    builder.mu_strategy = parsed;
1672                }
1673            }
1674        }
1675        if let Ok((v, found)) = self.options.get_string_value("mu_oracle", "") {
1676            if found {
1677                builder.mu_oracle = match v.as_str() {
1678                    "loqo" => crate::mu::adaptive::MuOracleKind::Loqo,
1679                    "probing" => crate::mu::adaptive::MuOracleKind::Probing,
1680                    _ => crate::mu::adaptive::MuOracleKind::QualityFunction,
1681                };
1682            }
1683        }
1684        if let Ok((v, found)) = self
1685            .options
1686            .get_string_value("adaptive_mu_globalization", "")
1687        {
1688            if found {
1689                use crate::mu::adaptive::AdaptiveMuGlobalization;
1690                builder.mu.adaptive_mu_globalization = match v.as_str() {
1691                    "kkt-error" => AdaptiveMuGlobalization::KktError,
1692                    "never-monotone-mode" => AdaptiveMuGlobalization::NeverMonotoneMode,
1693                    _ => AdaptiveMuGlobalization::ObjConstrFilter,
1694                };
1695            }
1696        }
1697        if let Ok((v, found)) = self.options.get_string_value("hessian_approximation", "") {
1698            if found {
1699                builder.hessian_approximation = match v.as_str() {
1700                    "limited-memory" => HessianApproxChoice::LimitedMemory,
1701                    _ => HessianApproxChoice::Exact,
1702                };
1703            }
1704        }
1705        // Limited-memory quasi-Newton update formula. Registered upstream
1706        // (`limited_memory_update_type`, IpLimMemQuasiNewtonUpdater.cpp) but
1707        // until now read nowhere on the IPM path — the updater was hard-wired
1708        // to Powell-damped BFGS. SR1 is honored too (the updater and the
1709        // low-rank/inertia path already handle its indefinite models).
1710        if let Ok((v, found)) = self
1711            .options
1712            .get_string_value("limited_memory_update_type", "")
1713        {
1714            if found {
1715                builder.limited_memory_update_type = match v.as_str() {
1716                    "sr1" => UpdateType::Sr1,
1717                    _ => UpdateType::Bfgs,
1718                };
1719            }
1720        }
1721        // Limited-memory history length (`limited_memory_max_history`).
1722        if let Ok((v, found)) = self
1723            .options
1724            .get_integer_value("limited_memory_max_history", "")
1725        {
1726            if found && v >= 0 {
1727                builder.limited_memory_max_history = v as Index;
1728            }
1729        }
1730        if let Ok((v, found)) = self.options.get_string_value("line_search_method", "") {
1731            if found {
1732                builder.line_search_method = match v.as_str() {
1733                    "cg-penalty" => LineSearchChoice::CgPenalty,
1734                    "penalty" => LineSearchChoice::Penalty,
1735                    _ => LineSearchChoice::Filter,
1736                };
1737            }
1738        }
1739        // `accept_every_trial_step` — direct user override. Parsed
1740        // after the Mehrotra cascade so an explicit `no` still wins.
1741        if let Ok((v, found)) = self.options.get_string_value("accept_every_trial_step", "") {
1742            if found {
1743                builder.line_search.accept_every_trial_step = v == "yes";
1744            }
1745        }
1746        // `alpha_for_y` — direct user override. Parsed after the
1747        // Mehrotra cascade so an explicit value still wins.
1748        if let Ok((v, found)) = self.options.get_string_value("alpha_for_y", "") {
1749            if found {
1750                use crate::line_search::backtracking::AlphaForY;
1751                builder.line_search.alpha_for_y = match v.as_str() {
1752                    "primal" => AlphaForY::Primal,
1753                    "bound-mult" | "bound_mult" => AlphaForY::BoundMult,
1754                    "full" => AlphaForY::Full,
1755                    "min" => AlphaForY::Min,
1756                    "max" => AlphaForY::Max,
1757                    "primal-and-full" | "dual-and-full" => AlphaForY::Primal,
1758                    _ => AlphaForY::Primal,
1759                };
1760            }
1761        }
1762        // `nlp_scaling_method` is consumed NLP-side in
1763        // `OrigIpoptNlp::determine_scaling_from_starting_point` (see the
1764        // `determine_scaling_from_starting_point` call earlier in this
1765        // method); there is no algorithm-side scaling strategy to wire.
1766
1767        // Unlike the other options here, we always honor the registry
1768        // value (not just when the user set it explicitly): the option
1769        // registry default is "ma57" but `AlgorithmBuilder::default`
1770        // has `linear_solver: Feral`, so gating on `found` would
1771        // silently route default runs through Feral while the banner
1772        // (and ipopt-compatible behavior) advertises MA57.
1773        if let Ok((v, _found)) = self.options.get_string_value("linear_solver", "") {
1774            builder.linear_solver = match v.as_str() {
1775                "ma57" => LinearSolverChoice::Ma57,
1776                _ => LinearSolverChoice::Feral,
1777            };
1778        }
1779
1780        // `linear_system_scaling` — symmetric scaling of the augmented
1781        // KKT matrix before factorization. Port of
1782        // `IpTSymLinearSolver.cpp:RegisterOptions` plumbing. Default
1783        // "none"; "ruiz" invokes the Ruiz-2001 symmetric ∞-norm
1784        // equilibration in `RuizTSymScalingMethod`. "mc19" and
1785        // "slack-based" are accepted by the registry but not yet
1786        // implemented at this layer; they fall back to no scaling
1787        // with a one-line stderr notice.
1788        if let Ok((v, found)) = self.options.get_string_value("linear_system_scaling", "") {
1789            if found {
1790                builder.linear_system_scaling = match v.as_str() {
1791                    "ruiz" => crate::alg_builder::LinearSystemScalingChoice::Ruiz,
1792                    "mc19" => crate::alg_builder::LinearSystemScalingChoice::Mc19,
1793                    _ => crate::alg_builder::LinearSystemScalingChoice::None,
1794                };
1795            }
1796        }
1797        if let Ok((v, found)) = self.options.get_bool_value("linear_scaling_on_demand", "") {
1798            if found {
1799                builder.linear_scaling_on_demand = v;
1800            }
1801        }
1802
1803        // Convergence tolerances (port of `IpOptErrorConvCheck.cpp`'s
1804        // `RegisterOptions` consumers). Defaults already match upstream
1805        // — only override when the user set the key explicitly.
1806        let read_num = |key: &str| -> Option<f64> {
1807            self.options
1808                .get_numeric_value(key, "")
1809                .ok()
1810                .and_then(|(v, f)| f.then_some(v))
1811        };
1812        let read_int = |key: &str| -> Option<i32> {
1813            self.options
1814                .get_integer_value(key, "")
1815                .ok()
1816                .and_then(|(v, f)| f.then_some(v))
1817        };
1818        if let Some(v) = read_num("tol") {
1819            builder.conv_check.tol = v;
1820        }
1821        if let Some(v) = read_num("dual_inf_tol") {
1822            builder.conv_check.dual_inf_tol = v;
1823        }
1824        if let Some(v) = read_num("constr_viol_tol") {
1825            builder.conv_check.constr_viol_tol = v;
1826        }
1827        if let Some(v) = read_num("compl_inf_tol") {
1828            builder.conv_check.compl_inf_tol = v;
1829        }
1830        if let Some(v) = read_int("max_iter") {
1831            builder.conv_check.max_iter = v;
1832        }
1833        if let Some(v) = read_num("max_cpu_time") {
1834            builder.conv_check.max_cpu_time = v;
1835        }
1836        if let Some(v) = read_num("max_wall_time") {
1837            builder.conv_check.max_wall_time = v;
1838        }
1839        if let Some(v) = read_num("acceptable_tol") {
1840            builder.conv_check.acceptable_tol = v;
1841        }
1842        if let Some(v) = read_num("acceptable_dual_inf_tol") {
1843            builder.conv_check.acceptable_dual_inf_tol = v;
1844        }
1845        if let Some(v) = read_num("acceptable_constr_viol_tol") {
1846            builder.conv_check.acceptable_constr_viol_tol = v;
1847        }
1848        if let Some(v) = read_num("acceptable_compl_inf_tol") {
1849            builder.conv_check.acceptable_compl_inf_tol = v;
1850        }
1851        if let Some(v) = read_num("acceptable_obj_change_tol") {
1852            builder.conv_check.acceptable_obj_change_tol = v;
1853        }
1854        if let Some(v) = read_int("acceptable_iter") {
1855            builder.conv_check.acceptable_iter = v;
1856        }
1857        if let Some(v) = read_num("infeas_stationarity_tol") {
1858            builder.conv_check.infeas_stationarity_tol = v;
1859        }
1860        if let Some(v) = read_num("infeas_viol_kappa") {
1861            builder.conv_check.infeas_viol_kappa = v;
1862        }
1863        if let Some(v) = read_int("infeas_max_streak") {
1864            builder.conv_check.infeas_max_streak = v;
1865        }
1866
1867        // Barrier-parameter (μ) options — consumers in
1868        // `IpMonotoneMuUpdate.cpp` / `IpAdaptiveMuUpdate.cpp`. Both
1869        // updaters share the same option names; the builder forwards
1870        // each into whichever strategy is assembled.
1871        if let Some(v) = read_num("mu_init") {
1872            builder.mu.mu_init = v;
1873        }
1874        if let Some(v) = read_num("mu_max") {
1875            builder.mu.mu_max = v;
1876        }
1877        if let Some(v) = read_num("mu_max_fact") {
1878            builder.mu.mu_max_fact = v;
1879        }
1880        if let Some(v) = read_num("mu_min") {
1881            builder.mu.mu_min = v;
1882        }
1883        if let Some(v) = read_num("mu_target") {
1884            builder.mu.mu_target = v;
1885        }
1886        if let Some(v) = read_num("mu_linear_decrease_factor") {
1887            builder.mu.mu_linear_decrease_factor = v;
1888        }
1889        if let Some(v) = read_num("mu_superlinear_decrease_power") {
1890            builder.mu.mu_superlinear_decrease_power = v;
1891        }
1892        if let Ok((v, found)) = self
1893            .options
1894            .get_string_value("mu_allow_fast_monotone_decrease", "")
1895        {
1896            if found {
1897                builder.mu.mu_allow_fast_monotone_decrease = v == "yes";
1898            }
1899        }
1900        if let Some(v) = read_num("barrier_tol_factor") {
1901            builder.mu.barrier_tol_factor = v;
1902        }
1903        if let Some(v) = read_num("sigma_max") {
1904            builder.mu.sigma_max = v;
1905        }
1906        if let Some(v) = read_num("sigma_min") {
1907            builder.mu.sigma_min = v;
1908        }
1909
1910        // Quality-function oracle knobs — consumers in
1911        // `IpQualityFunctionMuOracle.cpp:RegisterOptions`. Forwarded
1912        // to the oracle on every free-mode call.
1913        if let Ok((v, found)) = self
1914            .options
1915            .get_string_value("quality_function_norm_type", "")
1916        {
1917            if found {
1918                use crate::mu::oracle::quality_function::NormType;
1919                builder.mu.quality_function_norm_type = match v.as_str() {
1920                    "1-norm" => NormType::OneNorm,
1921                    "2-norm" => NormType::TwoNorm,
1922                    "max-norm" => NormType::MaxNorm,
1923                    _ => NormType::TwoNormSquared,
1924                };
1925            }
1926        }
1927        if let Ok((v, found)) = self
1928            .options
1929            .get_string_value("quality_function_centrality", "")
1930        {
1931            if found {
1932                use crate::mu::oracle::quality_function::CentralityType;
1933                builder.mu.quality_function_centrality = match v.as_str() {
1934                    "log" => CentralityType::LogCenter,
1935                    "reciprocal" => CentralityType::ReciprocalCenter,
1936                    "cubed-reciprocal" => CentralityType::CubedReciprocalCenter,
1937                    _ => CentralityType::None,
1938                };
1939            }
1940        }
1941        if let Ok((v, found)) = self
1942            .options
1943            .get_string_value("quality_function_balancing_term", "")
1944        {
1945            if found {
1946                use crate::mu::oracle::quality_function::BalancingTermType;
1947                builder.mu.quality_function_balancing_term = match v.as_str() {
1948                    "cubic" => BalancingTermType::CubicTerm,
1949                    _ => BalancingTermType::None,
1950                };
1951            }
1952        }
1953        if let Some(v) = read_int("quality_function_max_section_steps") {
1954            builder.mu.quality_function_max_section_steps = v;
1955        }
1956        if let Some(v) = read_num("quality_function_section_sigma_tol") {
1957            builder.mu.quality_function_section_sigma_tol = v;
1958        }
1959        if let Some(v) = read_num("quality_function_section_qf_tol") {
1960            builder.mu.quality_function_section_qf_tol = v;
1961        }
1962
1963        // `probing_iterate_quality_factor` — pounce-specific guard
1964        // (pounce#58) on the probing μ-oracle's input iterate. When
1965        // `curr_avrg_compl / curr_mu` exceeds this factor, the
1966        // μ-update layer signals restoration via
1967        // `IpoptData::request_resto` instead of letting probing
1968        // return `σ · mu_curr` ≫ previous μ. Default 1e4; set to ≤ 0
1969        // to disable. No upstream Ipopt counterpart.
1970        if let Some(v) = read_num("probing_iterate_quality_factor") {
1971            builder.mu.probing_iterate_quality_factor = v;
1972        }
1973
1974        // Adaptive-μ extras — consumers in
1975        // `IpAdaptiveMuUpdate.cpp:RegisterOptions`. Only active when
1976        // `mu_strategy=adaptive`.
1977        if let Some(v) = read_num("adaptive_mu_safeguard_factor") {
1978            builder.mu.adaptive_mu_safeguard_factor = v;
1979        }
1980        if let Some(v) = read_num("adaptive_mu_monotone_init_factor") {
1981            builder.mu.adaptive_mu_monotone_init_factor = v;
1982        }
1983        if let Ok((v, found)) = self
1984            .options
1985            .get_bool_value("adaptive_mu_restore_previous_iterate", "")
1986        {
1987            if found {
1988                builder.mu.adaptive_mu_restore_previous_iterate = v;
1989            }
1990        }
1991        if let Some(v) = read_int("adaptive_mu_kkterror_red_iters") {
1992            if v >= 0 {
1993                builder.mu.adaptive_mu_kkterror_red_iters = v as usize;
1994            }
1995        }
1996        if let Some(v) = read_num("adaptive_mu_kkterror_red_fact") {
1997            builder.mu.adaptive_mu_kkterror_red_fact = v;
1998        }
1999        if let Ok((v, found)) = self
2000            .options
2001            .get_string_value("adaptive_mu_kkt_norm_type", "")
2002        {
2003            if found {
2004                use crate::mu::adaptive::AdaptiveMuKktNorm;
2005                builder.mu.adaptive_mu_kkt_norm_type = match v.as_str() {
2006                    "1-norm" => AdaptiveMuKktNorm::OneNorm,
2007                    "2-norm" => AdaptiveMuKktNorm::TwoNorm,
2008                    "max-norm" => AdaptiveMuKktNorm::MaxNorm,
2009                    _ => AdaptiveMuKktNorm::TwoNormSquared,
2010                };
2011            }
2012        }
2013
2014        // Watchdog options — consumers in
2015        // `IpBacktrackingLineSearch.cpp:RegisterOptions`. Baked into
2016        // the `BacktrackingLineSearch` at build time.
2017        if let Some(v) = read_int("watchdog_shortened_iter_trigger") {
2018            builder.line_search.watchdog_shortened_iter_trigger = v;
2019        }
2020        if let Some(v) = read_int("watchdog_trial_iter_max") {
2021            builder.line_search.watchdog_trial_iter_max = v;
2022        }
2023        if let Some(v) = read_num("soft_resto_pderror_reduction_factor") {
2024            builder.line_search.soft_resto_pderror_reduction_factor = v;
2025        }
2026        if let Some(v) = read_int("max_soft_resto_iters") {
2027            builder.line_search.max_soft_resto_iters = v;
2028        }
2029
2030        // Iteration-output options — consumed by `OrigIterationOutput`.
2031        if let Some(v) = read_int("print_frequency_iter") {
2032            builder.output.print_frequency_iter = v;
2033        }
2034        if let Some(v) = read_num("print_frequency_time") {
2035            builder.output.print_frequency_time = v;
2036        }
2037        if let Ok((v, found)) = self.options.get_bool_value("print_info_string", "") {
2038            if found {
2039                builder.output.print_info_string = v;
2040            }
2041        }
2042        if let Ok((v, found)) = self.options.get_string_value("inf_pr_output", "") {
2043            if found {
2044                builder.output.inf_pr_output_internal = v == "internal";
2045            }
2046        }
2047
2048        // Warm-start options — consumed by `WarmStartIterateInitializer`
2049        // (port of `IpWarmStartIterateInitializer.cpp:RegisterOptions`).
2050        // `warm_start_init_point` is the toggle that picks between the
2051        // default (cold) and warm-start initializers; the remaining
2052        // knobs are baked onto the chosen initializer at build time.
2053        if let Ok((v, found)) = self.options.get_bool_value("warm_start_init_point", "") {
2054            if found {
2055                builder.warm_start_init_point = v;
2056            }
2057        }
2058        if let Ok((v, found)) = self.options.get_bool_value("warm_start_same_structure", "") {
2059            if found {
2060                builder.warm.same_structure = v;
2061            }
2062        }
2063        if let Some(v) = read_num("warm_start_bound_push") {
2064            builder.warm.bound_push = v;
2065        }
2066        if let Some(v) = read_num("warm_start_bound_frac") {
2067            builder.warm.bound_frac = v;
2068        }
2069        if let Some(v) = read_num("warm_start_slack_bound_push") {
2070            builder.warm.slack_bound_push = v;
2071        }
2072        if let Some(v) = read_num("warm_start_slack_bound_frac") {
2073            builder.warm.slack_bound_frac = v;
2074        }
2075        if let Some(v) = read_num("warm_start_mult_bound_push") {
2076            builder.warm.mult_bound_push = v;
2077        }
2078        if let Some(v) = read_num("warm_start_mult_init_max") {
2079            builder.warm.mult_init_max = v;
2080        }
2081        if let Some(v) = read_num("warm_start_target_mu") {
2082            builder.warm.target_mu = v;
2083        }
2084        if let Ok((v, found)) = self
2085            .options
2086            .get_string_value("warm_start_entire_iterate", "")
2087        {
2088            if found {
2089                builder.warm.entire_iterate = v == "yes";
2090            }
2091        }
2092
2093        // `DefaultIterateInitializer` knobs — parsed after the Mehrotra
2094        // cascade so explicit user values win
2095        // (mirrors upstream's `SetNumericValueIfUnset` semantics).
2096        if let Some(v) = read_num("bound_push") {
2097            builder.init.bound_push = v;
2098        }
2099        if let Some(v) = read_num("bound_frac") {
2100            builder.init.bound_frac = v;
2101        }
2102        if let Some(v) = read_num("slack_bound_push") {
2103            builder.init.slack_bound_push = v;
2104        }
2105        if let Some(v) = read_num("slack_bound_frac") {
2106            builder.init.slack_bound_frac = v;
2107        }
2108        if let Some(v) = read_num("constr_mult_init_max") {
2109            builder.init.constr_mult_init_max = v;
2110        }
2111        if let Some(v) = read_num("bound_mult_init_val") {
2112            builder.init.bound_mult_init_val = v;
2113        }
2114        if let Ok((v, found)) = self.options.get_string_value("bound_mult_init_method", "") {
2115            if found {
2116                builder.init.bound_mult_init_method = v;
2117            }
2118        }
2119        if let Ok((v, found)) = self
2120            .options
2121            .get_string_value("least_square_init_primal", "")
2122        {
2123            if found {
2124                builder.init.least_square_init_primal = v == "yes";
2125            }
2126        }
2127        builder
2128    }
2129}
2130
2131/// Map the integer `print_level` / `file_print_level` option to the
2132/// matching [`JournalLevel`] variant. Mirrors upstream's
2133/// `static_cast<EJournalLevel>(int_value)` with clamping.
2134/// The eight block dimensions of an iterate, in canonical order
2135/// (x, s, y_c, y_d, z_l, z_u, v_l, v_u). Used to guard the debugger's
2136/// warm-restart install against a structural mismatch between solves.
2137fn iterates_dims(c: &IteratesVector) -> [i32; 8] {
2138    [
2139        c.x.dim(),
2140        c.s.dim(),
2141        c.y_c.dim(),
2142        c.y_d.dim(),
2143        c.z_l.dim(),
2144        c.z_u.dim(),
2145        c.v_l.dim(),
2146        c.v_u.dim(),
2147    ]
2148}
2149
2150fn journal_level_from_int(v: i32) -> JournalLevel {
2151    match v.clamp(0, 12) {
2152        0 => JournalLevel::J_NONE,
2153        1 => JournalLevel::J_ERROR,
2154        2 => JournalLevel::J_STRONGWARNING,
2155        3 => JournalLevel::J_SUMMARY,
2156        4 => JournalLevel::J_WARNING,
2157        5 => JournalLevel::J_ITERSUMMARY,
2158        6 => JournalLevel::J_DETAILED,
2159        7 => JournalLevel::J_MOREDETAILED,
2160        8 => JournalLevel::J_VECTOR,
2161        9 => JournalLevel::J_MOREVECTOR,
2162        10 => JournalLevel::J_MATRIX,
2163        11 => JournalLevel::J_MOREMATRIX,
2164        _ => JournalLevel::J_ALL,
2165    }
2166}
2167
2168/// Default symmetric linear-solver factory, parameterized by the
2169/// pounce-extension FERAL knobs read off the application's
2170/// `OptionsList`.
2171///
2172/// FERAL (pure-Rust) is the shipping default. The HSL MA57 backend is
2173/// available when the `ma57` cargo feature is enabled; without it,
2174/// requesting `linear_solver = ma57` falls back to FERAL with a
2175/// warning printed by the journalist (see [`AlgorithmBuilder`]).
2176pub fn default_backend_factory(feral_cfg: pounce_feral::FeralConfig) -> LinearBackendFactory {
2177    Box::new(
2178        move |choice: LinearSolverChoice| -> Box<dyn SparseSymLinearSolverInterface> {
2179            match choice {
2180                LinearSolverChoice::Feral => Box::new(
2181                    pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone()),
2182                ),
2183                LinearSolverChoice::Ma57 => {
2184                    #[cfg(feature = "ma57")]
2185                    {
2186                        Box::new(pounce_hsl::Ma57SolverInterface::new())
2187                    }
2188                    #[cfg(not(feature = "ma57"))]
2189                    {
2190                        // ma57 feature not compiled in — fall back to FERAL.
2191                        Box::new(pounce_feral::FeralSolverInterface::with_config(
2192                            feral_cfg.clone(),
2193                        ))
2194                    }
2195                }
2196            }
2197        },
2198    )
2199}
2200
2201/// Sink-aware variant of [`default_backend_factory`]. Identical
2202/// dispatch, but the FERAL backend is constructed with a
2203/// `LinearSolverSummary` sink so [`IpoptApplication`] can read out
2204/// aggregate post-mortem stats (factor counts, fill ratio, extremal
2205/// pivots, final inertia) after the solve returns. MA57 ignores the
2206/// sink — the HSL backend doesn't carry the same instrumentation yet.
2207pub fn default_backend_factory_with_sink(
2208    feral_cfg: pounce_feral::FeralConfig,
2209    sink: Arc<Mutex<LinearSolverSummary>>,
2210) -> LinearBackendFactory {
2211    Box::new(
2212        move |choice: LinearSolverChoice| -> Box<dyn SparseSymLinearSolverInterface> {
2213            match choice {
2214                LinearSolverChoice::Feral => Box::new(
2215                    pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone())
2216                        .with_summary_sink(Arc::clone(&sink)),
2217                ),
2218                LinearSolverChoice::Ma57 => {
2219                    #[cfg(feature = "ma57")]
2220                    {
2221                        Box::new(pounce_hsl::Ma57SolverInterface::new())
2222                    }
2223                    #[cfg(not(feature = "ma57"))]
2224                    {
2225                        Box::new(
2226                            pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone())
2227                                .with_summary_sink(Arc::clone(&sink)),
2228                        )
2229                    }
2230                }
2231            }
2232        },
2233    )
2234}
2235
2236/// Read the `feral_*` extension options off `options`, falling
2237/// back to the env-var defaults baked into [`pounce_feral::FeralConfig::from_env`]
2238/// for any knob the caller did not set explicitly. The returned
2239/// config is what every default-factory invocation (main IPM and
2240/// restoration sub-IPM) consumes.
2241pub fn feral_config_from_options(
2242    options: &pounce_common::options_list::OptionsList,
2243) -> pounce_feral::FeralConfig {
2244    let mut cfg = pounce_feral::FeralConfig::from_env();
2245    // Tri-state: the `(_, true)` arm only fires when the user set the
2246    // option explicitly. Leaving it unset keeps `cfg.cascade_break` at
2247    // `None`, which inherits FERAL's `NumericParams::default()` (CB on
2248    // as of FERAL Phase B / pounce#55). `Some(false)` explicitly
2249    // disarms (reproduces pre-Phase-B behaviour, surfaces FERAL's
2250    // `DelayBudgetExceeded` on non-root cascade victims).
2251    if let Ok((v, true)) = options.get_bool_value("feral_cascade_break", "") {
2252        cfg.cascade_break = Some(v);
2253    }
2254    if let Ok((v, true)) = options.get_bool_value("feral_fma", "") {
2255        cfg.fma = v;
2256    }
2257    if let Ok((v, true)) = options.get_bool_value("feral_refine", "") {
2258        cfg.refine = v;
2259    }
2260    if let Ok((v, true)) = options.get_numeric_value("feral_singular_pivot_floor", "") {
2261        cfg.singular_pivot_floor = v;
2262    }
2263    if let Ok((v, true)) = options.get_numeric_value("feral_pivtol", "") {
2264        cfg.pivtol = v;
2265    }
2266    // Only override on explicit set so `from_env` (which itself
2267    // defaults to OrderingMethod::Auto) keeps governing unset cases.
2268    // Unrecognized tags are silently ignored — the registered enum
2269    // restricts inputs at the OptionsList layer.
2270    if let Ok((v, true)) = options.get_string_value("feral_ordering", "") {
2271        if let Some(m) = pounce_feral::parse_ordering_method(&v) {
2272            cfg.ordering = m;
2273        }
2274    }
2275    // Same explicit-set discipline as `feral_ordering`: `from_env`
2276    // defaults to ScalingStrategy::Auto (FERAL's current default), so
2277    // leaving the option unset preserves existing behaviour exactly.
2278    if let Ok((v, true)) = options.get_string_value("feral_scaling", "") {
2279        if let Some(s) = pounce_feral::parse_scaling_strategy(&v) {
2280            cfg.scaling = s;
2281        }
2282    }
2283    cfg
2284}
2285
2286/// Map upstream `SolverReturn` codes to `ApplicationReturnStatus`.
2287/// Mirrors the table in
2288/// `ref/Ipopt/AGENT_REFERENCE/MAIN_LOOP.md` ("exception → SolverReturn
2289/// map") and the corresponding switch in
2290/// `IpIpoptApplication.cpp:call_optimize`.
2291fn solver_return_to_app_status(s: SolverReturn) -> ApplicationReturnStatus {
2292    match s {
2293        SolverReturn::Success => ApplicationReturnStatus::SolveSucceeded,
2294        SolverReturn::StopAtAcceptablePoint => ApplicationReturnStatus::SolvedToAcceptableLevel,
2295        SolverReturn::FeasiblePointFound => ApplicationReturnStatus::FeasiblePointFound,
2296        SolverReturn::MaxiterExceeded => ApplicationReturnStatus::MaximumIterationsExceeded,
2297        SolverReturn::CpuTimeExceeded => ApplicationReturnStatus::MaximumCpuTimeExceeded,
2298        SolverReturn::WallTimeExceeded => ApplicationReturnStatus::MaximumWallTimeExceeded,
2299        SolverReturn::StopAtTinyStep => ApplicationReturnStatus::SearchDirectionBecomesTooSmall,
2300        SolverReturn::LocalInfeasibility => ApplicationReturnStatus::InfeasibleProblemDetected,
2301        SolverReturn::UserRequestedStop => ApplicationReturnStatus::UserRequestedStop,
2302        SolverReturn::DivergingIterates => ApplicationReturnStatus::DivergingIterates,
2303        SolverReturn::RestorationFailure => ApplicationReturnStatus::RestorationFailed,
2304        SolverReturn::ErrorInStepComputation => ApplicationReturnStatus::ErrorInStepComputation,
2305        SolverReturn::InvalidNumberDetected => ApplicationReturnStatus::InvalidNumberDetected,
2306        SolverReturn::TooFewDegreesOfFreedom => ApplicationReturnStatus::NotEnoughDegreesOfFreedom,
2307        SolverReturn::InvalidOption => ApplicationReturnStatus::InvalidOption,
2308        SolverReturn::OutOfMemory => ApplicationReturnStatus::InsufficientMemory,
2309        SolverReturn::InternalError | SolverReturn::Unassigned => {
2310            ApplicationReturnStatus::InternalError
2311        }
2312    }
2313}
2314
2315/// Best-effort evaluation of the objective at the algorithm's final
2316/// `x`. Returns the *scaled* objective (`f * obj_scale_factor`); used
2317/// to populate `SolveStatistics::final_scaled_objective`.
2318fn try_eval_curr_f(
2319    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2320    x: &Rc<dyn pounce_linalg::Vector>,
2321) -> Result<Number, ()> {
2322    let mut nlp_mut = nlp.borrow_mut();
2323    Ok(nlp_mut.eval_f(&**x))
2324}
2325
2326/// Trigger predicate for the Phase-3.5 ℓ₁ auto-fallback path. Returns
2327/// `true` when a status warrants a retry through the wrapper. Mirrors
2328/// ripopt#23's trigger set, extended per the audit's Refinement B
2329/// (pounce-side `Not_Enough_Degrees_Of_Freedom` is added because
2330/// pounce's DOF early-exit blocks NE-suffix problems that ripopt's
2331/// equivalent would let pass to the wrapper).
2332fn is_l1_fallback_trigger(status: ApplicationReturnStatus) -> bool {
2333    matches!(
2334        status,
2335        ApplicationReturnStatus::RestorationFailed
2336            | ApplicationReturnStatus::InfeasibleProblemDetected
2337            | ApplicationReturnStatus::SolvedToAcceptableLevel
2338            | ApplicationReturnStatus::MaximumIterationsExceeded
2339            | ApplicationReturnStatus::NotEnoughDegreesOfFreedom
2340    )
2341}
2342
2343/// Forward the final iterate back to the user's `TNLP::finalize_solution`.
2344/// We pull `x` (compressed in `x_var`-space) off the algorithm's
2345/// `data.curr`, lift it back to full TNLP indexing, and pass empty
2346/// multipliers for now (the algorithm's `y_c`, `y_d`, `z_l`, `z_u` are
2347/// in compressed split form — re-assembling them into the user's
2348/// `lambda` / `z_l` / `z_u` is mechanical but lives behind a
2349/// `OrigIpoptNlp::finalize_solution_*` accessor that's still being
2350/// fleshed out). On success returns the unscaled objective evaluated
2351/// on the user TNLP at the final iterate; returns `Err` if the final
2352/// iterate is missing.
2353fn finalize_via_orig_nlp(
2354    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2355    alg: &IpoptAlgorithm,
2356    solver_status: SolverReturn,
2357    _app_status: ApplicationReturnStatus,
2358    tnlp: &Rc<RefCell<dyn TNLP>>,
2359) -> Result<Number, ()> {
2360    let curr = alg.data.borrow().curr.clone().ok_or(())?;
2361    // Lift compressed x_var → full-x (length `info.n`) so the user
2362    // TNLP receives the same shape it provided. With `make_parameter`
2363    // the fixed components are spliced back in by the IpoptNlp.
2364    let nlp_borrow = nlp.borrow();
2365    let x_vec: Vec<Number> = nlp_borrow.lift_x_to_full(&*curr.x);
2366    let info = tnlp.borrow_mut().get_nlp_info().ok_or(())?;
2367    let n = info.n as usize;
2368    let m = info.m as usize;
2369    debug_assert_eq!(x_vec.len(), n);
2370    // Lift algorithm-side multipliers back into user-space (pounce#11).
2371    // Use the `finalize_solution_*` family (not the `pack_*` family): the
2372    // final solution duals must be reported in the user's *unscaled-
2373    // Lagrangian* convention `∇f + λ·∇g + z = 0`, which divides out the
2374    // `obj_scale_factor` the algorithm threads through `eval_h`. The `pack_*`
2375    // family deliberately omits that division because it feeds the scaled
2376    // `eval_h`; calling it here left every dual scaled by `obj_scale_factor`
2377    // whenever gradient-based scaling triggered (pounce#11 F1).
2378    // Backends without overrides return empty; fall back to zero stubs so the
2379    // user sees a length-consistent vector.
2380    let mut z_l = nlp_borrow.finalize_solution_z_l(&*curr.z_l);
2381    if z_l.is_empty() {
2382        z_l = vec![0.0; n];
2383    }
2384    let mut z_u = nlp_borrow.finalize_solution_z_u(&*curr.z_u);
2385    if z_u.is_empty() {
2386        z_u = vec![0.0; n];
2387    }
2388    let mut lambda = nlp_borrow.finalize_solution_lambda(&*curr.y_c, &*curr.y_d);
2389    if lambda.is_empty() {
2390        lambda = vec![0.0; m];
2391    }
2392    drop(nlp_borrow);
2393    // Compute g(x) via the user TNLP so the final residual is
2394    // populated for the user.
2395    let mut g_final = vec![0.0; m];
2396    let _ = tnlp.borrow_mut().eval_g(&x_vec, true, &mut g_final);
2397    let f_final = tnlp
2398        .borrow_mut()
2399        .eval_f(&x_vec, true)
2400        .unwrap_or(Number::NAN);
2401    tnlp.borrow_mut().finalize_solution(
2402        Solution {
2403            status: solver_status,
2404            x: &x_vec,
2405            z_l: &z_l,
2406            z_u: &z_u,
2407            g: &g_final,
2408            lambda: &lambda,
2409            obj_value: f_final,
2410        },
2411        &TnlpIpoptData::default(),
2412        &TnlpIpoptCq::default(),
2413    );
2414    Ok(f_final)
2415}
2416
2417/// Bind SQP suboptions registered in `upstream_options.rs`
2418/// (`sqp_globalization`, `sqp_hessian`, `sqp_max_iter`, `sqp_tol`,
2419/// `sqp_constr_viol_tol`, `sqp_dual_inf_tol`, `sqp_l1_penalty`,
2420/// `sqp_bt_reduction`, `sqp_bt_min_alpha`, `sqp_print_level`,
2421/// `sqp_lbfgs_max_history`) onto
2422/// `opts`. Used by [`IpoptApplication::algorithm_builder_snapshot`]
2423/// before constructing an SQP algorithm.
2424fn apply_sqp_options(options: &OptionsList, opts: &mut crate::sqp::SqpOptions) {
2425    use crate::sqp::{SqpGlobalization, SqpHessianSource};
2426
2427    if let Ok((s, true)) = options.get_string_value("sqp_globalization", "") {
2428        opts.globalization = match s.as_str() {
2429            "filter" => SqpGlobalization::Filter,
2430            "l1-elastic" => SqpGlobalization::L1Elastic,
2431            _ => opts.globalization,
2432        };
2433    }
2434    if let Ok((s, true)) = options.get_string_value("sqp_hessian", "") {
2435        opts.hessian = match s.as_str() {
2436            "exact" => SqpHessianSource::Exact,
2437            "damped-bfgs" => SqpHessianSource::DampedBfgs,
2438            "lbfgs" => SqpHessianSource::Lbfgs,
2439            _ => opts.hessian,
2440        };
2441    }
2442    if let Ok((v, true)) = options.get_integer_value("sqp_max_iter", "") {
2443        if v >= 0 {
2444            opts.max_iter = v as u32;
2445        }
2446    }
2447    if let Ok((v, true)) = options.get_numeric_value("sqp_tol", "") {
2448        opts.tol = v;
2449    }
2450    if let Ok((v, true)) = options.get_numeric_value("sqp_constr_viol_tol", "") {
2451        opts.constr_viol_tol = v;
2452    }
2453    if let Ok((v, true)) = options.get_numeric_value("sqp_dual_inf_tol", "") {
2454        opts.dual_inf_tol = v;
2455    }
2456    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty", "") {
2457        opts.l1_penalty = v;
2458    }
2459    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty_safety", "") {
2460        opts.l1_penalty_safety = v;
2461    }
2462    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty_max", "") {
2463        opts.l1_penalty_max = v;
2464    }
2465    if let Ok((v, true)) = options.get_numeric_value("sqp_bt_reduction", "") {
2466        opts.bt_reduction = v;
2467    }
2468    if let Ok((v, true)) = options.get_numeric_value("sqp_bt_min_alpha", "") {
2469        opts.bt_min_alpha = v;
2470    }
2471    if let Ok((v, true)) = options.get_integer_value("sqp_print_level", "") {
2472        opts.print_level = v.clamp(0, u8::MAX as i32) as u8;
2473    }
2474    if let Ok((v, true)) = options.get_integer_value("sqp_lbfgs_max_history", "") {
2475        if v >= 1 {
2476            opts.lbfgs_max_history = v as u32;
2477        }
2478    }
2479}
2480
2481/// Populate the active-set SQP **QP-subproblem** options
2482/// ([`pounce_qp::QpOptions`]) from the `sqp_qp_*` option family.
2483///
2484/// Sister to [`apply_sqp_options`], which handles the SQP *outer-loop*
2485/// options ([`crate::sqp::SqpOptions`]); this one feeds the inner QP
2486/// solver that `SqpAlgorithm` delegates each subproblem to. Consulted
2487/// only on the `ActiveSetSqp` path. Each knob is forwarded only when
2488/// the user explicitly set it (the `true` flag), so the `pounce_qp`
2489/// defaults stand otherwise.
2490fn apply_qp_subproblem_options(options: &OptionsList, opts: &mut pounce_qp::QpOptions) {
2491    use pounce_qp::AntiCyclingChoice;
2492
2493    if let Ok((v, true)) = options.get_integer_value("sqp_qp_max_iter", "") {
2494        if v >= 0 {
2495            opts.max_iter = v as u32;
2496        }
2497    }
2498    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_feas_tol", "") {
2499        opts.feas_tol = v;
2500    }
2501    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_opt_tol", "") {
2502        opts.opt_tol = v;
2503    }
2504    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_elastic_gamma", "") {
2505        opts.elastic_gamma = v;
2506    }
2507    if let Ok((s, true)) = options.get_string_value("sqp_qp_anti_cycling", "") {
2508        opts.anti_cycling = match s.as_str() {
2509            "expand" => AntiCyclingChoice::Expand,
2510            "bland" => AntiCyclingChoice::Bland,
2511            "none" => AntiCyclingChoice::None,
2512            _ => opts.anti_cycling,
2513        };
2514    }
2515}
2516
2517/// SQP-side analog of [`finalize_via_orig_nlp`]. Hands the SQP
2518/// solution iterate to the user TNLP via the standard
2519/// `finalize_solution` callback. Multiplier lifting goes through
2520/// the same OrigIpoptNlp hooks so the user sees the same shape
2521/// regardless of which algorithm produced the iterate.
2522///
2523/// Returns the user-space objective value on success.
2524fn finalize_via_sqp(
2525    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2526    res: &crate::sqp::SqpResult,
2527    solver_status: pounce_nlp::SolverReturn,
2528    tnlp: &Rc<RefCell<dyn TNLP>>,
2529) -> Result<Number, ()> {
2530    use pounce_linalg::dense_vector::DenseVectorSpace;
2531
2532    let info = tnlp.borrow_mut().get_nlp_info().ok_or(())?;
2533    let n = info.n as usize;
2534    let m = info.m as usize;
2535
2536    // Wrap SQP slices in DenseVectors so we can pass them through
2537    // the OrigIpoptNlp lift_x_to_full / pack_*_for_user hooks.
2538    let nlp_borrow = nlp.borrow();
2539    let n_alg = nlp_borrow.n() as usize;
2540    let m_eq = nlp_borrow.m_eq() as usize;
2541    let m_ineq = nlp_borrow.m_ineq() as usize;
2542    debug_assert_eq!(res.x.len(), n_alg);
2543    debug_assert_eq!(res.lambda_g.len(), m_eq + m_ineq);
2544    debug_assert_eq!(res.lambda_x.len(), n_alg);
2545
2546    let x_space = DenseVectorSpace::new(n_alg as Index);
2547    let c_space = DenseVectorSpace::new(m_eq as Index);
2548    let d_space = DenseVectorSpace::new(m_ineq as Index);
2549
2550    let mut x_dv = x_space.make_new_dense();
2551    x_dv.set_values(&res.x);
2552    let x_vec: Vec<Number> = nlp_borrow.lift_x_to_full(&x_dv);
2553    debug_assert_eq!(x_vec.len(), n);
2554
2555    // λ_x is packed signed (z_l − z_u). Split for lift.
2556    let mut z_l_compressed = x_space.make_new_dense();
2557    let mut z_u_compressed = x_space.make_new_dense();
2558    let zl_vals: Vec<Number> = res.lambda_x.iter().map(|v| v.max(0.0)).collect();
2559    let zu_vals: Vec<Number> = res.lambda_x.iter().map(|v| (-v).max(0.0)).collect();
2560    z_l_compressed.set_values(&zl_vals);
2561    z_u_compressed.set_values(&zu_vals);
2562    // `finalize_solution_*` (not `pack_*`): report unscaled-Lagrangian duals,
2563    // dividing out `obj_scale_factor` — see `finalize_via_orig_nlp` (F1).
2564    let mut z_l = nlp_borrow.finalize_solution_z_l(&z_l_compressed);
2565    if z_l.is_empty() {
2566        z_l = vec![0.0; n];
2567    }
2568    let mut z_u = nlp_borrow.finalize_solution_z_u(&z_u_compressed);
2569    if z_u.is_empty() {
2570        z_u = vec![0.0; n];
2571    }
2572
2573    // λ_g is [y_c; y_d]; split into the c/d blocks for lift.
2574    let mut y_c_dv = c_space.make_new_dense();
2575    let mut y_d_dv = d_space.make_new_dense();
2576    if m_eq > 0 {
2577        y_c_dv.set_values(&res.lambda_g[..m_eq]);
2578    }
2579    if m_ineq > 0 {
2580        y_d_dv.set_values(&res.lambda_g[m_eq..]);
2581    }
2582    let mut lambda = nlp_borrow.finalize_solution_lambda(&y_c_dv, &y_d_dv);
2583    if lambda.is_empty() {
2584        lambda = vec![0.0; m];
2585    }
2586    drop(nlp_borrow);
2587
2588    let mut g_final = vec![0.0; m];
2589    let _ = tnlp.borrow_mut().eval_g(&x_vec, true, &mut g_final);
2590    let f_final = tnlp
2591        .borrow_mut()
2592        .eval_f(&x_vec, true)
2593        .unwrap_or(Number::NAN);
2594    tnlp.borrow_mut().finalize_solution(
2595        pounce_nlp::tnlp::Solution {
2596            status: solver_status,
2597            x: &x_vec,
2598            z_l: &z_l,
2599            z_u: &z_u,
2600            g: &g_final,
2601            lambda: &lambda,
2602            obj_value: f_final,
2603        },
2604        &TnlpIpoptData::default(),
2605        &TnlpIpoptCq::default(),
2606    );
2607    Ok(f_final)
2608}
2609
2610#[cfg(test)]
2611mod tests {
2612    use super::*;
2613    use pounce_nlp::tnlp::{
2614        BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest,
2615        StartingPoint,
2616    };
2617
2618    struct Hs071Stub;
2619    impl TNLP for Hs071Stub {
2620        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
2621            // HS071 dimensions: n=4, m=2, dense Jacobian (8 nz),
2622            // dense lower-triangular Hessian (10 nz).
2623            Some(NlpInfo {
2624                n: 4,
2625                m: 2,
2626                nnz_jac_g: 8,
2627                nnz_h_lag: 10,
2628                index_style: IndexStyle::C,
2629            })
2630        }
2631        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
2632            b.x_l.copy_from_slice(&[1.0; 4]);
2633            b.x_u.copy_from_slice(&[5.0; 4]);
2634            b.g_l.copy_from_slice(&[25.0, 40.0]);
2635            b.g_u.copy_from_slice(&[2.0e19, 40.0]);
2636            true
2637        }
2638        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
2639            sp.x.copy_from_slice(&[1.0, 5.0, 5.0, 1.0]);
2640            true
2641        }
2642        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
2643            Some(x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2])
2644        }
2645        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
2646            grad.fill(0.0);
2647            true
2648        }
2649        fn eval_g(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
2650            g.fill(0.0);
2651            true
2652        }
2653        fn eval_jac_g(
2654            &mut self,
2655            _x: Option<&[Number]>,
2656            _new_x: bool,
2657            mode: SparsityRequest<'_>,
2658        ) -> bool {
2659            if let SparsityRequest::Structure { irow, jcol } = mode {
2660                irow.copy_from_slice(&[0, 0, 0, 0, 1, 1, 1, 1]);
2661                jcol.copy_from_slice(&[0, 1, 2, 3, 0, 1, 2, 3]);
2662            }
2663            true
2664        }
2665        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
2666    }
2667
2668    #[test]
2669    fn application_default_does_not_select_sqp() {
2670        let mut app = IpoptApplication::new();
2671        app.initialize().unwrap();
2672        assert!(!app.is_sqp_algorithm_selected());
2673    }
2674
2675    #[test]
2676    fn application_routes_to_sqp_when_algorithm_option_set() {
2677        let mut app = IpoptApplication::new();
2678        app.initialize().unwrap();
2679        app.initialize_with_options_str("algorithm active-set-sqp\n")
2680            .unwrap();
2681        assert!(app.is_sqp_algorithm_selected());
2682    }
2683
2684    /// Convex equality NLP fixture for end-to-end SQP testing
2685    /// through `IpoptApplication`:
2686    ///
2687    ///     min ½(x₁² + x₂²) − x₁ − 2x₂  s.t.  x₁ + x₂ = 1
2688    ///
2689    /// Closed form: x* = (0, 1), obj = -1.5, λ_g = 1.
2690    struct ConvexEqTnlp {
2691        finalize_called: std::rc::Rc<std::cell::RefCell<Option<(Vec<Number>, Number)>>>,
2692    }
2693    impl TNLP for ConvexEqTnlp {
2694        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
2695            Some(NlpInfo {
2696                n: 2,
2697                m: 1,
2698                nnz_jac_g: 2,
2699                nnz_h_lag: 2,
2700                index_style: IndexStyle::C,
2701            })
2702        }
2703        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
2704            b.x_l.copy_from_slice(&[-2.0e19; 2]);
2705            b.x_u.copy_from_slice(&[2.0e19; 2]);
2706            b.g_l.copy_from_slice(&[1.0]);
2707            b.g_u.copy_from_slice(&[1.0]);
2708            true
2709        }
2710        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
2711            sp.x.copy_from_slice(&[0.0, 0.0]);
2712            true
2713        }
2714        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
2715            Some(0.5 * (x[0] * x[0] + x[1] * x[1]) - x[0] - 2.0 * x[1])
2716        }
2717        fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
2718            grad[0] = x[0] - 1.0;
2719            grad[1] = x[1] - 2.0;
2720            true
2721        }
2722        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
2723            g[0] = x[0] + x[1];
2724            true
2725        }
2726        fn eval_jac_g(
2727            &mut self,
2728            _x: Option<&[Number]>,
2729            _new_x: bool,
2730            mode: SparsityRequest<'_>,
2731        ) -> bool {
2732            match mode {
2733                SparsityRequest::Structure { irow, jcol } => {
2734                    irow.copy_from_slice(&[0, 0]);
2735                    jcol.copy_from_slice(&[0, 1]);
2736                }
2737                SparsityRequest::Values { values, .. } => {
2738                    values.copy_from_slice(&[1.0, 1.0]);
2739                }
2740            }
2741            true
2742        }
2743        fn eval_h(
2744            &mut self,
2745            _x: Option<&[Number]>,
2746            _new_x: bool,
2747            _obj_factor: Number,
2748            _lambda: Option<&[Number]>,
2749            _new_lambda: bool,
2750            mode: SparsityRequest<'_>,
2751        ) -> bool {
2752            match mode {
2753                SparsityRequest::Structure { irow, jcol } => {
2754                    irow.copy_from_slice(&[0, 1]);
2755                    jcol.copy_from_slice(&[0, 1]);
2756                }
2757                SparsityRequest::Values { values, .. } => {
2758                    values.copy_from_slice(&[1.0, 1.0]);
2759                }
2760            }
2761            true
2762        }
2763        fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
2764            *self.finalize_called.borrow_mut() = Some((sol.x.to_vec(), sol.obj_value));
2765        }
2766    }
2767
2768    #[test]
2769    fn application_sqp_path_solves_convex_eq_nlp_and_finalizes() {
2770        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2771        let tnlp = std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2772            finalize_called: std::rc::Rc::clone(&finalize_slot),
2773        }));
2774
2775        let mut app = IpoptApplication::new();
2776        app.initialize().unwrap();
2777        app.initialize_with_options_str("algorithm active-set-sqp\n")
2778            .unwrap();
2779        let status = app.optimize_tnlp(tnlp);
2780        assert_eq!(status, ApplicationReturnStatus::SolveSucceeded);
2781
2782        // The TNLP's finalize_solution must have been invoked.
2783        let recv = finalize_slot.borrow().clone();
2784        let (x_recv, obj_recv) = recv.expect("finalize_solution was not called");
2785        assert_eq!(x_recv.len(), 2);
2786        assert!((x_recv[0] - 0.0).abs() < 1e-6, "x[0] = {}", x_recv[0]);
2787        assert!((x_recv[1] - 1.0).abs() < 1e-6, "x[1] = {}", x_recv[1]);
2788        assert!(
2789            (obj_recv - (-1.5)).abs() < 1e-6,
2790            "obj = {} but expected -1.5",
2791            obj_recv
2792        );
2793    }
2794
2795    #[test]
2796    fn application_routes_to_sqp_case_insensitively() {
2797        let mut app = IpoptApplication::new();
2798        app.initialize().unwrap();
2799        app.initialize_with_options_str("algorithm Active-Set-SQP\n")
2800            .unwrap();
2801        // get_string_value may return the value as-stored (no
2802        // normalization); the dispatch must handle case
2803        // insensitively per the c11 design choice.
2804        assert!(app.is_sqp_algorithm_selected());
2805    }
2806
2807    #[test]
2808    fn application_constructs_and_loads_options() {
2809        let mut app = IpoptApplication::new();
2810        app.initialize().unwrap();
2811        // ipopt.opt-style file: an integer-typed option registered by
2812        // the Interfaces layer.
2813        app.initialize_with_options_str("print_level 5\nfile_print_level 7\n")
2814            .unwrap();
2815        let (level, found) = app.options().get_integer_value("print_level", "").unwrap();
2816        assert!(found);
2817        assert_eq!(level, 5);
2818    }
2819
2820    #[test]
2821    fn application_sqp_suboptions_propagate_to_builder() {
2822        // All SQP suboptions are read by algorithm_builder_snapshot
2823        // and baked into the builder's `sqp` field.
2824        let mut app = IpoptApplication::new();
2825        app.initialize().unwrap();
2826        app.initialize_with_options_str(
2827            "algorithm active-set-sqp\n\
2828             sqp_globalization l1-elastic\n\
2829             sqp_hessian lbfgs\n\
2830             sqp_max_iter 17\n\
2831             sqp_tol 1e-7\n\
2832             sqp_constr_viol_tol 1e-5\n\
2833             sqp_dual_inf_tol 1e-3\n\
2834             sqp_l1_penalty 2.5\n\
2835             sqp_bt_reduction 0.25\n\
2836             sqp_bt_min_alpha 1e-10\n\
2837             sqp_print_level 2\n\
2838             sqp_lbfgs_max_history 12\n",
2839        )
2840        .unwrap();
2841        let snap = app.algorithm_builder_snapshot();
2842        assert_eq!(
2843            snap.sqp.globalization,
2844            crate::sqp::SqpGlobalization::L1Elastic
2845        );
2846        assert_eq!(snap.sqp.hessian, crate::sqp::SqpHessianSource::Lbfgs);
2847        assert_eq!(snap.sqp.max_iter, 17);
2848        assert!((snap.sqp.tol - 1e-7).abs() < 1e-18);
2849        assert!((snap.sqp.constr_viol_tol - 1e-5).abs() < 1e-18);
2850        assert!((snap.sqp.dual_inf_tol - 1e-3).abs() < 1e-18);
2851        assert!((snap.sqp.l1_penalty - 2.5).abs() < 1e-18);
2852        assert!((snap.sqp.bt_reduction - 0.25).abs() < 1e-18);
2853        assert!((snap.sqp.bt_min_alpha - 1e-10).abs() < 1e-18);
2854        assert_eq!(snap.sqp.print_level, 2);
2855        assert_eq!(snap.sqp.lbfgs_max_history, 12);
2856    }
2857
2858    #[test]
2859    fn application_limited_memory_options_propagate_to_builder() {
2860        use crate::hess::lim_mem_quasi_newton::UpdateType;
2861
2862        // Default: no options set -> bit-exact with Ipopt's default
2863        // (bfgs, history 6). This is what the IPM path runs unless the
2864        // user opts in, so it must not drift.
2865        let mut app = IpoptApplication::new();
2866        app.initialize().unwrap();
2867        let def = app.algorithm_builder_from_options();
2868        assert_eq!(def.limited_memory_update_type, UpdateType::Bfgs);
2869        assert_eq!(def.limited_memory_max_history, 6);
2870
2871        // `limited_memory_update_type=sr1` and a custom history length
2872        // must reach the builder (these were registered upstream but
2873        // read nowhere on the IPM path before — see #131). Honoring
2874        // them is what lets SR1 break the monotone L-BFGS stall.
2875        let mut app = IpoptApplication::new();
2876        app.initialize().unwrap();
2877        app.initialize_with_options_str(
2878            "hessian_approximation limited-memory\n\
2879             limited_memory_update_type sr1\n\
2880             limited_memory_max_history 9\n",
2881        )
2882        .unwrap();
2883        let snap = app.algorithm_builder_from_options();
2884        assert_eq!(snap.limited_memory_update_type, UpdateType::Sr1);
2885        assert_eq!(snap.limited_memory_max_history, 9);
2886    }
2887
2888    #[test]
2889    fn application_sqp_warm_start_round_trip() {
2890        // Drive the convex-equality TNLP through the SQP path
2891        // twice. The first solve produces a working set; the
2892        // second is warm-started from it. The second must converge
2893        // with zero QP solves (the first KKT check declares
2894        // optimality immediately).
2895        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2896        let tnlp_rc: std::rc::Rc<std::cell::RefCell<dyn TNLP>> =
2897            std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2898                finalize_called: std::rc::Rc::clone(&finalize_slot),
2899            }));
2900
2901        let mut app = IpoptApplication::new();
2902        app.initialize().unwrap();
2903        app.initialize_with_options_str("algorithm active-set-sqp\n")
2904            .unwrap();
2905
2906        // Cold solve.
2907        let status_a = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2908        assert_eq!(status_a, ApplicationReturnStatus::SolveSucceeded);
2909        let ws = app.last_sqp_working_set().cloned();
2910        assert!(ws.is_some(), "cold solve must yield a working set");
2911
2912        // Build the warm-start iterate from the converged finalize
2913        // payload (just x; pad multipliers to 0 since the test
2914        // problem is convex).
2915        let (x_recv, _) = finalize_slot.borrow().clone().unwrap();
2916        let warm = crate::sqp::SqpIterates {
2917            x: x_recv,
2918            lambda_g: vec![1.0],
2919            lambda_x: vec![0.0, 0.0],
2920            working: ws,
2921        };
2922        app.set_sqp_warm_start(warm);
2923
2924        // Warm solve.
2925        let status_b = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2926        assert_eq!(status_b, ApplicationReturnStatus::SolveSucceeded);
2927        assert!(app.last_sqp_working_set().is_some());
2928    }
2929
2930    #[test]
2931    fn application_sqp_warm_start_auto_clears_after_use() {
2932        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2933        let tnlp_rc: std::rc::Rc<std::cell::RefCell<dyn TNLP>> =
2934            std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2935                finalize_called: std::rc::Rc::clone(&finalize_slot),
2936            }));
2937        let mut app = IpoptApplication::new();
2938        app.initialize().unwrap();
2939        app.initialize_with_options_str("algorithm active-set-sqp\n")
2940            .unwrap();
2941        app.set_sqp_warm_start(crate::sqp::SqpIterates {
2942            x: vec![0.0, 1.0],
2943            lambda_g: vec![1.0],
2944            lambda_x: vec![0.0, 0.0],
2945            working: None,
2946        });
2947        assert!(app.sqp_warm_start.is_some());
2948        let _ = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2949        assert!(
2950            app.sqp_warm_start.is_none(),
2951            "warm-start input must be auto-cleared after use"
2952        );
2953    }
2954
2955    #[test]
2956    fn application_sqp_suboptions_default_when_unset() {
2957        // Without any sqp_* settings, the snapshot should equal
2958        // SqpOptions::default().
2959        let mut app = IpoptApplication::new();
2960        app.initialize().unwrap();
2961        let snap = app.algorithm_builder_snapshot();
2962        let d = crate::sqp::SqpOptions::default();
2963        assert_eq!(snap.sqp.globalization, d.globalization);
2964        assert_eq!(snap.sqp.hessian, d.hessian);
2965        assert_eq!(snap.sqp.max_iter, d.max_iter);
2966        assert!((snap.sqp.tol - d.tol).abs() < 1e-18);
2967        assert!((snap.sqp.constr_viol_tol - d.constr_viol_tol).abs() < 1e-18);
2968        assert!((snap.sqp.dual_inf_tol - d.dual_inf_tol).abs() < 1e-18);
2969        assert!((snap.sqp.l1_penalty - d.l1_penalty).abs() < 1e-18);
2970        assert!((snap.sqp.bt_reduction - d.bt_reduction).abs() < 1e-18);
2971        assert!((snap.sqp.bt_min_alpha - d.bt_min_alpha).abs() < 1e-18);
2972        assert_eq!(snap.sqp.print_level, d.print_level);
2973        assert_eq!(snap.sqp.lbfgs_max_history, d.lbfgs_max_history);
2974    }
2975
2976    #[test]
2977    fn application_reports_problem_dimensions() {
2978        let app = IpoptApplication::new();
2979        let mut tnlp = Hs071Stub;
2980        let info = app.problem_dimensions(&mut tnlp).unwrap();
2981        assert_eq!(info.n, 4);
2982        assert_eq!(info.m, 2);
2983        assert_eq!(info.nnz_jac_g, 8);
2984        assert_eq!(info.nnz_h_lag, 10);
2985    }
2986}