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