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        }
719        let (app_status, solver_status) = match res.status {
720            crate::sqp::SqpStatus::Optimal => (
721                ApplicationReturnStatus::SolveSucceeded,
722                pounce_nlp::SolverReturn::Success,
723            ),
724            crate::sqp::SqpStatus::MaxIter => (
725                ApplicationReturnStatus::MaximumIterationsExceeded,
726                pounce_nlp::SolverReturn::MaxiterExceeded,
727            ),
728            crate::sqp::SqpStatus::InfeasibleSubproblem => (
729                ApplicationReturnStatus::InfeasibleProblemDetected,
730                pounce_nlp::SolverReturn::LocalInfeasibility,
731            ),
732            crate::sqp::SqpStatus::LineSearchFailed => (
733                ApplicationReturnStatus::SearchDirectionBecomesTooSmall,
734                pounce_nlp::SolverReturn::ErrorInStepComputation,
735            ),
736        };
737
738        // Forward to the user TNLP's finalize_solution. We pass
739        // the SQP iterate and recovered multipliers via the
740        // OrigIpoptNlp's lifting hooks. Failure here is silent
741        // (we still return the algorithm's status) — the user
742        // sees the right ApplicationReturnStatus regardless.
743        let _ = finalize_via_sqp(&nlp_rc, &res, solver_status, &tnlp);
744
745        app_status
746    }
747
748    /// Build a *copy* of the algorithm builder configured per the
749    /// current options. The SQP path uses this so it gets a
750    /// fresh builder without mutating the application's state.
751    fn algorithm_builder_snapshot(&self) -> AlgorithmBuilder {
752        let mut builder = AlgorithmBuilder::default();
753        apply_sqp_options(&self.options, &mut builder.sqp);
754        apply_qp_subproblem_options(&self.options, &mut builder.sqp_qp);
755        builder
756    }
757
758    /// Construct a LinearBackendFactory honoring the
759    /// `linear_solver` option. Default FERAL; HSL MA57 when
760    /// built with the `ma57` feature.
761    fn make_backend_factory(&self) -> LinearBackendFactory {
762        Box::new(
763            |_choice| -> Box<dyn pounce_linsol::SparseSymLinearSolverInterface> {
764                Box::new(pounce_feral::FeralSolverInterface::new())
765            },
766        )
767    }
768
769    /// Phase 3.5 auto-fallback driver.
770    ///
771    /// Runs the standard solve (no wrapper) first. If it ends in a
772    /// trigger-class status (`Restoration_Failed`, `Infeasible_Problem_Detected`,
773    /// `Solved_To_Acceptable_Level`, `Maximum_Iterations_Exceeded`, or
774    /// `Not_Enough_Degrees_Of_Freedom`), retries transparently with
775    /// the ℓ₁ wrapper enabled. Promotes the retry's status only if
776    /// it returns `Solve_Succeeded`; otherwise returns the original
777    /// status.
778    ///
779    /// Caveat: the user TNLP's `finalize_solution` runs once per
780    /// attempt. When the retry doesn't promote, the user's captured
781    /// fields hold the retry's iterate (the ℓ₁-best least-infeasible
782    /// point) even though the returned status is the original's.
783    /// Documented on the option's help text; tightening this is a
784    /// Phase-4 follow-up.
785    fn run_with_l1_fallback(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
786        // First attempt: the standard IPM solve, no ℓ₁ wrapper. Only
787        // reached for `m > 0`, so `optimize_constrained` is exact.
788        let first_status = self.optimize_constrained(Rc::clone(&tnlp));
789        if !is_l1_fallback_trigger(first_status) {
790            return first_status;
791        }
792        // Trigger fired. Flip the wrapper option for the retry and
793        // restore it after — keeps the user's option-table view of the
794        // session exactly as they left it.
795        let prev = self
796            .options
797            .get_string_value("l1_exact_penalty_barrier", "")
798            .ok();
799        let _ = self
800            .options
801            .set_string_value("l1_exact_penalty_barrier", "yes", true, false);
802        let retry_status = self
803            .run_l1_penalty_outer_loop(Rc::clone(&tnlp))
804            .unwrap_or(ApplicationReturnStatus::InternalError);
805        let _ = self.options.set_string_value(
806            "l1_exact_penalty_barrier",
807            prev.as_ref().map(|(v, _)| v.as_str()).unwrap_or("no"),
808            true,
809            false,
810        );
811        if matches!(retry_status, ApplicationReturnStatus::SolveSucceeded) {
812            retry_status
813        } else {
814            first_status
815        }
816    }
817
818    /// μ-strategy auto-fallback driver (pounce#138).
819    ///
820    /// Runs the standard solve first. If it stalls short of optimal in a
821    /// way a μ-strategy flip can plausibly fix — `Solved_To_Acceptable_Level`
822    /// or `Maximum_Iterations_Exceeded`, the two signatures seen on the
823    /// princetonlib instances where the dual infeasibility parks above
824    /// `tol` while constraint violation and complementarity are already
825    /// deeply converged — it flips `mu_strategy` (adaptive↔monotone) and
826    /// solves once more. The retry's status is promoted only if it returns
827    /// `Solve_Succeeded`; otherwise the original status is returned.
828    ///
829    /// (maxcut/price stall at acceptable-level under adaptive; fermat2_vareps
830    /// stalls at `max_iter` — hence both triggers. flosp2tm is μ-independent
831    /// and correctly does not promote.)
832    ///
833    /// The flip direction is taken from the option's current value:
834    /// `adaptive` → `monotone`, anything else (including absent, which
835    /// the builder treats as monotone) → `adaptive`. The option table is
836    /// restored to the user's original view afterward.
837    ///
838    /// Caveat (shared with the ℓ₁ fallback): the user TNLP's
839    /// `finalize_solution` runs once per attempt, so when the retry
840    /// doesn't promote the captured fields hold the retry's iterate.
841    fn run_with_mu_strategy_fallback(
842        &mut self,
843        tnlp: Rc<RefCell<dyn TNLP>>,
844    ) -> ApplicationReturnStatus {
845        let first_status = self.optimize_constrained(Rc::clone(&tnlp));
846        if !matches!(
847            first_status,
848            ApplicationReturnStatus::SolvedToAcceptableLevel
849                | ApplicationReturnStatus::MaximumIterationsExceeded
850        ) {
851            return first_status;
852        }
853        // Flip the strategy for one retry. The parser maps "adaptive" →
854        // Adaptive and every other value (incl. unset) → Monotone, so the
855        // opposite of an explicit "adaptive" is "monotone" and the
856        // opposite of anything else is "adaptive".
857        let prev = self.options.get_string_value("mu_strategy", "").ok();
858        let was_adaptive = prev
859            .as_ref()
860            .map(|(v, found)| *found && v == "adaptive")
861            .unwrap_or(false);
862        let flipped = if was_adaptive { "monotone" } else { "adaptive" };
863        let _ = self
864            .options
865            .set_string_value("mu_strategy", flipped, true, false);
866        let retry_status = self.optimize_constrained(Rc::clone(&tnlp));
867        // Restore the user's original option-table view.
868        let _ = self.options.set_string_value(
869            "mu_strategy",
870            prev.as_ref()
871                .filter(|(_, found)| *found)
872                .map(|(v, _)| v.as_str())
873                .unwrap_or("monotone"),
874            true,
875            false,
876        );
877        if matches!(retry_status, ApplicationReturnStatus::SolveSucceeded) {
878            retry_status
879        } else {
880            first_status
881        }
882    }
883
884    /// Phase-3 ℓ₁-exact penalty-barrier outer loop.
885    ///
886    /// Builds an [`L1PenaltyBarrierTnlp`] wrapper around the user
887    /// TNLP, runs the constrained IPM at the current ρ, escalates ρ
888    /// per Byrd-Nocedal-Waltz steering, and terminates on any of:
889    ///   - slack sum collapses (`Σ(p+n) ≤ l1_slack_tol`)
890    ///   - inner solve returns non-Optimal (escalation won't fix
891    ///     numerical / restoration failure at this ρ)
892    ///   - ρ already at `l1_penalty_max`
893    ///   - `l1_penalty_max_outer_iter` reached
894    ///
895    /// After the loop, if the inner status is `SolveSucceeded` or
896    /// `SolvedToAcceptableLevel` but slacks didn't collapse, override
897    /// to `Infeasible_Problem_Detected` — the returned point is the
898    /// ℓ₁-best least-infeasible iterate, which is informative even
899    /// though the original constraints are not satisfied.
900    ///
901    /// Returns `Some(status)` if the wrapper ran the solve, `None` if
902    /// wrapper construction failed (caller should fall through to the
903    /// standard dispatch path).
904    fn run_l1_penalty_outer_loop(
905        &mut self,
906        tnlp: Rc<RefCell<dyn TNLP>>,
907    ) -> Option<ApplicationReturnStatus> {
908        let rho_init = self.l1_penalty_init();
909        let rho_max = self.l1_penalty_max().max(rho_init);
910        let factor = self.l1_penalty_increase_factor().max(1.0);
911        let tau = self.l1_steering_factor();
912        let slack_tol = self.l1_slack_tol();
913        let max_outer = self.l1_penalty_max_outer_iter().max(1);
914
915        let mut wrapper = pounce_l1penalty::L1PenaltyBarrierTnlp::new(Rc::clone(&tnlp), rho_init)?;
916        if wrapper.m_eq() == 0 {
917            // Nothing to slack — let the standard dispatch path handle
918            // this TNLP unmodified.
919            return None;
920        }
921        wrapper.set_defer_inner_finalize(true);
922        let wrapper_rc = Rc::new(RefCell::new(wrapper));
923
924        let mut rho = rho_init;
925        let mut last_status = ApplicationReturnStatus::InternalError;
926        for _outer in 0..max_outer {
927            wrapper_rc.borrow_mut().set_rho(rho);
928            let dyn_tnlp: Rc<RefCell<dyn TNLP>> = wrapper_rc.clone();
929            last_status = self.optimize_constrained(dyn_tnlp);
930
931            let w = wrapper_rc.borrow();
932            if !w.has_solution() {
933                // Inner solve aborted before producing an iterate.
934                drop(w);
935                break;
936            }
937            let slack_sum = w.last_slack_sum();
938            let y_eq_inf = w.last_y_eq_inf_norm();
939            drop(w);
940
941            // Termination decisions.
942            let inner_ok = matches!(
943                last_status,
944                ApplicationReturnStatus::SolveSucceeded
945                    | ApplicationReturnStatus::SolvedToAcceptableLevel
946            );
947            if !inner_ok {
948                break;
949            }
950            if slack_sum.is_finite() && slack_sum <= slack_tol {
951                break;
952            }
953            if rho >= rho_max {
954                break;
955            }
956            // BNW steering: ρ_new = max(ρ·factor, τ·‖y_eq‖∞ + ε)
957            let geom = rho * factor;
958            let steer = tau * y_eq_inf + 1.0e-12;
959            rho = geom.max(steer).min(rho_max);
960        }
961
962        // Forward to the user's inner.finalize_solution exactly once.
963        let w = wrapper_rc.borrow();
964        if w.has_solution() {
965            let x_trunc: Vec<Number> = w.last_x_trunc().to_vec();
966            let lambda: Vec<Number> = w.last_lambda().to_vec();
967            let z_l: Vec<Number> = w.last_z_l_trunc().to_vec();
968            let z_u: Vec<Number> = w.last_z_u_trunc().to_vec();
969            let solver_status = w.last_status().unwrap_or(SolverReturn::InternalError);
970            let slack_sum = w.last_slack_sum();
971            drop(w);
972
973            // Honest-infeasibility upgrade (Phase 3): if the inner
974            // solve says SolveSucceeded / SolvedToAcceptableLevel but
975            // the slacks did not collapse, the original problem is
976            // locally infeasible at the returned point. Override the
977            // application status; the user-visible Solution.status is
978            // updated below to the matching SolverReturn so the inner
979            // TNLP sees a consistent picture.
980            let infeasible_certificate = matches!(
981                last_status,
982                ApplicationReturnStatus::SolveSucceeded
983                    | ApplicationReturnStatus::SolvedToAcceptableLevel
984            ) && slack_sum.is_finite()
985                && slack_sum > slack_tol;
986            let final_app_status = if infeasible_certificate {
987                ApplicationReturnStatus::InfeasibleProblemDetected
988            } else {
989                last_status
990            };
991            let final_solver_status = if infeasible_certificate {
992                SolverReturn::LocalInfeasibility
993            } else {
994                solver_status
995            };
996
997            // Recompute f(x*) and c(x*) on the inner.
998            let f_inner = tnlp
999                .borrow_mut()
1000                .eval_f(&x_trunc, true)
1001                .unwrap_or(Number::NAN);
1002            let m = tnlp
1003                .borrow_mut()
1004                .get_nlp_info()
1005                .map(|i| i.m as usize)
1006                .unwrap_or(0);
1007            let mut g_inner = vec![0.0; m];
1008            if m > 0 {
1009                let _ = tnlp.borrow_mut().eval_g(&x_trunc, false, &mut g_inner);
1010            }
1011            tnlp.borrow_mut().finalize_solution(
1012                Solution {
1013                    status: final_solver_status,
1014                    x: &x_trunc,
1015                    z_l: &z_l,
1016                    z_u: &z_u,
1017                    g: &g_inner,
1018                    lambda: &lambda,
1019                    obj_value: f_inner,
1020                },
1021                &TnlpIpoptData::default(),
1022                &TnlpIpoptCq::default(),
1023            );
1024            return Some(final_app_status);
1025        }
1026        // No solution captured at all — pass the inner status through.
1027        Some(last_status)
1028    }
1029
1030    /// Constrained-NLP path: build adapter → OrigIpoptNlp → algorithm
1031    /// bundle, run `optimize`, populate statistics, and call
1032    /// `finalize_solution` on the user's TNLP.
1033    fn optimize_constrained(&mut self, tnlp: Rc<RefCell<dyn TNLP>>) -> ApplicationReturnStatus {
1034        let t_start = Instant::now();
1035
1036        // `print_user_options yes` — dump the OptionsList before the
1037        // solve. Mirrors `IpoptApplication::call_optimize` (upstream
1038        // calls `Jnlst().Printf(.., "%s", options_->PrintUserOptions())`).
1039        let print_opts = self
1040            .options
1041            .get_bool_value("print_user_options", "")
1042            .ok()
1043            .and_then(|(v, f)| f.then_some(v))
1044            .unwrap_or(false);
1045        if print_opts {
1046            print!(
1047                "\nList of user-set options:\n\n{}",
1048                self.options.print_user_options()
1049            );
1050        }
1051
1052        // `print_options_documentation yes` — dump the full registry
1053        // (every option with type, default, valid range/strings, and
1054        // long description) before the solve. Honors
1055        // `print_options_mode` (`text` / `latex` / `doxygen`; only
1056        // `text` is implemented today, the others fall through with a
1057        // one-line note) and `print_advanced_options`. Mirrors
1058        // upstream `IpoptApplication::call_optimize`'s
1059        // `print_options_documentation` branch and `Common/IpRegOptions.cpp`
1060        // `OutputOptionDocumentation`.
1061        let print_doc = self
1062            .options
1063            .get_bool_value("print_options_documentation", "")
1064            .ok()
1065            .and_then(|(v, f)| f.then_some(v))
1066            .unwrap_or(false);
1067        if print_doc {
1068            let mode = self
1069                .options
1070                .get_string_value("print_options_mode", "")
1071                .ok()
1072                .map(|(v, _)| PrintOptionsMode::from_tag(&v))
1073                .unwrap_or(PrintOptionsMode::Text);
1074            let advanced = self
1075                .options
1076                .get_bool_value("print_advanced_options", "")
1077                .ok()
1078                .map(|(v, _)| v)
1079                .unwrap_or(false);
1080            print!(
1081                "\n# Pounce options registry\n\n{}",
1082                self.reg_options.print_options_documentation(mode, advanced)
1083            );
1084        }
1085
1086        // Mint a fresh `TimingStatistics` for this solve — shared (via
1087        // `Rc`) with the data and the NLP below so every `eval_*` and
1088        // every iterate-phase records into the same accumulator. The
1089        // application keeps its own `Rc` so callers can read totals out
1090        // via [`Self::timing_stats`].
1091        let timing = Rc::new(TimingStatistics::new());
1092        *self.timing.borrow_mut() = Rc::clone(&timing);
1093        timing.overall_alg.start();
1094
1095        // Reset the linear-solver summary sink so back-to-back solves
1096        // don't bleed factor counters / extremal pivots into each
1097        // other. Surviving the lock failure with a debug-assert keeps
1098        // a poisoned mutex from sinking a release build that doesn't
1099        // even consume the summary.
1100        if let Ok(mut guard) = self.linsol_summary_sink.lock() {
1101            *guard = LinearSolverSummary::default();
1102        } else {
1103            debug_assert!(false, "linsol summary sink mutex poisoned");
1104        }
1105
1106        // Build adapter + Nlp. Honor `fixed_variable_treatment` (default
1107        // `make_parameter`; pounce additionally implements `relax_bounds`,
1108        // which the adapter also auto-selects as a fallback when
1109        // `make_parameter` would leave `n_x_var < n_c` — mirrors upstream
1110        // `IpTNLPAdapter.cpp:623-633`).
1111        let lo_inf = self
1112            .options
1113            .get_numeric_value("nlp_lower_bound_inf", "")
1114            .ok()
1115            .and_then(|(v, f)| f.then_some(v))
1116            .unwrap_or(DEFAULT_NLP_LOWER_BOUND_INF);
1117        let up_inf = self
1118            .options
1119            .get_numeric_value("nlp_upper_bound_inf", "")
1120            .ok()
1121            .and_then(|(v, f)| f.then_some(v))
1122            .unwrap_or(DEFAULT_NLP_UPPER_BOUND_INF);
1123        let fixed_treatment = match self
1124            .options
1125            .get_string_value("fixed_variable_treatment", "")
1126            .ok()
1127            .and_then(|(v, f)| f.then_some(v))
1128            .as_deref()
1129        {
1130            Some("relax_bounds") => FixedVarTreatment::RelaxBounds,
1131            // `make_constraint` / `make_parameter_nodual` not yet
1132            // implemented; fall back to `make_parameter` (auto-retry to
1133            // `relax_bounds` will still kick in if DOF runs short).
1134            _ => FixedVarTreatment::MakeParameter,
1135        };
1136        let adapter = match TNLPAdapter::new_with_options(
1137            Rc::clone(&tnlp),
1138            lo_inf,
1139            up_inf,
1140            fixed_treatment,
1141        ) {
1142            Ok(a) => Rc::new(RefCell::new(a)),
1143            Err(_) => {
1144                timing.overall_alg.end();
1145                return ApplicationReturnStatus::InvalidProblemDefinition;
1146            }
1147        };
1148        // Carry the user's constant `obj_scaling_factor` (default 1.0;
1149        // negative ⇒ maximize) into the NLP. Until pounce#128's
1150        // follow-up this option was registered but never read, so it
1151        // was silently a no-op — maximization diverged because the
1152        // algorithm minimized the unscaled objective.
1153        let obj_scaling_factor = self
1154            .options
1155            .get_numeric_value("obj_scaling_factor", "")
1156            .ok()
1157            .and_then(|(v, f)| f.then_some(v))
1158            .unwrap_or(1.0);
1159        let mut orig_nlp = match OrigIpoptNlp::new(
1160            Rc::clone(&adapter),
1161            Rc::new(ConstObjScaling(obj_scaling_factor)),
1162        ) {
1163            Ok(n) => n,
1164            Err(_) => {
1165                timing.overall_alg.end();
1166                return ApplicationReturnStatus::InternalError;
1167            }
1168        };
1169        orig_nlp.set_timing_stats(Rc::clone(&timing));
1170
1171        // Mirror upstream `OrigIpoptNLP::InitializeStructures` (IpOrigIpoptNLP.cpp:299):
1172        // bail out with NotEnoughDegreesOfFreedom when there are fewer free
1173        // variables than equality constraints. Without this gate, square /
1174        // over-determined systems push the algorithm into restoration on
1175        // iter 0 and exit Restoration_Failed instead of the cleaner DOF code.
1176        let n_x_var = orig_nlp.x_space().dim();
1177        let n_c = orig_nlp.c_space().dim();
1178        if n_x_var > 0 && n_x_var < n_c {
1179            timing.overall_alg.end();
1180            return ApplicationReturnStatus::NotEnoughDegreesOfFreedom;
1181        }
1182
1183        // Relax `x_L / x_U / d_L / d_U` by `bound_relax_factor` (default
1184        // 1e-8), capped by `constr_viol_tol` (default 1e-4). Matches
1185        // `OrigIpoptNLP::InitializeStructures` lines 343-358.
1186        let bound_relax_factor = self
1187            .options
1188            .get_numeric_value("bound_relax_factor", "")
1189            .ok()
1190            .and_then(|(v, f)| f.then_some(v))
1191            .unwrap_or(1e-8);
1192        let constr_viol_tol = self
1193            .options
1194            .get_numeric_value("constr_viol_tol", "")
1195            .ok()
1196            .and_then(|(v, f)| f.then_some(v))
1197            .unwrap_or(1e-4);
1198        orig_nlp.relax_bounds(bound_relax_factor, constr_viol_tol);
1199
1200        // Apply automatic NLP scaling per `nlp_scaling_method` option
1201        // (port of `OrigIpoptNLP::InitializeStructures` →
1202        // `NLPScalingObject::DetermineScaling`). Default is
1203        // `gradient-based` to match upstream Ipopt 3.14.
1204        let scaling_method = self
1205            .options
1206            .get_string_value("nlp_scaling_method", "")
1207            .ok()
1208            .and_then(|(v, f)| f.then_some(v))
1209            .unwrap_or_else(|| "gradient-based".to_string());
1210        let scaling_method = match scaling_method.as_str() {
1211            "none" => ScalingMethod::None,
1212            "gradient-based" => ScalingMethod::GradientBased,
1213            "user-scaling" => ScalingMethod::UserScaling,
1214            // `equilibration-based` is registered upstream but not yet
1215            // implemented in pounce; fall back to gradient-based (the
1216            // upstream default) to keep behavior predictable.
1217            _ => ScalingMethod::GradientBased,
1218        };
1219        let max_gradient = self
1220            .options
1221            .get_numeric_value("nlp_scaling_max_gradient", "")
1222            .ok()
1223            .and_then(|(v, f)| f.then_some(v))
1224            .unwrap_or(100.0);
1225        let min_value = self
1226            .options
1227            .get_numeric_value("nlp_scaling_min_value", "")
1228            .ok()
1229            .and_then(|(v, f)| f.then_some(v))
1230            .unwrap_or(1e-8);
1231        let obj_target_gradient = self
1232            .options
1233            .get_numeric_value("nlp_scaling_obj_target_gradient", "")
1234            .ok()
1235            .and_then(|(v, f)| f.then_some(v))
1236            .unwrap_or(0.0);
1237        let constr_target_gradient = self
1238            .options
1239            .get_numeric_value("nlp_scaling_constr_target_gradient", "")
1240            .ok()
1241            .and_then(|(v, f)| f.then_some(v))
1242            .unwrap_or(0.0);
1243        orig_nlp.determine_scaling_from_starting_point(
1244            scaling_method,
1245            max_gradient,
1246            min_value,
1247            obj_target_gradient,
1248            constr_target_gradient,
1249        );
1250
1251        let nlp_handle: Rc<RefCell<dyn IpoptNlp>> = Rc::new(RefCell::new(orig_nlp));
1252
1253        // Build the algorithm strategy bundle. Read coarse knobs from
1254        // the OptionsList where we have them; fall through to defaults
1255        // otherwise. The full upstream parsing surface (mu_strategy,
1256        // hessian_approximation, line_search_method, ...) is wired by
1257        // `AlgBuilder::RegisterOptions` in upstream — that registry
1258        // hookup lands as a follow-up; default builder is correct for
1259        // HS71-class problems.
1260        let builder = self.algorithm_builder_from_options();
1261
1262        // Linear-solver backend. The default factory is option-aware
1263        // — it reads the `feral_*` extension options off the same
1264        // `OptionsList` that drove the IPM-level builder above so
1265        // per-problem `.opt` files can flip backend knobs without
1266        // rebuilding pounce.
1267        let feral_cfg = feral_config_from_options(&self.options);
1268        let factory = self.linear_backend_factory.take().unwrap_or_else(|| {
1269            default_backend_factory_with_sink(feral_cfg, Arc::clone(&self.linsol_summary_sink))
1270        });
1271        let bundle = builder.build_with_backend(factory);
1272
1273        // Wire the data / cq pair around the NLP. Install the shared
1274        // `TimingStatistics` so the algorithm's iterate phases
1275        // (output, convergence, hessian, μ, search-direction,
1276        // line-search, accept) all record into the same accumulator
1277        // the application exposes via `timing_stats()`.
1278        let data: crate::ipopt_data::IpoptDataHandle = Rc::new(RefCell::new(AlgIpoptData::new()));
1279        data.borrow_mut().timing = Rc::clone(&timing);
1280        let cq: crate::ipopt_cq::IpoptCqHandle = Rc::new(RefCell::new(
1281            IpoptCalculatedQuantities::new(Rc::clone(&data), Rc::clone(&nlp_handle)),
1282        ));
1283        // Correction size for very small slacks (default mach_eps^{3/4});
1284        // drives the safe-slack bound-adjustment mechanism.
1285        if let Ok((v, true)) = self.options.get_numeric_value("slack_move", "") {
1286            cq.borrow_mut().slack_move = v;
1287        }
1288
1289        // Seed `data.curr` with a zero-valued iterate of the correct
1290        // dimensions. The `IterateInitializer` consumes these as its
1291        // template (it overwrites `x`, `s`, multipliers in place); we
1292        // just need the dim metadata.
1293        {
1294            let nlp_borrow = nlp_handle.borrow();
1295            let n_x = nlp_borrow.n();
1296            let n_s = nlp_borrow.m_ineq();
1297            let n_yc = nlp_borrow.m_eq();
1298            let n_yd = nlp_borrow.m_ineq();
1299            let n_zl = nlp_borrow.x_l().dim();
1300            let n_zu = nlp_borrow.x_u().dim();
1301            let n_vl = nlp_borrow.d_l().dim();
1302            let n_vu = nlp_borrow.d_u().dim();
1303            drop(nlp_borrow);
1304            let iv = IteratesVector::new(
1305                Rc::new(DenseVectorSpace::new(n_x).make_new_dense()),
1306                Rc::new(DenseVectorSpace::new(n_s).make_new_dense()),
1307                Rc::new(DenseVectorSpace::new(n_yc).make_new_dense()),
1308                Rc::new(DenseVectorSpace::new(n_yd).make_new_dense()),
1309                Rc::new(DenseVectorSpace::new(n_zl).make_new_dense()),
1310                Rc::new(DenseVectorSpace::new(n_zu).make_new_dense()),
1311                Rc::new(DenseVectorSpace::new(n_vl).make_new_dense()),
1312                Rc::new(DenseVectorSpace::new(n_vu).make_new_dense()),
1313            );
1314            data.borrow_mut().set_curr(iv);
1315        }
1316
1317        // Full primal-dual warm restart (debugger `resolve`): if a
1318        // captured iterate is queued, install it onto `data.curr` over
1319        // the placeholder so the `WarmStartIterateInitializer`'s
1320        // re-optimize branch (x already initialized) keeps it and only
1321        // clamps multipliers / sets target_mu — no cold re-seed from the
1322        // NLP. Skipped (with a warning) if the dimensions don't line up,
1323        // e.g. an option changed the problem structure between solves.
1324        if let Some(snap) = self.warm_start_iterate.take() {
1325            let dims_match = {
1326                let borrow = data.borrow();
1327                borrow
1328                    .curr
1329                    .as_ref()
1330                    .map(|c| iterates_dims(c) == iterates_dims(snap.iterates()))
1331                    .unwrap_or(false)
1332            };
1333            if dims_match {
1334                data.borrow_mut().set_curr(snap.iterates().clone());
1335                data.borrow_mut().curr_mu = snap.mu();
1336            } else {
1337                tracing::warn!(
1338                    target: "pounce::warm_start",
1339                    "debugger warm-restart iterate dimensions differ from the fresh \
1340                     solve; ignoring the captured iterate and seeding normally"
1341                );
1342            }
1343        }
1344
1345        let max_iter = self
1346            .options
1347            .get_integer_value("max_iter", "")
1348            .ok()
1349            .and_then(|(v, f)| f.then_some(v))
1350            .unwrap_or(3000);
1351        let tol = self
1352            .options
1353            .get_numeric_value("tol", "")
1354            .ok()
1355            .and_then(|(v, f)| f.then_some(v))
1356            .unwrap_or(1e-8);
1357        data.borrow_mut().tol = tol;
1358
1359        let mut alg = IpoptAlgorithm::new(data, cq, bundle)
1360            .with_nlp(Rc::clone(&nlp_handle))
1361            .with_tnlp(Rc::clone(&tnlp));
1362        // Mint a fresh restoration factory per inner solve if a
1363        // provider is configured (pounce#10 Phase 3). Falls back to
1364        // the legacy one-shot `restoration_factory` slot when no
1365        // provider is set, preserving single-shot caller behavior.
1366        if let Some(provider) = self.restoration_factory_provider.as_mut() {
1367            self.restoration_factory = Some(provider());
1368        }
1369        if let Some(factory) = self.restoration_factory.as_mut() {
1370            alg = alg.with_restoration(factory());
1371        }
1372        if let Some(diag) = self.diagnostics.as_ref() {
1373            alg = alg.with_diagnostics(Rc::clone(diag));
1374        }
1375        // Move the interactive debugger hook (if any) into the main
1376        // algorithm. Taken — not cloned — so it drives exactly this
1377        // solve; a subsequent solve must reinstall it.
1378        if let Some(hook) = self.debug_hook.take() {
1379            alg = alg.with_debug_hook(hook);
1380        }
1381        alg.max_iter = max_iter;
1382        // Honor `print_level == 0`: suppress the per-iteration table
1383        // that the algorithm writes straight to stdout. (The Phase-7
1384        // journalist surface respects `print_level` already; this is
1385        // the legacy direct-print site that needs the same gate.)
1386        if let Ok((v, found)) = self.options.get_integer_value("print_level", "") {
1387            if found && v <= 0 {
1388                alg.print_iter_output = false;
1389                // The nested restoration IPM is built inside the
1390                // restoration driver, not by `IpoptAlgorithm::new`, so
1391                // it never sees this gate unless we forward it.
1392                if let Some(resto) = alg.restoration.as_mut() {
1393                    resto.set_print_iter_output(false);
1394                }
1395            }
1396        }
1397
1398        // Per-iteration history (pounce#71): when requested, capture the
1399        // `pounce::iteration` events emitted during the solve into an
1400        // `IterRecord` trajectory via the observability collector layer.
1401        // This replaces the old in-loop `iter_history` accumulation; it
1402        // requires the collector to be installed in the active
1403        // subscriber (the CLI / Python / C frontends install it via
1404        // `pounce_observability::init_subscriber`; tests call
1405        // `init_for_tests`). The collector scopes out restoration
1406        // sub-solve iterations via the `restoration` span, so the
1407        // trajectory matches the previous behavior (outer iters only).
1408        let iter_capture = self
1409            .record_iter_history
1410            .then(pounce_observability::IterCaptureGuard::start);
1411
1412        let solver_status = alg.optimize();
1413
1414        let captured_iters = iter_capture.map(|g| g.finish()).unwrap_or_default();
1415        // Close the overall-algorithm timer on the success path. The
1416        // early-return arms above end it themselves before bailing out;
1417        // this one matches upstream `IpoptApplication::call_optimize`
1418        // (which calls `EndCpuTime()` on overall_alg right after
1419        // `Optimize` returns, regardless of solver_status).
1420        timing.overall_alg.end();
1421
1422        // Drain counters / iter count off the algorithm.
1423        {
1424            let mut stats = self.statistics.borrow_mut();
1425            {
1426                let d = alg.data.borrow();
1427                stats.iteration_count = d.iter_count;
1428                // Converged barrier parameter μ — threaded forward into a
1429                // warm-started corrector's `mu_init` / `warm_start_target_mu`
1430                // for predictor–corrector path following (pounce#86).
1431                stats.final_mu = d.curr_mu;
1432            }
1433            stats.total_wallclock_time_secs = t_start.elapsed().as_secs_f64();
1434            // Restoration-phase audit counters (pounce#12). Zero on
1435            // problems where restoration never fires; populated by
1436            // `IpoptAlgorithm::invoke_restoration`.
1437            stats.restoration_calls = alg.resto_calls;
1438            stats.restoration_inner_iters = alg.resto_inner_iters;
1439            stats.restoration_outer_iters = alg.resto_outer_iters;
1440            stats.restoration_wall_secs = alg.resto_wall_secs;
1441            stats.iterations = captured_iters;
1442            // Capture the final *scaled* objective at the algorithm's
1443            // (compressed `x_var`-space) iterate via the NLP: the
1444            // algorithm-side `eval_f` returns `f * obj_scale_factor`.
1445            // `final_objective` is seeded with it only as a best-effort
1446            // fallback; the success path below overwrites it with the
1447            // true unscaled objective from `finalize_via_orig_nlp`
1448            // (which evaluates the user TNLP directly).
1449            let curr_x = alg.data.borrow().curr.as_ref().map(|c| c.x.clone());
1450            if let Some(x) = curr_x {
1451                if let Ok(f) = try_eval_curr_f(&nlp_handle, &x) {
1452                    stats.final_objective = f;
1453                    stats.final_scaled_objective = f;
1454                }
1455            }
1456            // Final residuals straight off the cq cache. These mirror
1457            // the values upstream prints in its end-of-run summary
1458            // ("Dual infeasibility / Constraint violation /
1459            // Complementarity / Overall NLP error").
1460            let cq = alg.cq.borrow();
1461            stats.final_dual_inf = cq.curr_dual_infeasibility_max();
1462            stats.final_constr_viol = cq.curr_primal_infeasibility_max();
1463            // Infinity-norm complementarity, max over all four bound
1464            // blocks (s_xl·z_l, s_xu·z_u, s_sl·v_l, s_su·v_u). The
1465            // empty-bound blocks return `0` from amax(), so the max is
1466            // safe even when only one side has bounds.
1467            let compl = cq
1468                .curr_compl_x_l()
1469                .amax()
1470                .max(cq.curr_compl_x_u().amax())
1471                .max(cq.curr_compl_s_l().amax())
1472                .max(cq.curr_compl_s_u().amax());
1473            stats.final_compl = compl;
1474            stats.final_kkt_error = cq.curr_nlp_error();
1475        }
1476
1477        // Map SolverReturn → ApplicationReturnStatus per
1478        // MAIN_LOOP.md's exception table.
1479        let app_status = solver_return_to_app_status(solver_status);
1480
1481        // On convergence, fire the user-supplied callback (post-optimal
1482        // sensitivity hook, pounce#16) before flowing back through
1483        // `finalize_via_orig_nlp`. Borrowed handles into the converged
1484        // KKT state stay alive for the duration of the closure.
1485        if matches!(
1486            app_status,
1487            ApplicationReturnStatus::SolveSucceeded
1488                | ApplicationReturnStatus::SolvedToAcceptableLevel
1489        ) {
1490            if let Some(cb) = self.on_converged.as_mut() {
1491                if let Some(sd) = alg.search_dir.as_mut() {
1492                    let pd = sd.pd_solver_rc();
1493                    cb(&alg.data, &alg.cq, &nlp_handle, pd);
1494                }
1495            }
1496        }
1497
1498        // Finalize: forward the final iterate to the user's TNLP. The
1499        // returned objective is evaluated on the *user* TNLP at the
1500        // unscaled iterate, so it overrides the scaled best-effort
1501        // value stashed in `final_objective` above (the algorithm-side
1502        // `eval_f` returns `f * obj_scale_factor`).
1503        match finalize_via_orig_nlp(&nlp_handle, &alg, solver_status, app_status, &tnlp) {
1504            Ok(f_unscaled) => {
1505                self.statistics.borrow_mut().final_objective = f_unscaled;
1506            }
1507            Err(()) => {
1508                // Couldn't finalize; keep the scaled fallback and
1509                // surface the original status.
1510            }
1511        }
1512
1513        // End-of-solve timing report. Gated on `print_timing_statistics`
1514        // (default "no"); mirrors upstream's
1515        // `IpoptApplication::call_optimize` →
1516        // `IpTimingStatistics::PrintAllValues` call site. The report
1517        // goes to stdout (for parity with the banner / iter-row output
1518        // path) and is also fanned out to the journalist so an
1519        // `output_file` attached via `Initialize` picks it up.
1520        let print_timing = self
1521            .options
1522            .get_bool_value("print_timing_statistics", "")
1523            .ok()
1524            .and_then(|(v, f)| f.then_some(v))
1525            .unwrap_or(false);
1526        if print_timing {
1527            let report = timing.report();
1528            print!("{}", report);
1529            use pounce_common::journalist::{JournalCategory, JournalLevel};
1530            self.journalist.print(
1531                JournalLevel::J_SUMMARY,
1532                JournalCategory::J_TIMING_STATISTICS,
1533                &report,
1534            );
1535        }
1536
1537        app_status
1538    }
1539
1540    /// Build an [`AlgorithmBuilder`] populated from the app's
1541    /// [`OptionsList`]. Public so callers wiring the restoration
1542    /// factory can hand the *inner* IPM a builder that mirrors the
1543    /// outer's `mu_strategy`/`mu_oracle`/line-search choices —
1544    /// matching upstream `IpAlgBuilder::BuildRestoIpoptAlgorithm`,
1545    /// which reads the same `mu_strategy` option with prefix `"resto."
1546    /// + prefix` and falls back to the outer setting.
1547    pub fn algorithm_builder_from_options(&self) -> AlgorithmBuilder {
1548        let mut builder = AlgorithmBuilder::new();
1549
1550        // `mehrotra_algorithm` is parsed first so its cascading
1551        // defaults (mu_strategy=adaptive, mu_oracle=probing) can be
1552        // overridden by an explicit user setting of those keys
1553        // below. Mirrors `IpAlgBuilder.cpp:Mehrotra`.
1554        let mut mehrotra_on = false;
1555        if let Ok((v, found)) = self.options.get_string_value("mehrotra_algorithm", "") {
1556            if found && v == "yes" {
1557                mehrotra_on = true;
1558                builder.mehrotra_algorithm = true;
1559                builder.mu_strategy = MuStrategyChoice::Adaptive;
1560                builder.mu_oracle = crate::mu::adaptive::MuOracleKind::Probing;
1561                // `accept_every_trial_step` short-circuits the alpha
1562                // loop / filter — Mehrotra steps would otherwise be
1563                // rejected by the filter on LP-shaped problems because
1564                // the barrier objective is non-monotone along the
1565                // corrector. Mirrors upstream `IpAlgBuilder.cpp:Mehrotra`.
1566                builder.line_search.accept_every_trial_step = true;
1567                // Aggressive iterate-push defaults (`SetNumericValueIfUnset`
1568                // in upstream). The explicit user parses below will
1569                // overwrite these if the user set them explicitly.
1570                builder.init.bound_push = 10.0;
1571                builder.init.bound_frac = 0.2;
1572                builder.init.slack_bound_push = 10.0;
1573                builder.init.slack_bound_frac = 0.2;
1574                builder.init.bound_mult_init_val = 10.0;
1575                builder.init.constr_mult_init_max = 0.0;
1576                // `alpha_for_y=bound_mult` — Mehrotra wants the
1577                // equality multipliers to advance with the dual
1578                // alpha so they stay in step with z/v. Mirrors
1579                // upstream `IpIpoptAlg.cpp:InitializeImpl`.
1580                builder.line_search.alpha_for_y =
1581                    crate::line_search::backtracking::AlphaForY::BoundMult;
1582                // `adaptive_mu_globalization=never-monotone-mode` —
1583                // upstream `IpIpoptAlg.cpp:148-154` enforces this:
1584                // Mehrotra disables the globalization switch entirely
1585                // (no fallback to monotone mode when convergence
1586                // stalls). Required for the unsafeguarded Mehrotra
1587                // path to function.
1588                builder.mu.adaptive_mu_globalization =
1589                    crate::mu::adaptive::AdaptiveMuGlobalization::NeverMonotoneMode;
1590                // `least_square_init_primal=yes` — upstream
1591                // `IpIpoptAlg.cpp:182` enables this for the Mehrotra
1592                // cascade. Replaces the user's starting `x` with the
1593                // min-norm primal that satisfies the linearized
1594                // equality+inequality constraints. Critical on
1595                // LP-shaped problems where the user's starting point
1596                // can be wildly infeasible (e.g. nuffield2_trap).
1597                builder.init.least_square_init_primal = true;
1598            }
1599        }
1600
1601        if let Ok((v, found)) = self.options.get_string_value("mu_strategy", "") {
1602            if found {
1603                let parsed = match v.as_str() {
1604                    "adaptive" => MuStrategyChoice::Adaptive,
1605                    _ => MuStrategyChoice::Monotone,
1606                };
1607                if mehrotra_on && matches!(parsed, MuStrategyChoice::Monotone) {
1608                    // Upstream Ipopt refuses this combination: Mehrotra
1609                    // needs an affine step every iter, which only the
1610                    // adaptive path computes. Keep adaptive and warn.
1611                    tracing::warn!(target: "pounce::algorithm",
1612                        "pounce: mehrotra_algorithm=yes requires \
1613                         mu_strategy=adaptive; ignoring \
1614                         mu_strategy=monotone."
1615                    );
1616                } else {
1617                    builder.mu_strategy = parsed;
1618                }
1619            }
1620        }
1621        if let Ok((v, found)) = self.options.get_string_value("mu_oracle", "") {
1622            if found {
1623                builder.mu_oracle = match v.as_str() {
1624                    "loqo" => crate::mu::adaptive::MuOracleKind::Loqo,
1625                    "probing" => crate::mu::adaptive::MuOracleKind::Probing,
1626                    _ => crate::mu::adaptive::MuOracleKind::QualityFunction,
1627                };
1628            }
1629        }
1630        if let Ok((v, found)) = self
1631            .options
1632            .get_string_value("adaptive_mu_globalization", "")
1633        {
1634            if found {
1635                use crate::mu::adaptive::AdaptiveMuGlobalization;
1636                builder.mu.adaptive_mu_globalization = match v.as_str() {
1637                    "kkt-error" => AdaptiveMuGlobalization::KktError,
1638                    "never-monotone-mode" => AdaptiveMuGlobalization::NeverMonotoneMode,
1639                    _ => AdaptiveMuGlobalization::ObjConstrFilter,
1640                };
1641            }
1642        }
1643        if let Ok((v, found)) = self.options.get_string_value("hessian_approximation", "") {
1644            if found {
1645                builder.hessian_approximation = match v.as_str() {
1646                    "limited-memory" => HessianApproxChoice::LimitedMemory,
1647                    _ => HessianApproxChoice::Exact,
1648                };
1649            }
1650        }
1651        // Limited-memory quasi-Newton update formula. Registered upstream
1652        // (`limited_memory_update_type`, IpLimMemQuasiNewtonUpdater.cpp) but
1653        // until now read nowhere on the IPM path — the updater was hard-wired
1654        // to Powell-damped BFGS. SR1 is honored too (the updater and the
1655        // low-rank/inertia path already handle its indefinite models).
1656        if let Ok((v, found)) = self
1657            .options
1658            .get_string_value("limited_memory_update_type", "")
1659        {
1660            if found {
1661                builder.limited_memory_update_type = match v.as_str() {
1662                    "sr1" => UpdateType::Sr1,
1663                    _ => UpdateType::Bfgs,
1664                };
1665            }
1666        }
1667        // Limited-memory history length (`limited_memory_max_history`).
1668        if let Ok((v, found)) = self
1669            .options
1670            .get_integer_value("limited_memory_max_history", "")
1671        {
1672            if found && v >= 0 {
1673                builder.limited_memory_max_history = v as Index;
1674            }
1675        }
1676        if let Ok((v, found)) = self.options.get_string_value("line_search_method", "") {
1677            if found {
1678                builder.line_search_method = match v.as_str() {
1679                    "cg-penalty" => LineSearchChoice::CgPenalty,
1680                    "penalty" => LineSearchChoice::Penalty,
1681                    _ => LineSearchChoice::Filter,
1682                };
1683            }
1684        }
1685        // `accept_every_trial_step` — direct user override. Parsed
1686        // after the Mehrotra cascade so an explicit `no` still wins.
1687        if let Ok((v, found)) = self.options.get_string_value("accept_every_trial_step", "") {
1688            if found {
1689                builder.line_search.accept_every_trial_step = v == "yes";
1690            }
1691        }
1692        // `alpha_for_y` — direct user override. Parsed after the
1693        // Mehrotra cascade so an explicit value still wins.
1694        if let Ok((v, found)) = self.options.get_string_value("alpha_for_y", "") {
1695            if found {
1696                use crate::line_search::backtracking::AlphaForY;
1697                builder.line_search.alpha_for_y = match v.as_str() {
1698                    "primal" => AlphaForY::Primal,
1699                    "bound-mult" | "bound_mult" => AlphaForY::BoundMult,
1700                    "full" => AlphaForY::Full,
1701                    "min" => AlphaForY::Min,
1702                    "max" => AlphaForY::Max,
1703                    "primal-and-full" | "dual-and-full" => AlphaForY::Primal,
1704                    _ => AlphaForY::Primal,
1705                };
1706            }
1707        }
1708        // `nlp_scaling_method` is consumed NLP-side in
1709        // `OrigIpoptNlp::determine_scaling_from_starting_point` (see the
1710        // `determine_scaling_from_starting_point` call earlier in this
1711        // method); there is no algorithm-side scaling strategy to wire.
1712
1713        // Unlike the other options here, we always honor the registry
1714        // value (not just when the user set it explicitly): the option
1715        // registry default is "ma57" but `AlgorithmBuilder::default`
1716        // has `linear_solver: Feral`, so gating on `found` would
1717        // silently route default runs through Feral while the banner
1718        // (and ipopt-compatible behavior) advertises MA57.
1719        if let Ok((v, _found)) = self.options.get_string_value("linear_solver", "") {
1720            builder.linear_solver = match v.as_str() {
1721                "ma57" => LinearSolverChoice::Ma57,
1722                _ => LinearSolverChoice::Feral,
1723            };
1724        }
1725
1726        // `linear_system_scaling` — symmetric scaling of the augmented
1727        // KKT matrix before factorization. Port of
1728        // `IpTSymLinearSolver.cpp:RegisterOptions` plumbing. Default
1729        // "none"; "ruiz" invokes the Ruiz-2001 symmetric ∞-norm
1730        // equilibration in `RuizTSymScalingMethod`. "mc19" and
1731        // "slack-based" are accepted by the registry but not yet
1732        // implemented at this layer; they fall back to no scaling
1733        // with a one-line stderr notice.
1734        if let Ok((v, found)) = self.options.get_string_value("linear_system_scaling", "") {
1735            if found {
1736                builder.linear_system_scaling = match v.as_str() {
1737                    "ruiz" => crate::alg_builder::LinearSystemScalingChoice::Ruiz,
1738                    "mc19" => crate::alg_builder::LinearSystemScalingChoice::Mc19,
1739                    _ => crate::alg_builder::LinearSystemScalingChoice::None,
1740                };
1741            }
1742        }
1743        if let Ok((v, found)) = self.options.get_bool_value("linear_scaling_on_demand", "") {
1744            if found {
1745                builder.linear_scaling_on_demand = v;
1746            }
1747        }
1748
1749        // Convergence tolerances (port of `IpOptErrorConvCheck.cpp`'s
1750        // `RegisterOptions` consumers). Defaults already match upstream
1751        // — only override when the user set the key explicitly.
1752        let read_num = |key: &str| -> Option<f64> {
1753            self.options
1754                .get_numeric_value(key, "")
1755                .ok()
1756                .and_then(|(v, f)| f.then_some(v))
1757        };
1758        let read_int = |key: &str| -> Option<i32> {
1759            self.options
1760                .get_integer_value(key, "")
1761                .ok()
1762                .and_then(|(v, f)| f.then_some(v))
1763        };
1764        if let Some(v) = read_num("tol") {
1765            builder.conv_check.tol = v;
1766        }
1767        if let Some(v) = read_num("dual_inf_tol") {
1768            builder.conv_check.dual_inf_tol = v;
1769        }
1770        if let Some(v) = read_num("constr_viol_tol") {
1771            builder.conv_check.constr_viol_tol = v;
1772        }
1773        if let Some(v) = read_num("compl_inf_tol") {
1774            builder.conv_check.compl_inf_tol = v;
1775        }
1776        if let Some(v) = read_int("max_iter") {
1777            builder.conv_check.max_iter = v;
1778        }
1779        if let Some(v) = read_num("max_cpu_time") {
1780            builder.conv_check.max_cpu_time = v;
1781        }
1782        if let Some(v) = read_num("max_wall_time") {
1783            builder.conv_check.max_wall_time = v;
1784        }
1785        if let Some(v) = read_num("acceptable_tol") {
1786            builder.conv_check.acceptable_tol = v;
1787        }
1788        if let Some(v) = read_num("acceptable_dual_inf_tol") {
1789            builder.conv_check.acceptable_dual_inf_tol = v;
1790        }
1791        if let Some(v) = read_num("acceptable_constr_viol_tol") {
1792            builder.conv_check.acceptable_constr_viol_tol = v;
1793        }
1794        if let Some(v) = read_num("acceptable_compl_inf_tol") {
1795            builder.conv_check.acceptable_compl_inf_tol = v;
1796        }
1797        if let Some(v) = read_num("acceptable_obj_change_tol") {
1798            builder.conv_check.acceptable_obj_change_tol = v;
1799        }
1800        if let Some(v) = read_int("acceptable_iter") {
1801            builder.conv_check.acceptable_iter = v;
1802        }
1803        if let Some(v) = read_num("infeas_stationarity_tol") {
1804            builder.conv_check.infeas_stationarity_tol = v;
1805        }
1806        if let Some(v) = read_num("infeas_viol_kappa") {
1807            builder.conv_check.infeas_viol_kappa = v;
1808        }
1809        if let Some(v) = read_int("infeas_max_streak") {
1810            builder.conv_check.infeas_max_streak = v;
1811        }
1812
1813        // Barrier-parameter (μ) options — consumers in
1814        // `IpMonotoneMuUpdate.cpp` / `IpAdaptiveMuUpdate.cpp`. Both
1815        // updaters share the same option names; the builder forwards
1816        // each into whichever strategy is assembled.
1817        if let Some(v) = read_num("mu_init") {
1818            builder.mu.mu_init = v;
1819        }
1820        if let Some(v) = read_num("mu_max") {
1821            builder.mu.mu_max = v;
1822        }
1823        if let Some(v) = read_num("mu_max_fact") {
1824            builder.mu.mu_max_fact = v;
1825        }
1826        if let Some(v) = read_num("mu_min") {
1827            builder.mu.mu_min = v;
1828        }
1829        if let Some(v) = read_num("mu_target") {
1830            builder.mu.mu_target = v;
1831        }
1832        if let Some(v) = read_num("mu_linear_decrease_factor") {
1833            builder.mu.mu_linear_decrease_factor = v;
1834        }
1835        if let Some(v) = read_num("mu_superlinear_decrease_power") {
1836            builder.mu.mu_superlinear_decrease_power = v;
1837        }
1838        if let Ok((v, found)) = self
1839            .options
1840            .get_string_value("mu_allow_fast_monotone_decrease", "")
1841        {
1842            if found {
1843                builder.mu.mu_allow_fast_monotone_decrease = v == "yes";
1844            }
1845        }
1846        if let Some(v) = read_num("barrier_tol_factor") {
1847            builder.mu.barrier_tol_factor = v;
1848        }
1849        if let Some(v) = read_num("sigma_max") {
1850            builder.mu.sigma_max = v;
1851        }
1852        if let Some(v) = read_num("sigma_min") {
1853            builder.mu.sigma_min = v;
1854        }
1855
1856        // Quality-function oracle knobs — consumers in
1857        // `IpQualityFunctionMuOracle.cpp:RegisterOptions`. Forwarded
1858        // to the oracle on every free-mode call.
1859        if let Ok((v, found)) = self
1860            .options
1861            .get_string_value("quality_function_norm_type", "")
1862        {
1863            if found {
1864                use crate::mu::oracle::quality_function::NormType;
1865                builder.mu.quality_function_norm_type = match v.as_str() {
1866                    "1-norm" => NormType::OneNorm,
1867                    "2-norm" => NormType::TwoNorm,
1868                    "max-norm" => NormType::MaxNorm,
1869                    _ => NormType::TwoNormSquared,
1870                };
1871            }
1872        }
1873        if let Ok((v, found)) = self
1874            .options
1875            .get_string_value("quality_function_centrality", "")
1876        {
1877            if found {
1878                use crate::mu::oracle::quality_function::CentralityType;
1879                builder.mu.quality_function_centrality = match v.as_str() {
1880                    "log" => CentralityType::LogCenter,
1881                    "reciprocal" => CentralityType::ReciprocalCenter,
1882                    "cubed-reciprocal" => CentralityType::CubedReciprocalCenter,
1883                    _ => CentralityType::None,
1884                };
1885            }
1886        }
1887        if let Ok((v, found)) = self
1888            .options
1889            .get_string_value("quality_function_balancing_term", "")
1890        {
1891            if found {
1892                use crate::mu::oracle::quality_function::BalancingTermType;
1893                builder.mu.quality_function_balancing_term = match v.as_str() {
1894                    "cubic" => BalancingTermType::CubicTerm,
1895                    _ => BalancingTermType::None,
1896                };
1897            }
1898        }
1899        if let Some(v) = read_int("quality_function_max_section_steps") {
1900            builder.mu.quality_function_max_section_steps = v;
1901        }
1902        if let Some(v) = read_num("quality_function_section_sigma_tol") {
1903            builder.mu.quality_function_section_sigma_tol = v;
1904        }
1905        if let Some(v) = read_num("quality_function_section_qf_tol") {
1906            builder.mu.quality_function_section_qf_tol = v;
1907        }
1908
1909        // `probing_iterate_quality_factor` — pounce-specific guard
1910        // (pounce#58) on the probing μ-oracle's input iterate. When
1911        // `curr_avrg_compl / curr_mu` exceeds this factor, the
1912        // μ-update layer signals restoration via
1913        // `IpoptData::request_resto` instead of letting probing
1914        // return `σ · mu_curr` ≫ previous μ. Default 1e4; set to ≤ 0
1915        // to disable. No upstream Ipopt counterpart.
1916        if let Some(v) = read_num("probing_iterate_quality_factor") {
1917            builder.mu.probing_iterate_quality_factor = v;
1918        }
1919
1920        // Adaptive-μ extras — consumers in
1921        // `IpAdaptiveMuUpdate.cpp:RegisterOptions`. Only active when
1922        // `mu_strategy=adaptive`.
1923        if let Some(v) = read_num("adaptive_mu_safeguard_factor") {
1924            builder.mu.adaptive_mu_safeguard_factor = v;
1925        }
1926        if let Some(v) = read_num("adaptive_mu_monotone_init_factor") {
1927            builder.mu.adaptive_mu_monotone_init_factor = v;
1928        }
1929        if let Ok((v, found)) = self
1930            .options
1931            .get_bool_value("adaptive_mu_restore_previous_iterate", "")
1932        {
1933            if found {
1934                builder.mu.adaptive_mu_restore_previous_iterate = v;
1935            }
1936        }
1937        if let Some(v) = read_int("adaptive_mu_kkterror_red_iters") {
1938            if v >= 0 {
1939                builder.mu.adaptive_mu_kkterror_red_iters = v as usize;
1940            }
1941        }
1942        if let Some(v) = read_num("adaptive_mu_kkterror_red_fact") {
1943            builder.mu.adaptive_mu_kkterror_red_fact = v;
1944        }
1945        if let Ok((v, found)) = self
1946            .options
1947            .get_string_value("adaptive_mu_kkt_norm_type", "")
1948        {
1949            if found {
1950                use crate::mu::adaptive::AdaptiveMuKktNorm;
1951                builder.mu.adaptive_mu_kkt_norm_type = match v.as_str() {
1952                    "1-norm" => AdaptiveMuKktNorm::OneNorm,
1953                    "2-norm" => AdaptiveMuKktNorm::TwoNorm,
1954                    "max-norm" => AdaptiveMuKktNorm::MaxNorm,
1955                    _ => AdaptiveMuKktNorm::TwoNormSquared,
1956                };
1957            }
1958        }
1959
1960        // Watchdog options — consumers in
1961        // `IpBacktrackingLineSearch.cpp:RegisterOptions`. Baked into
1962        // the `BacktrackingLineSearch` at build time.
1963        if let Some(v) = read_int("watchdog_shortened_iter_trigger") {
1964            builder.line_search.watchdog_shortened_iter_trigger = v;
1965        }
1966        if let Some(v) = read_int("watchdog_trial_iter_max") {
1967            builder.line_search.watchdog_trial_iter_max = v;
1968        }
1969        if let Some(v) = read_num("soft_resto_pderror_reduction_factor") {
1970            builder.line_search.soft_resto_pderror_reduction_factor = v;
1971        }
1972        if let Some(v) = read_int("max_soft_resto_iters") {
1973            builder.line_search.max_soft_resto_iters = v;
1974        }
1975
1976        // Iteration-output options — consumed by `OrigIterationOutput`.
1977        if let Some(v) = read_int("print_frequency_iter") {
1978            builder.output.print_frequency_iter = v;
1979        }
1980        if let Some(v) = read_num("print_frequency_time") {
1981            builder.output.print_frequency_time = v;
1982        }
1983        if let Ok((v, found)) = self.options.get_bool_value("print_info_string", "") {
1984            if found {
1985                builder.output.print_info_string = v;
1986            }
1987        }
1988        if let Ok((v, found)) = self.options.get_string_value("inf_pr_output", "") {
1989            if found {
1990                builder.output.inf_pr_output_internal = v == "internal";
1991            }
1992        }
1993
1994        // Warm-start options — consumed by `WarmStartIterateInitializer`
1995        // (port of `IpWarmStartIterateInitializer.cpp:RegisterOptions`).
1996        // `warm_start_init_point` is the toggle that picks between the
1997        // default (cold) and warm-start initializers; the remaining
1998        // knobs are baked onto the chosen initializer at build time.
1999        if let Ok((v, found)) = self.options.get_bool_value("warm_start_init_point", "") {
2000            if found {
2001                builder.warm_start_init_point = v;
2002            }
2003        }
2004        if let Ok((v, found)) = self.options.get_bool_value("warm_start_same_structure", "") {
2005            if found {
2006                builder.warm.same_structure = v;
2007            }
2008        }
2009        if let Some(v) = read_num("warm_start_bound_push") {
2010            builder.warm.bound_push = v;
2011        }
2012        if let Some(v) = read_num("warm_start_bound_frac") {
2013            builder.warm.bound_frac = v;
2014        }
2015        if let Some(v) = read_num("warm_start_slack_bound_push") {
2016            builder.warm.slack_bound_push = v;
2017        }
2018        if let Some(v) = read_num("warm_start_slack_bound_frac") {
2019            builder.warm.slack_bound_frac = v;
2020        }
2021        if let Some(v) = read_num("warm_start_mult_bound_push") {
2022            builder.warm.mult_bound_push = v;
2023        }
2024        if let Some(v) = read_num("warm_start_mult_init_max") {
2025            builder.warm.mult_init_max = v;
2026        }
2027        if let Some(v) = read_num("warm_start_target_mu") {
2028            builder.warm.target_mu = v;
2029        }
2030        if let Ok((v, found)) = self
2031            .options
2032            .get_string_value("warm_start_entire_iterate", "")
2033        {
2034            if found {
2035                builder.warm.entire_iterate = v == "yes";
2036            }
2037        }
2038
2039        // `DefaultIterateInitializer` knobs — parsed after the Mehrotra
2040        // cascade so explicit user values win
2041        // (mirrors upstream's `SetNumericValueIfUnset` semantics).
2042        if let Some(v) = read_num("bound_push") {
2043            builder.init.bound_push = v;
2044        }
2045        if let Some(v) = read_num("bound_frac") {
2046            builder.init.bound_frac = v;
2047        }
2048        if let Some(v) = read_num("slack_bound_push") {
2049            builder.init.slack_bound_push = v;
2050        }
2051        if let Some(v) = read_num("slack_bound_frac") {
2052            builder.init.slack_bound_frac = v;
2053        }
2054        if let Some(v) = read_num("constr_mult_init_max") {
2055            builder.init.constr_mult_init_max = v;
2056        }
2057        if let Some(v) = read_num("bound_mult_init_val") {
2058            builder.init.bound_mult_init_val = v;
2059        }
2060        if let Ok((v, found)) = self.options.get_string_value("bound_mult_init_method", "") {
2061            if found {
2062                builder.init.bound_mult_init_method = v;
2063            }
2064        }
2065        if let Ok((v, found)) = self
2066            .options
2067            .get_string_value("least_square_init_primal", "")
2068        {
2069            if found {
2070                builder.init.least_square_init_primal = v == "yes";
2071            }
2072        }
2073        builder
2074    }
2075}
2076
2077/// Map the integer `print_level` / `file_print_level` option to the
2078/// matching [`JournalLevel`] variant. Mirrors upstream's
2079/// `static_cast<EJournalLevel>(int_value)` with clamping.
2080/// The eight block dimensions of an iterate, in canonical order
2081/// (x, s, y_c, y_d, z_l, z_u, v_l, v_u). Used to guard the debugger's
2082/// warm-restart install against a structural mismatch between solves.
2083fn iterates_dims(c: &IteratesVector) -> [i32; 8] {
2084    [
2085        c.x.dim(),
2086        c.s.dim(),
2087        c.y_c.dim(),
2088        c.y_d.dim(),
2089        c.z_l.dim(),
2090        c.z_u.dim(),
2091        c.v_l.dim(),
2092        c.v_u.dim(),
2093    ]
2094}
2095
2096fn journal_level_from_int(v: i32) -> JournalLevel {
2097    match v.clamp(0, 12) {
2098        0 => JournalLevel::J_NONE,
2099        1 => JournalLevel::J_ERROR,
2100        2 => JournalLevel::J_STRONGWARNING,
2101        3 => JournalLevel::J_SUMMARY,
2102        4 => JournalLevel::J_WARNING,
2103        5 => JournalLevel::J_ITERSUMMARY,
2104        6 => JournalLevel::J_DETAILED,
2105        7 => JournalLevel::J_MOREDETAILED,
2106        8 => JournalLevel::J_VECTOR,
2107        9 => JournalLevel::J_MOREVECTOR,
2108        10 => JournalLevel::J_MATRIX,
2109        11 => JournalLevel::J_MOREMATRIX,
2110        _ => JournalLevel::J_ALL,
2111    }
2112}
2113
2114/// Default symmetric linear-solver factory, parameterized by the
2115/// pounce-extension FERAL knobs read off the application's
2116/// `OptionsList`.
2117///
2118/// FERAL (pure-Rust) is the shipping default. The HSL MA57 backend is
2119/// available when the `ma57` cargo feature is enabled; without it,
2120/// requesting `linear_solver = ma57` falls back to FERAL with a
2121/// warning printed by the journalist (see [`AlgorithmBuilder`]).
2122pub fn default_backend_factory(feral_cfg: pounce_feral::FeralConfig) -> LinearBackendFactory {
2123    Box::new(
2124        move |choice: LinearSolverChoice| -> Box<dyn SparseSymLinearSolverInterface> {
2125            match choice {
2126                LinearSolverChoice::Feral => Box::new(
2127                    pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone()),
2128                ),
2129                LinearSolverChoice::Ma57 => {
2130                    #[cfg(feature = "ma57")]
2131                    {
2132                        Box::new(pounce_hsl::Ma57SolverInterface::new())
2133                    }
2134                    #[cfg(not(feature = "ma57"))]
2135                    {
2136                        // ma57 feature not compiled in — fall back to FERAL.
2137                        Box::new(pounce_feral::FeralSolverInterface::with_config(
2138                            feral_cfg.clone(),
2139                        ))
2140                    }
2141                }
2142            }
2143        },
2144    )
2145}
2146
2147/// Sink-aware variant of [`default_backend_factory`]. Identical
2148/// dispatch, but the FERAL backend is constructed with a
2149/// `LinearSolverSummary` sink so [`IpoptApplication`] can read out
2150/// aggregate post-mortem stats (factor counts, fill ratio, extremal
2151/// pivots, final inertia) after the solve returns. MA57 ignores the
2152/// sink — the HSL backend doesn't carry the same instrumentation yet.
2153pub fn default_backend_factory_with_sink(
2154    feral_cfg: pounce_feral::FeralConfig,
2155    sink: Arc<Mutex<LinearSolverSummary>>,
2156) -> LinearBackendFactory {
2157    Box::new(
2158        move |choice: LinearSolverChoice| -> Box<dyn SparseSymLinearSolverInterface> {
2159            match choice {
2160                LinearSolverChoice::Feral => Box::new(
2161                    pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone())
2162                        .with_summary_sink(Arc::clone(&sink)),
2163                ),
2164                LinearSolverChoice::Ma57 => {
2165                    #[cfg(feature = "ma57")]
2166                    {
2167                        Box::new(pounce_hsl::Ma57SolverInterface::new())
2168                    }
2169                    #[cfg(not(feature = "ma57"))]
2170                    {
2171                        Box::new(
2172                            pounce_feral::FeralSolverInterface::with_config(feral_cfg.clone())
2173                                .with_summary_sink(Arc::clone(&sink)),
2174                        )
2175                    }
2176                }
2177            }
2178        },
2179    )
2180}
2181
2182/// Read the `feral_*` extension options off `options`, falling
2183/// back to the env-var defaults baked into [`pounce_feral::FeralConfig::from_env`]
2184/// for any knob the caller did not set explicitly. The returned
2185/// config is what every default-factory invocation (main IPM and
2186/// restoration sub-IPM) consumes.
2187pub fn feral_config_from_options(
2188    options: &pounce_common::options_list::OptionsList,
2189) -> pounce_feral::FeralConfig {
2190    let mut cfg = pounce_feral::FeralConfig::from_env();
2191    // Tri-state: the `(_, true)` arm only fires when the user set the
2192    // option explicitly. Leaving it unset keeps `cfg.cascade_break` at
2193    // `None`, which inherits FERAL's `NumericParams::default()` (CB on
2194    // as of FERAL Phase B / pounce#55). `Some(false)` explicitly
2195    // disarms (reproduces pre-Phase-B behaviour, surfaces FERAL's
2196    // `DelayBudgetExceeded` on non-root cascade victims).
2197    if let Ok((v, true)) = options.get_bool_value("feral_cascade_break", "") {
2198        cfg.cascade_break = Some(v);
2199    }
2200    if let Ok((v, true)) = options.get_bool_value("feral_fma", "") {
2201        cfg.fma = v;
2202    }
2203    if let Ok((v, true)) = options.get_bool_value("feral_refine", "") {
2204        cfg.refine = v;
2205    }
2206    if let Ok((v, true)) = options.get_numeric_value("feral_singular_pivot_floor", "") {
2207        cfg.singular_pivot_floor = v;
2208    }
2209    if let Ok((v, true)) = options.get_numeric_value("feral_pivtol", "") {
2210        cfg.pivtol = v;
2211    }
2212    // Only override on explicit set so `from_env` (which itself
2213    // defaults to OrderingMethod::Auto) keeps governing unset cases.
2214    // Unrecognized tags are silently ignored — the registered enum
2215    // restricts inputs at the OptionsList layer.
2216    if let Ok((v, true)) = options.get_string_value("feral_ordering", "") {
2217        if let Some(m) = pounce_feral::parse_ordering_method(&v) {
2218            cfg.ordering = m;
2219        }
2220    }
2221    // Same explicit-set discipline as `feral_ordering`: `from_env`
2222    // defaults to ScalingStrategy::Auto (FERAL's current default), so
2223    // leaving the option unset preserves existing behaviour exactly.
2224    if let Ok((v, true)) = options.get_string_value("feral_scaling", "") {
2225        if let Some(s) = pounce_feral::parse_scaling_strategy(&v) {
2226            cfg.scaling = s;
2227        }
2228    }
2229    cfg
2230}
2231
2232/// Map upstream `SolverReturn` codes to `ApplicationReturnStatus`.
2233/// Mirrors the table in
2234/// `ref/Ipopt/AGENT_REFERENCE/MAIN_LOOP.md` ("exception → SolverReturn
2235/// map") and the corresponding switch in
2236/// `IpIpoptApplication.cpp:call_optimize`.
2237fn solver_return_to_app_status(s: SolverReturn) -> ApplicationReturnStatus {
2238    match s {
2239        SolverReturn::Success => ApplicationReturnStatus::SolveSucceeded,
2240        SolverReturn::StopAtAcceptablePoint => ApplicationReturnStatus::SolvedToAcceptableLevel,
2241        SolverReturn::FeasiblePointFound => ApplicationReturnStatus::FeasiblePointFound,
2242        SolverReturn::MaxiterExceeded => ApplicationReturnStatus::MaximumIterationsExceeded,
2243        SolverReturn::CpuTimeExceeded => ApplicationReturnStatus::MaximumCpuTimeExceeded,
2244        SolverReturn::WallTimeExceeded => ApplicationReturnStatus::MaximumWallTimeExceeded,
2245        SolverReturn::StopAtTinyStep => ApplicationReturnStatus::SearchDirectionBecomesTooSmall,
2246        SolverReturn::LocalInfeasibility => ApplicationReturnStatus::InfeasibleProblemDetected,
2247        SolverReturn::UserRequestedStop => ApplicationReturnStatus::UserRequestedStop,
2248        SolverReturn::DivergingIterates => ApplicationReturnStatus::DivergingIterates,
2249        SolverReturn::RestorationFailure => ApplicationReturnStatus::RestorationFailed,
2250        SolverReturn::ErrorInStepComputation => ApplicationReturnStatus::ErrorInStepComputation,
2251        SolverReturn::InvalidNumberDetected => ApplicationReturnStatus::InvalidNumberDetected,
2252        SolverReturn::TooFewDegreesOfFreedom => ApplicationReturnStatus::NotEnoughDegreesOfFreedom,
2253        SolverReturn::InvalidOption => ApplicationReturnStatus::InvalidOption,
2254        SolverReturn::OutOfMemory => ApplicationReturnStatus::InsufficientMemory,
2255        SolverReturn::InternalError | SolverReturn::Unassigned => {
2256            ApplicationReturnStatus::InternalError
2257        }
2258    }
2259}
2260
2261/// Best-effort evaluation of the objective at the algorithm's final
2262/// `x`. Returns the *scaled* objective (`f * obj_scale_factor`); used
2263/// to populate `SolveStatistics::final_scaled_objective`.
2264fn try_eval_curr_f(
2265    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2266    x: &Rc<dyn pounce_linalg::Vector>,
2267) -> Result<Number, ()> {
2268    let mut nlp_mut = nlp.borrow_mut();
2269    Ok(nlp_mut.eval_f(&**x))
2270}
2271
2272/// Trigger predicate for the Phase-3.5 ℓ₁ auto-fallback path. Returns
2273/// `true` when a status warrants a retry through the wrapper. Mirrors
2274/// ripopt#23's trigger set, extended per the audit's Refinement B
2275/// (pounce-side `Not_Enough_Degrees_Of_Freedom` is added because
2276/// pounce's DOF early-exit blocks NE-suffix problems that ripopt's
2277/// equivalent would let pass to the wrapper).
2278fn is_l1_fallback_trigger(status: ApplicationReturnStatus) -> bool {
2279    matches!(
2280        status,
2281        ApplicationReturnStatus::RestorationFailed
2282            | ApplicationReturnStatus::InfeasibleProblemDetected
2283            | ApplicationReturnStatus::SolvedToAcceptableLevel
2284            | ApplicationReturnStatus::MaximumIterationsExceeded
2285            | ApplicationReturnStatus::NotEnoughDegreesOfFreedom
2286    )
2287}
2288
2289/// Forward the final iterate back to the user's `TNLP::finalize_solution`.
2290/// We pull `x` (compressed in `x_var`-space) off the algorithm's
2291/// `data.curr`, lift it back to full TNLP indexing, and pass empty
2292/// multipliers for now (the algorithm's `y_c`, `y_d`, `z_l`, `z_u` are
2293/// in compressed split form — re-assembling them into the user's
2294/// `lambda` / `z_l` / `z_u` is mechanical but lives behind a
2295/// `OrigIpoptNlp::finalize_solution_*` accessor that's still being
2296/// fleshed out). On success returns the unscaled objective evaluated
2297/// on the user TNLP at the final iterate; returns `Err` if the final
2298/// iterate is missing.
2299fn finalize_via_orig_nlp(
2300    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2301    alg: &IpoptAlgorithm,
2302    solver_status: SolverReturn,
2303    _app_status: ApplicationReturnStatus,
2304    tnlp: &Rc<RefCell<dyn TNLP>>,
2305) -> Result<Number, ()> {
2306    let curr = alg.data.borrow().curr.clone().ok_or(())?;
2307    // Lift compressed x_var → full-x (length `info.n`) so the user
2308    // TNLP receives the same shape it provided. With `make_parameter`
2309    // the fixed components are spliced back in by the IpoptNlp.
2310    let nlp_borrow = nlp.borrow();
2311    let x_vec: Vec<Number> = nlp_borrow.lift_x_to_full(&*curr.x);
2312    let info = tnlp.borrow_mut().get_nlp_info().ok_or(())?;
2313    let n = info.n as usize;
2314    let m = info.m as usize;
2315    debug_assert_eq!(x_vec.len(), n);
2316    // Lift algorithm-side multipliers back into user-space (pounce#11).
2317    // Use the `finalize_solution_*` family (not the `pack_*` family): the
2318    // final solution duals must be reported in the user's *unscaled-
2319    // Lagrangian* convention `∇f + λ·∇g + z = 0`, which divides out the
2320    // `obj_scale_factor` the algorithm threads through `eval_h`. The `pack_*`
2321    // family deliberately omits that division because it feeds the scaled
2322    // `eval_h`; calling it here left every dual scaled by `obj_scale_factor`
2323    // whenever gradient-based scaling triggered (pounce#11 F1).
2324    // Backends without overrides return empty; fall back to zero stubs so the
2325    // user sees a length-consistent vector.
2326    let mut z_l = nlp_borrow.finalize_solution_z_l(&*curr.z_l);
2327    if z_l.is_empty() {
2328        z_l = vec![0.0; n];
2329    }
2330    let mut z_u = nlp_borrow.finalize_solution_z_u(&*curr.z_u);
2331    if z_u.is_empty() {
2332        z_u = vec![0.0; n];
2333    }
2334    let mut lambda = nlp_borrow.finalize_solution_lambda(&*curr.y_c, &*curr.y_d);
2335    if lambda.is_empty() {
2336        lambda = vec![0.0; m];
2337    }
2338    drop(nlp_borrow);
2339    // Compute g(x) via the user TNLP so the final residual is
2340    // populated for the user.
2341    let mut g_final = vec![0.0; m];
2342    let _ = tnlp.borrow_mut().eval_g(&x_vec, true, &mut g_final);
2343    let f_final = tnlp
2344        .borrow_mut()
2345        .eval_f(&x_vec, true)
2346        .unwrap_or(Number::NAN);
2347    tnlp.borrow_mut().finalize_solution(
2348        Solution {
2349            status: solver_status,
2350            x: &x_vec,
2351            z_l: &z_l,
2352            z_u: &z_u,
2353            g: &g_final,
2354            lambda: &lambda,
2355            obj_value: f_final,
2356        },
2357        &TnlpIpoptData::default(),
2358        &TnlpIpoptCq::default(),
2359    );
2360    Ok(f_final)
2361}
2362
2363/// Bind SQP suboptions registered in `upstream_options.rs`
2364/// (`sqp_globalization`, `sqp_hessian`, `sqp_max_iter`, `sqp_tol`,
2365/// `sqp_constr_viol_tol`, `sqp_dual_inf_tol`, `sqp_l1_penalty`,
2366/// `sqp_bt_reduction`, `sqp_bt_min_alpha`, `sqp_print_level`,
2367/// `sqp_lbfgs_max_history`) onto
2368/// `opts`. Used by [`IpoptApplication::algorithm_builder_snapshot`]
2369/// before constructing an SQP algorithm.
2370fn apply_sqp_options(options: &OptionsList, opts: &mut crate::sqp::SqpOptions) {
2371    use crate::sqp::{SqpGlobalization, SqpHessianSource};
2372
2373    if let Ok((s, true)) = options.get_string_value("sqp_globalization", "") {
2374        opts.globalization = match s.as_str() {
2375            "filter" => SqpGlobalization::Filter,
2376            "l1-elastic" => SqpGlobalization::L1Elastic,
2377            _ => opts.globalization,
2378        };
2379    }
2380    if let Ok((s, true)) = options.get_string_value("sqp_hessian", "") {
2381        opts.hessian = match s.as_str() {
2382            "exact" => SqpHessianSource::Exact,
2383            "damped-bfgs" => SqpHessianSource::DampedBfgs,
2384            "lbfgs" => SqpHessianSource::Lbfgs,
2385            _ => opts.hessian,
2386        };
2387    }
2388    if let Ok((v, true)) = options.get_integer_value("sqp_max_iter", "") {
2389        if v >= 0 {
2390            opts.max_iter = v as u32;
2391        }
2392    }
2393    if let Ok((v, true)) = options.get_numeric_value("sqp_tol", "") {
2394        opts.tol = v;
2395    }
2396    if let Ok((v, true)) = options.get_numeric_value("sqp_constr_viol_tol", "") {
2397        opts.constr_viol_tol = v;
2398    }
2399    if let Ok((v, true)) = options.get_numeric_value("sqp_dual_inf_tol", "") {
2400        opts.dual_inf_tol = v;
2401    }
2402    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty", "") {
2403        opts.l1_penalty = v;
2404    }
2405    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty_safety", "") {
2406        opts.l1_penalty_safety = v;
2407    }
2408    if let Ok((v, true)) = options.get_numeric_value("sqp_l1_penalty_max", "") {
2409        opts.l1_penalty_max = v;
2410    }
2411    if let Ok((v, true)) = options.get_numeric_value("sqp_bt_reduction", "") {
2412        opts.bt_reduction = v;
2413    }
2414    if let Ok((v, true)) = options.get_numeric_value("sqp_bt_min_alpha", "") {
2415        opts.bt_min_alpha = v;
2416    }
2417    if let Ok((v, true)) = options.get_integer_value("sqp_print_level", "") {
2418        opts.print_level = v.clamp(0, u8::MAX as i32) as u8;
2419    }
2420    if let Ok((v, true)) = options.get_integer_value("sqp_lbfgs_max_history", "") {
2421        if v >= 1 {
2422            opts.lbfgs_max_history = v as u32;
2423        }
2424    }
2425}
2426
2427/// Populate the active-set SQP **QP-subproblem** options
2428/// ([`pounce_qp::QpOptions`]) from the `sqp_qp_*` option family.
2429///
2430/// Sister to [`apply_sqp_options`], which handles the SQP *outer-loop*
2431/// options ([`crate::sqp::SqpOptions`]); this one feeds the inner QP
2432/// solver that `SqpAlgorithm` delegates each subproblem to. Consulted
2433/// only on the `ActiveSetSqp` path. Each knob is forwarded only when
2434/// the user explicitly set it (the `true` flag), so the `pounce_qp`
2435/// defaults stand otherwise.
2436fn apply_qp_subproblem_options(options: &OptionsList, opts: &mut pounce_qp::QpOptions) {
2437    use pounce_qp::AntiCyclingChoice;
2438
2439    if let Ok((v, true)) = options.get_integer_value("sqp_qp_max_iter", "") {
2440        if v >= 0 {
2441            opts.max_iter = v as u32;
2442        }
2443    }
2444    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_feas_tol", "") {
2445        opts.feas_tol = v;
2446    }
2447    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_opt_tol", "") {
2448        opts.opt_tol = v;
2449    }
2450    if let Ok((v, true)) = options.get_numeric_value("sqp_qp_elastic_gamma", "") {
2451        opts.elastic_gamma = v;
2452    }
2453    if let Ok((s, true)) = options.get_string_value("sqp_qp_anti_cycling", "") {
2454        opts.anti_cycling = match s.as_str() {
2455            "expand" => AntiCyclingChoice::Expand,
2456            "bland" => AntiCyclingChoice::Bland,
2457            "none" => AntiCyclingChoice::None,
2458            _ => opts.anti_cycling,
2459        };
2460    }
2461}
2462
2463/// SQP-side analog of [`finalize_via_orig_nlp`]. Hands the SQP
2464/// solution iterate to the user TNLP via the standard
2465/// `finalize_solution` callback. Multiplier lifting goes through
2466/// the same OrigIpoptNlp hooks so the user sees the same shape
2467/// regardless of which algorithm produced the iterate.
2468///
2469/// Returns the user-space objective value on success.
2470fn finalize_via_sqp(
2471    nlp: &Rc<RefCell<dyn IpoptNlp>>,
2472    res: &crate::sqp::SqpResult,
2473    solver_status: pounce_nlp::SolverReturn,
2474    tnlp: &Rc<RefCell<dyn TNLP>>,
2475) -> Result<Number, ()> {
2476    use pounce_linalg::dense_vector::DenseVectorSpace;
2477
2478    let info = tnlp.borrow_mut().get_nlp_info().ok_or(())?;
2479    let n = info.n as usize;
2480    let m = info.m as usize;
2481
2482    // Wrap SQP slices in DenseVectors so we can pass them through
2483    // the OrigIpoptNlp lift_x_to_full / pack_*_for_user hooks.
2484    let nlp_borrow = nlp.borrow();
2485    let n_alg = nlp_borrow.n() as usize;
2486    let m_eq = nlp_borrow.m_eq() as usize;
2487    let m_ineq = nlp_borrow.m_ineq() as usize;
2488    debug_assert_eq!(res.x.len(), n_alg);
2489    debug_assert_eq!(res.lambda_g.len(), m_eq + m_ineq);
2490    debug_assert_eq!(res.lambda_x.len(), n_alg);
2491
2492    let x_space = DenseVectorSpace::new(n_alg as Index);
2493    let c_space = DenseVectorSpace::new(m_eq as Index);
2494    let d_space = DenseVectorSpace::new(m_ineq as Index);
2495
2496    let mut x_dv = x_space.make_new_dense();
2497    x_dv.set_values(&res.x);
2498    let x_vec: Vec<Number> = nlp_borrow.lift_x_to_full(&x_dv);
2499    debug_assert_eq!(x_vec.len(), n);
2500
2501    // λ_x is packed signed (z_l − z_u). Split for lift.
2502    let mut z_l_compressed = x_space.make_new_dense();
2503    let mut z_u_compressed = x_space.make_new_dense();
2504    let zl_vals: Vec<Number> = res.lambda_x.iter().map(|v| v.max(0.0)).collect();
2505    let zu_vals: Vec<Number> = res.lambda_x.iter().map(|v| (-v).max(0.0)).collect();
2506    z_l_compressed.set_values(&zl_vals);
2507    z_u_compressed.set_values(&zu_vals);
2508    // `finalize_solution_*` (not `pack_*`): report unscaled-Lagrangian duals,
2509    // dividing out `obj_scale_factor` — see `finalize_via_orig_nlp` (F1).
2510    let mut z_l = nlp_borrow.finalize_solution_z_l(&z_l_compressed);
2511    if z_l.is_empty() {
2512        z_l = vec![0.0; n];
2513    }
2514    let mut z_u = nlp_borrow.finalize_solution_z_u(&z_u_compressed);
2515    if z_u.is_empty() {
2516        z_u = vec![0.0; n];
2517    }
2518
2519    // λ_g is [y_c; y_d]; split into the c/d blocks for lift.
2520    let mut y_c_dv = c_space.make_new_dense();
2521    let mut y_d_dv = d_space.make_new_dense();
2522    if m_eq > 0 {
2523        y_c_dv.set_values(&res.lambda_g[..m_eq]);
2524    }
2525    if m_ineq > 0 {
2526        y_d_dv.set_values(&res.lambda_g[m_eq..]);
2527    }
2528    let mut lambda = nlp_borrow.finalize_solution_lambda(&y_c_dv, &y_d_dv);
2529    if lambda.is_empty() {
2530        lambda = vec![0.0; m];
2531    }
2532    drop(nlp_borrow);
2533
2534    let mut g_final = vec![0.0; m];
2535    let _ = tnlp.borrow_mut().eval_g(&x_vec, true, &mut g_final);
2536    let f_final = tnlp
2537        .borrow_mut()
2538        .eval_f(&x_vec, true)
2539        .unwrap_or(Number::NAN);
2540    tnlp.borrow_mut().finalize_solution(
2541        pounce_nlp::tnlp::Solution {
2542            status: solver_status,
2543            x: &x_vec,
2544            z_l: &z_l,
2545            z_u: &z_u,
2546            g: &g_final,
2547            lambda: &lambda,
2548            obj_value: f_final,
2549        },
2550        &TnlpIpoptData::default(),
2551        &TnlpIpoptCq::default(),
2552    );
2553    Ok(f_final)
2554}
2555
2556#[cfg(test)]
2557mod tests {
2558    use super::*;
2559    use pounce_nlp::tnlp::{
2560        BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest,
2561        StartingPoint,
2562    };
2563
2564    struct Hs071Stub;
2565    impl TNLP for Hs071Stub {
2566        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
2567            // HS071 dimensions: n=4, m=2, dense Jacobian (8 nz),
2568            // dense lower-triangular Hessian (10 nz).
2569            Some(NlpInfo {
2570                n: 4,
2571                m: 2,
2572                nnz_jac_g: 8,
2573                nnz_h_lag: 10,
2574                index_style: IndexStyle::C,
2575            })
2576        }
2577        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
2578            b.x_l.copy_from_slice(&[1.0; 4]);
2579            b.x_u.copy_from_slice(&[5.0; 4]);
2580            b.g_l.copy_from_slice(&[25.0, 40.0]);
2581            b.g_u.copy_from_slice(&[2.0e19, 40.0]);
2582            true
2583        }
2584        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
2585            sp.x.copy_from_slice(&[1.0, 5.0, 5.0, 1.0]);
2586            true
2587        }
2588        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
2589            Some(x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2])
2590        }
2591        fn eval_grad_f(&mut self, _x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
2592            grad.fill(0.0);
2593            true
2594        }
2595        fn eval_g(&mut self, _x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
2596            g.fill(0.0);
2597            true
2598        }
2599        fn eval_jac_g(
2600            &mut self,
2601            _x: Option<&[Number]>,
2602            _new_x: bool,
2603            mode: SparsityRequest<'_>,
2604        ) -> bool {
2605            if let SparsityRequest::Structure { irow, jcol } = mode {
2606                irow.copy_from_slice(&[0, 0, 0, 0, 1, 1, 1, 1]);
2607                jcol.copy_from_slice(&[0, 1, 2, 3, 0, 1, 2, 3]);
2608            }
2609            true
2610        }
2611        fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
2612    }
2613
2614    #[test]
2615    fn application_default_does_not_select_sqp() {
2616        let mut app = IpoptApplication::new();
2617        app.initialize().unwrap();
2618        assert!(!app.is_sqp_algorithm_selected());
2619    }
2620
2621    #[test]
2622    fn application_routes_to_sqp_when_algorithm_option_set() {
2623        let mut app = IpoptApplication::new();
2624        app.initialize().unwrap();
2625        app.initialize_with_options_str("algorithm active-set-sqp\n")
2626            .unwrap();
2627        assert!(app.is_sqp_algorithm_selected());
2628    }
2629
2630    /// Convex equality NLP fixture for end-to-end SQP testing
2631    /// through `IpoptApplication`:
2632    ///
2633    ///     min ½(x₁² + x₂²) − x₁ − 2x₂  s.t.  x₁ + x₂ = 1
2634    ///
2635    /// Closed form: x* = (0, 1), obj = -1.5, λ_g = 1.
2636    struct ConvexEqTnlp {
2637        finalize_called: std::rc::Rc<std::cell::RefCell<Option<(Vec<Number>, Number)>>>,
2638    }
2639    impl TNLP for ConvexEqTnlp {
2640        fn get_nlp_info(&mut self) -> Option<NlpInfo> {
2641            Some(NlpInfo {
2642                n: 2,
2643                m: 1,
2644                nnz_jac_g: 2,
2645                nnz_h_lag: 2,
2646                index_style: IndexStyle::C,
2647            })
2648        }
2649        fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
2650            b.x_l.copy_from_slice(&[-2.0e19; 2]);
2651            b.x_u.copy_from_slice(&[2.0e19; 2]);
2652            b.g_l.copy_from_slice(&[1.0]);
2653            b.g_u.copy_from_slice(&[1.0]);
2654            true
2655        }
2656        fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
2657            sp.x.copy_from_slice(&[0.0, 0.0]);
2658            true
2659        }
2660        fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
2661            Some(0.5 * (x[0] * x[0] + x[1] * x[1]) - x[0] - 2.0 * x[1])
2662        }
2663        fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, grad: &mut [Number]) -> bool {
2664            grad[0] = x[0] - 1.0;
2665            grad[1] = x[1] - 2.0;
2666            true
2667        }
2668        fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
2669            g[0] = x[0] + x[1];
2670            true
2671        }
2672        fn eval_jac_g(
2673            &mut self,
2674            _x: Option<&[Number]>,
2675            _new_x: bool,
2676            mode: SparsityRequest<'_>,
2677        ) -> bool {
2678            match mode {
2679                SparsityRequest::Structure { irow, jcol } => {
2680                    irow.copy_from_slice(&[0, 0]);
2681                    jcol.copy_from_slice(&[0, 1]);
2682                }
2683                SparsityRequest::Values { values, .. } => {
2684                    values.copy_from_slice(&[1.0, 1.0]);
2685                }
2686            }
2687            true
2688        }
2689        fn eval_h(
2690            &mut self,
2691            _x: Option<&[Number]>,
2692            _new_x: bool,
2693            _obj_factor: Number,
2694            _lambda: Option<&[Number]>,
2695            _new_lambda: bool,
2696            mode: SparsityRequest<'_>,
2697        ) -> bool {
2698            match mode {
2699                SparsityRequest::Structure { irow, jcol } => {
2700                    irow.copy_from_slice(&[0, 1]);
2701                    jcol.copy_from_slice(&[0, 1]);
2702                }
2703                SparsityRequest::Values { values, .. } => {
2704                    values.copy_from_slice(&[1.0, 1.0]);
2705                }
2706            }
2707            true
2708        }
2709        fn finalize_solution(&mut self, sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {
2710            *self.finalize_called.borrow_mut() = Some((sol.x.to_vec(), sol.obj_value));
2711        }
2712    }
2713
2714    #[test]
2715    fn application_sqp_path_solves_convex_eq_nlp_and_finalizes() {
2716        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2717        let tnlp = std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2718            finalize_called: std::rc::Rc::clone(&finalize_slot),
2719        }));
2720
2721        let mut app = IpoptApplication::new();
2722        app.initialize().unwrap();
2723        app.initialize_with_options_str("algorithm active-set-sqp\n")
2724            .unwrap();
2725        let status = app.optimize_tnlp(tnlp);
2726        assert_eq!(status, ApplicationReturnStatus::SolveSucceeded);
2727
2728        // The TNLP's finalize_solution must have been invoked.
2729        let recv = finalize_slot.borrow().clone();
2730        let (x_recv, obj_recv) = recv.expect("finalize_solution was not called");
2731        assert_eq!(x_recv.len(), 2);
2732        assert!((x_recv[0] - 0.0).abs() < 1e-6, "x[0] = {}", x_recv[0]);
2733        assert!((x_recv[1] - 1.0).abs() < 1e-6, "x[1] = {}", x_recv[1]);
2734        assert!(
2735            (obj_recv - (-1.5)).abs() < 1e-6,
2736            "obj = {} but expected -1.5",
2737            obj_recv
2738        );
2739    }
2740
2741    #[test]
2742    fn application_routes_to_sqp_case_insensitively() {
2743        let mut app = IpoptApplication::new();
2744        app.initialize().unwrap();
2745        app.initialize_with_options_str("algorithm Active-Set-SQP\n")
2746            .unwrap();
2747        // get_string_value may return the value as-stored (no
2748        // normalization); the dispatch must handle case
2749        // insensitively per the c11 design choice.
2750        assert!(app.is_sqp_algorithm_selected());
2751    }
2752
2753    #[test]
2754    fn application_constructs_and_loads_options() {
2755        let mut app = IpoptApplication::new();
2756        app.initialize().unwrap();
2757        // ipopt.opt-style file: an integer-typed option registered by
2758        // the Interfaces layer.
2759        app.initialize_with_options_str("print_level 5\nfile_print_level 7\n")
2760            .unwrap();
2761        let (level, found) = app.options().get_integer_value("print_level", "").unwrap();
2762        assert!(found);
2763        assert_eq!(level, 5);
2764    }
2765
2766    #[test]
2767    fn application_sqp_suboptions_propagate_to_builder() {
2768        // All SQP suboptions are read by algorithm_builder_snapshot
2769        // and baked into the builder's `sqp` field.
2770        let mut app = IpoptApplication::new();
2771        app.initialize().unwrap();
2772        app.initialize_with_options_str(
2773            "algorithm active-set-sqp\n\
2774             sqp_globalization l1-elastic\n\
2775             sqp_hessian lbfgs\n\
2776             sqp_max_iter 17\n\
2777             sqp_tol 1e-7\n\
2778             sqp_constr_viol_tol 1e-5\n\
2779             sqp_dual_inf_tol 1e-3\n\
2780             sqp_l1_penalty 2.5\n\
2781             sqp_bt_reduction 0.25\n\
2782             sqp_bt_min_alpha 1e-10\n\
2783             sqp_print_level 2\n\
2784             sqp_lbfgs_max_history 12\n",
2785        )
2786        .unwrap();
2787        let snap = app.algorithm_builder_snapshot();
2788        assert_eq!(
2789            snap.sqp.globalization,
2790            crate::sqp::SqpGlobalization::L1Elastic
2791        );
2792        assert_eq!(snap.sqp.hessian, crate::sqp::SqpHessianSource::Lbfgs);
2793        assert_eq!(snap.sqp.max_iter, 17);
2794        assert!((snap.sqp.tol - 1e-7).abs() < 1e-18);
2795        assert!((snap.sqp.constr_viol_tol - 1e-5).abs() < 1e-18);
2796        assert!((snap.sqp.dual_inf_tol - 1e-3).abs() < 1e-18);
2797        assert!((snap.sqp.l1_penalty - 2.5).abs() < 1e-18);
2798        assert!((snap.sqp.bt_reduction - 0.25).abs() < 1e-18);
2799        assert!((snap.sqp.bt_min_alpha - 1e-10).abs() < 1e-18);
2800        assert_eq!(snap.sqp.print_level, 2);
2801        assert_eq!(snap.sqp.lbfgs_max_history, 12);
2802    }
2803
2804    #[test]
2805    fn application_limited_memory_options_propagate_to_builder() {
2806        use crate::hess::lim_mem_quasi_newton::UpdateType;
2807
2808        // Default: no options set -> bit-exact with Ipopt's default
2809        // (bfgs, history 6). This is what the IPM path runs unless the
2810        // user opts in, so it must not drift.
2811        let mut app = IpoptApplication::new();
2812        app.initialize().unwrap();
2813        let def = app.algorithm_builder_from_options();
2814        assert_eq!(def.limited_memory_update_type, UpdateType::Bfgs);
2815        assert_eq!(def.limited_memory_max_history, 6);
2816
2817        // `limited_memory_update_type=sr1` and a custom history length
2818        // must reach the builder (these were registered upstream but
2819        // read nowhere on the IPM path before — see #131). Honoring
2820        // them is what lets SR1 break the monotone L-BFGS stall.
2821        let mut app = IpoptApplication::new();
2822        app.initialize().unwrap();
2823        app.initialize_with_options_str(
2824            "hessian_approximation limited-memory\n\
2825             limited_memory_update_type sr1\n\
2826             limited_memory_max_history 9\n",
2827        )
2828        .unwrap();
2829        let snap = app.algorithm_builder_from_options();
2830        assert_eq!(snap.limited_memory_update_type, UpdateType::Sr1);
2831        assert_eq!(snap.limited_memory_max_history, 9);
2832    }
2833
2834    #[test]
2835    fn application_sqp_warm_start_round_trip() {
2836        // Drive the convex-equality TNLP through the SQP path
2837        // twice. The first solve produces a working set; the
2838        // second is warm-started from it. The second must converge
2839        // with zero QP solves (the first KKT check declares
2840        // optimality immediately).
2841        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2842        let tnlp_rc: std::rc::Rc<std::cell::RefCell<dyn TNLP>> =
2843            std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2844                finalize_called: std::rc::Rc::clone(&finalize_slot),
2845            }));
2846
2847        let mut app = IpoptApplication::new();
2848        app.initialize().unwrap();
2849        app.initialize_with_options_str("algorithm active-set-sqp\n")
2850            .unwrap();
2851
2852        // Cold solve.
2853        let status_a = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2854        assert_eq!(status_a, ApplicationReturnStatus::SolveSucceeded);
2855        let ws = app.last_sqp_working_set().cloned();
2856        assert!(ws.is_some(), "cold solve must yield a working set");
2857
2858        // Build the warm-start iterate from the converged finalize
2859        // payload (just x; pad multipliers to 0 since the test
2860        // problem is convex).
2861        let (x_recv, _) = finalize_slot.borrow().clone().unwrap();
2862        let warm = crate::sqp::SqpIterates {
2863            x: x_recv,
2864            lambda_g: vec![1.0],
2865            lambda_x: vec![0.0, 0.0],
2866            working: ws,
2867        };
2868        app.set_sqp_warm_start(warm);
2869
2870        // Warm solve.
2871        let status_b = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2872        assert_eq!(status_b, ApplicationReturnStatus::SolveSucceeded);
2873        assert!(app.last_sqp_working_set().is_some());
2874    }
2875
2876    #[test]
2877    fn application_sqp_warm_start_auto_clears_after_use() {
2878        let finalize_slot = std::rc::Rc::new(std::cell::RefCell::new(None));
2879        let tnlp_rc: std::rc::Rc<std::cell::RefCell<dyn TNLP>> =
2880            std::rc::Rc::new(std::cell::RefCell::new(ConvexEqTnlp {
2881                finalize_called: std::rc::Rc::clone(&finalize_slot),
2882            }));
2883        let mut app = IpoptApplication::new();
2884        app.initialize().unwrap();
2885        app.initialize_with_options_str("algorithm active-set-sqp\n")
2886            .unwrap();
2887        app.set_sqp_warm_start(crate::sqp::SqpIterates {
2888            x: vec![0.0, 1.0],
2889            lambda_g: vec![1.0],
2890            lambda_x: vec![0.0, 0.0],
2891            working: None,
2892        });
2893        assert!(app.sqp_warm_start.is_some());
2894        let _ = app.optimize_tnlp(std::rc::Rc::clone(&tnlp_rc));
2895        assert!(
2896            app.sqp_warm_start.is_none(),
2897            "warm-start input must be auto-cleared after use"
2898        );
2899    }
2900
2901    #[test]
2902    fn application_sqp_suboptions_default_when_unset() {
2903        // Without any sqp_* settings, the snapshot should equal
2904        // SqpOptions::default().
2905        let mut app = IpoptApplication::new();
2906        app.initialize().unwrap();
2907        let snap = app.algorithm_builder_snapshot();
2908        let d = crate::sqp::SqpOptions::default();
2909        assert_eq!(snap.sqp.globalization, d.globalization);
2910        assert_eq!(snap.sqp.hessian, d.hessian);
2911        assert_eq!(snap.sqp.max_iter, d.max_iter);
2912        assert!((snap.sqp.tol - d.tol).abs() < 1e-18);
2913        assert!((snap.sqp.constr_viol_tol - d.constr_viol_tol).abs() < 1e-18);
2914        assert!((snap.sqp.dual_inf_tol - d.dual_inf_tol).abs() < 1e-18);
2915        assert!((snap.sqp.l1_penalty - d.l1_penalty).abs() < 1e-18);
2916        assert!((snap.sqp.bt_reduction - d.bt_reduction).abs() < 1e-18);
2917        assert!((snap.sqp.bt_min_alpha - d.bt_min_alpha).abs() < 1e-18);
2918        assert_eq!(snap.sqp.print_level, d.print_level);
2919        assert_eq!(snap.sqp.lbfgs_max_history, d.lbfgs_max_history);
2920    }
2921
2922    #[test]
2923    fn application_reports_problem_dimensions() {
2924        let app = IpoptApplication::new();
2925        let mut tnlp = Hs071Stub;
2926        let info = app.problem_dimensions(&mut tnlp).unwrap();
2927        assert_eq!(info.n, 4);
2928        assert_eq!(info.m, 2);
2929        assert_eq!(info.nnz_jac_g, 8);
2930        assert_eq!(info.nnz_h_lag, 10);
2931    }
2932}