pounce_algorithm/alg_builder.rs
1//! Algorithm builder — port of `Algorithm/IpAlgBuilder.{hpp,cpp}`.
2//!
3//! Reads `OptionsList`, walks the dependency order documented in
4//! `ref/Ipopt/AGENT_REFERENCE/ARCHITECTURE.md` §"BuildBasicAlgorithm",
5//! and assembles the strategy objects needed by `IpoptAlgorithm`:
6//!
7//! * `SymLinearSolver` (MA57 / MUMPS / FERAL) → `AugSystemSolver`
8//! (`StdAugSystemSolver`) → `PdSystemSolver` (`PdFullSpaceSolver`)
9//! → `SearchDirCalculator` (`PdSearchDirCalc`).
10//! * `BacktrackingLsAcceptor` (filter / penalty / cg-penalty) →
11//! `BacktrackingLineSearch`.
12//! * `MuUpdate` (monotone / adaptive[+oracle]).
13//! * `ConvCheck` (`OptErrorConvCheck`).
14//! * `IterateInitializer` (default / warm-start) and
15//! `EqMultCalculator` (`LeastSquareMults`).
16//! * `HessianUpdater` (exact / limited-memory).
17//! * `IterationOutput` (`OrigIterationOutput`).
18//! * `NLPScalingObject` (none / user / gradient-based / equilibration-based).
19//!
20//! Phase 7 ships the option-driven dispatch surface; the assembled
21//! `IpoptAlgorithm` lands once each strategy's arithmetic does.
22
23use crate::conv_check::opt_error::OptErrorConvCheck;
24use crate::eq_mult::least_square::LeastSquareMults;
25use crate::hess::exact::ExactHessianUpdater;
26use crate::hess::lim_mem_quasi_newton::{InitialApprox, LimMemQuasiNewtonUpdater, UpdateType};
27use crate::init::default::DefaultIterateInitializer;
28use crate::init::warm_start::WarmStartIterateInitializer;
29use crate::kkt::aug_system_solver::AugSystemSolver;
30use crate::kkt::low_rank_aug_system_solver::LowRankAugSystemSolver;
31use crate::kkt::pd_full_space_solver::PdFullSpaceSolver;
32use crate::kkt::pd_search_dir_calc::PdSearchDirCalc;
33use crate::kkt::perturbation_handler::PdPerturbationHandler;
34use crate::kkt::std_aug_system_solver::StdAugSystemSolver;
35use crate::line_search::backtracking::BacktrackingLineSearch;
36use crate::line_search::filter_acceptor::FilterLsAcceptor;
37use crate::line_search::ls_acceptor::BacktrackingLsAcceptor;
38use crate::line_search::penalty_acceptor::PenaltyLsAcceptor;
39use crate::mu::adaptive::{AdaptiveMuUpdate, MuOracleKind};
40use crate::mu::monotone::MonotoneMuUpdate;
41use crate::output::orig::OrigIterationOutput;
42use pounce_common::types::{Index, Number};
43use pounce_linsol::{SparseSymLinearSolverInterface, TSymLinearSolver};
44use std::cell::RefCell;
45use std::rc::Rc;
46
47/// Backend factory — the application supplies one before calling
48/// [`AlgorithmBuilder::build`]. Mirrors upstream's
49/// `SymLinearSolverFactory` knob in `IpAlgBuilder.cpp`. The default
50/// factory wires in FERAL; MA57 is selectable when the `ma57` cargo
51/// feature is enabled.
52pub type LinearBackendFactory =
53 Box<dyn FnMut(LinearSolverChoice) -> Box<dyn SparseSymLinearSolverInterface>>;
54
55/// Top-level algorithm choice. `InteriorPoint` is pounce's default
56/// (the existing `IpoptAlgorithm`); `ActiveSetSqp` is the
57/// Phase 5b SQP driver in `crate::sqp::SqpAlgorithm`, which uses
58/// `pounce-qp` for QP subproblem solves and reuses
59/// `FilterLsAcceptor` for globalization.
60#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
61pub enum AlgorithmChoice {
62 #[default]
63 InteriorPoint,
64 ActiveSetSqp,
65}
66
67#[derive(Debug, Clone, Copy, PartialEq, Eq)]
68pub enum LinearSolverChoice {
69 Ma57,
70 Feral,
71}
72
73/// Symmetric scaling method applied to the augmented KKT system by
74/// [`TSymLinearSolver`]. Mirrors the `linear_system_scaling` option
75/// in `IpAlgBuilder.cpp:302-318` and the `RuizTSymScalingMethod` /
76/// `Mc19TSymScalingMethod` strategies in upstream Ipopt.
77///
78/// * `None` (default) — no scaling; `TSymLinearSolver` runs with a
79/// null scaling method. Matches upstream's default.
80/// * `Ruiz` — iterative symmetric ∞-norm equilibration (Ruiz, 2001).
81/// Implemented in `pounce_linsol::RuizTSymScalingMethod`.
82/// * `Mc19` — Curtis-Reid (HSL MC19) scaling. Not yet implemented;
83/// falls back to `None` with a warning.
84#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
85pub enum LinearSystemScalingChoice {
86 #[default]
87 None,
88 Ruiz,
89 Mc19,
90 /// `slack-based` — `IpSlackBasedTSymScalingMethod`. Unlike the
91 /// others this one is a function of the iterate, not of the matrix,
92 /// so the algorithm pushes the `s`-block factors down each
93 /// iteration (see `IpoptAlgorithm::push_slack_scaling`).
94 SlackBased,
95}
96
97#[derive(Debug, Clone, Copy, PartialEq, Eq)]
98pub enum MuStrategyChoice {
99 Monotone,
100 Adaptive,
101}
102
103#[derive(Debug, Clone, Copy, PartialEq, Eq)]
104pub enum HessianApproxChoice {
105 Exact,
106 LimitedMemory,
107 /// Partitioned quasi-Newton: one small dense block per element
108 /// function (the objective and each constraint row), assembled into
109 /// a genuine sparse `SymTMatrix`. See
110 /// [`crate::hess::partitioned_quasi_newton`].
111 Partitioned,
112 /// Sparse finite-difference Lagrangian Hessian, recovered by graph
113 /// coloring from the analytic Jacobian. See
114 /// [`crate::hess::fd_hessian`].
115 FiniteDifference,
116}
117
118#[derive(Debug, Clone, Copy, PartialEq, Eq)]
119pub enum LineSearchChoice {
120 Filter,
121 CgPenalty,
122 Penalty,
123}
124
125/// Assembled strategy bundle. Phase 7 ships the structural bundle;
126/// `IpoptAlgorithm::new` reads from this when it lands.
127pub struct AlgorithmBundle {
128 pub mu_update: Box<dyn crate::mu::r#trait::MuUpdate>,
129 pub conv_check: Box<dyn crate::conv_check::r#trait::ConvCheck>,
130 pub init: Box<dyn crate::init::r#trait::IterateInitializer>,
131 pub eq_mult: Box<dyn crate::eq_mult::r#trait::EqMultCalculator>,
132 pub hess: Box<dyn crate::hess::r#trait::HessianUpdater>,
133 pub line_search: BacktrackingLineSearch,
134 pub iter_output: Box<dyn crate::output::r#trait::IterationOutput>,
135 /// `Some` when the builder was given a [`LinearBackendFactory`];
136 /// `None` for the bare structural bundle that pre-Phase-6 unit
137 /// tests still rely on.
138 pub search_dir: Option<PdSearchDirCalc>,
139}
140
141/// Knobs read off `OptionsList` and baked into the assembled
142/// `OptErrorConvCheck`. Defaults mirror
143/// `IpOptErrorConvCheck.cpp:RegisterOptions`.
144#[derive(Debug, Clone)]
145pub struct ConvCheckOptions {
146 pub tol: Number,
147 pub dual_inf_tol: Number,
148 pub constr_viol_tol: Number,
149 pub compl_inf_tol: Number,
150 pub acceptable_tol: Number,
151 pub acceptable_dual_inf_tol: Number,
152 pub acceptable_constr_viol_tol: Number,
153 pub acceptable_compl_inf_tol: Number,
154 pub acceptable_obj_change_tol: Number,
155 pub acceptable_iter: Index,
156 pub max_iter: Index,
157 pub max_cpu_time: Number,
158 pub max_wall_time: Number,
159 pub infeas_stationarity_tol: Number,
160 pub infeas_viol_kappa: Number,
161 pub infeas_max_streak: Index,
162 /// Objective-scale floor below which a strict termination certificate is
163 /// refused while the unscaled KKT error is still above `acceptable_tol`
164 /// (gh #200). `0` disables the mechanism.
165 pub obj_scale_certificate_threshold: Number,
166 /// Safety factor on the per-row floor the **strict** gate uses to decide
167 /// when a constraint residual is finer than the row can represent
168 /// (gh #528). `0` disables the floor, restoring upstream Ipopt's
169 /// bare-absolute primal term.
170 pub primal_noise_floor_kappa: Number,
171 /// Fraction of `acceptable_tol` the KKT error and the objective may drift
172 /// across the acceptable-level streak's window while the streak still
173 /// counts as settled (gh #533). `0` disables the progress test, leaving
174 /// acceptable-level termination the bare consecutive-count criterion.
175 pub acceptable_progress_kappa: Number,
176 /// Safety factor on the scale-relative floor under `dual_inf_tol` the
177 /// **strict** gate judges the dual infeasibility against (gh #532). `0`
178 /// disables the floor, restoring upstream Ipopt's bare-absolute bound.
179 pub dual_inf_scale_kappa: Number,
180}
181
182impl Default for ConvCheckOptions {
183 fn default() -> Self {
184 Self {
185 tol: 1e-8,
186 dual_inf_tol: 1.0,
187 constr_viol_tol: 1e-4,
188 compl_inf_tol: 1e-4,
189 acceptable_tol: 1e-6,
190 acceptable_dual_inf_tol: 1e10,
191 acceptable_constr_viol_tol: 1e-2,
192 acceptable_compl_inf_tol: 1e-2,
193 acceptable_obj_change_tol: 1e20,
194 acceptable_iter: 15,
195 max_iter: 3000,
196 max_cpu_time: 1e6,
197 max_wall_time: 1e6,
198 infeas_stationarity_tol: 1e-8,
199 infeas_viol_kappa: 1e2,
200 infeas_max_streak: 5,
201 obj_scale_certificate_threshold: 1e-4,
202 primal_noise_floor_kappa: 64.0,
203 acceptable_progress_kappa: 1e-1,
204 dual_inf_scale_kappa: 1.0,
205 }
206 }
207}
208
209#[derive(Debug, Clone)]
210pub struct AlgorithmBuilder {
211 /// Top-level algorithm dispatch. Default `InteriorPoint` ⇒
212 /// `build_with_backend` returns the existing `AlgorithmBundle`
213 /// (consumed by `IpoptAlgorithm`). `ActiveSetSqp` ⇒ caller
214 /// must use `build_sqp_with_backend` to assemble the Phase 5b
215 /// `SqpAlgorithm`. The two builder methods sit side by side
216 /// because the assembled algorithm shape differs (IPM bundle
217 /// vs SQP struct).
218 pub algorithm: AlgorithmChoice,
219 pub linear_solver: LinearSolverChoice,
220 /// Symmetric scaling method for the augmented KKT system. Wired
221 /// into [`TSymLinearSolver`] by [`Self::build_with_backend`].
222 /// Mirrors upstream `linear_system_scaling` (`IpAlgBuilder.cpp:538-560`).
223 pub linear_system_scaling: LinearSystemScalingChoice,
224 /// Lazy-vs-eager scaling toggle (`linear_scaling_on_demand`,
225 /// `IpTSymLinearSolver.cpp:50-58`). Only consulted when
226 /// `linear_system_scaling != None`. Upstream default is `true`
227 /// (compute scaling only on the first solve that fails / shows
228 /// poor conditioning); pounce mirrors that. Set to `false` to
229 /// scale every factorization.
230 pub linear_scaling_on_demand: bool,
231 pub mu_strategy: MuStrategyChoice,
232 /// Selector forwarded to [`AdaptiveMuUpdate`] when
233 /// `mu_strategy = Adaptive`. Ignored for `Monotone`. Defaults to
234 /// `QualityFunction` per upstream's `RegisterOptions` default.
235 pub mu_oracle: MuOracleKind,
236 pub hessian_approximation: HessianApproxChoice,
237
238 /// Element update formula for
239 /// [`HessianApproxChoice::Partitioned`] (`partitioned_update_type`).
240 /// SR1 by default: a single constraint is not convex, so damped
241 /// BFGS would force every `∇²c_j` model PSD and then scale it by a
242 /// multiplier of either sign.
243 pub partitioned_update_type: UpdateType,
244 /// Whether the caller named `partitioned_update_type` explicitly, so
245 /// the block mode's BFGS default does not override them.
246 pub partitioned_update_type_was_set: bool,
247 /// Widest element that keeps a dense block under
248 /// [`HessianApproxChoice::Partitioned`]; wider elements degrade to a
249 /// diagonal approximation (`partitioned_max_element`).
250 pub partitioned_max_element: usize,
251 /// Variables the **objective** is nonlinear in, in the compressed
252 /// `x_var` space — `TNLPAdapter::objective_nonlinear_vars`. Consumed
253 /// by both the partitioned updater (as its objective element's
254 /// support) and the finite-difference updater (as the objective's
255 /// contribution to a Jacobian-derived Hessian pattern, which the
256 /// constraint Jacobian cannot supply). `None` leaves each to fall
257 /// back on the first `∇f`'s nonzeros, which is value-derived; see
258 /// that method for what it costs.
259 pub objective_nonlinear_vars: Option<Vec<Index>>,
260 /// `partitioned_curvature_cap` — multiple of an element's implied
261 /// curvature that one update may reach. See
262 /// [`crate::hess::partitioned_quasi_newton`].
263 pub partitioned_curvature_cap: Number,
264 /// How the Lagrangian is split into elements under
265 /// [`HessianApproxChoice::Partitioned`] (`partitioned_elements`).
266 pub partitioned_elements: crate::hess::partitioned_quasi_newton::ElementMode,
267 /// Target primal-block width when `partitioned_elements` is
268 /// `blocks` (`partitioned_block_size`).
269 pub partitioned_block_size: usize,
270 /// Where [`HessianApproxChoice::FiniteDifference`] takes its
271 /// sparsity pattern from (`fd_hessian_pattern`).
272 pub fd_hessian_pattern: crate::hess::fd_hessian::FdPatternSource,
273 /// How finite-difference probe groups are formed
274 /// (`fd_hessian_coloring`).
275 pub fd_hessian_coloring: crate::hess::fd_hessian::FdColoring,
276 /// Relative movement in `x` AND `y` below which the previous Hessian
277 /// is reused (`fd_hessian_reuse_tol`). `0` rebuilds every iteration.
278 pub fd_hessian_reuse_tol: Number,
279 pub limited_memory_update_type: UpdateType,
280 /// History length for the limited-memory quasi-Newton approximation
281 /// (`limited_memory_max_history`). Defaults to upstream's 6.
282 pub limited_memory_max_history: i32,
283 /// `limited_memory_init_val_max` / `_min` — the clamp on the initial
284 /// Hessian scalar σ before the rank-2 updates. Upstream defaults 1e8
285 /// / 1e-8, which `LimMemQuasiNewtonUpdater` has carried as hard-coded
286 /// fields and consumed in `initial_hessian_scalar` all along; only
287 /// the read sites were missing (gh#483, #191 round 2).
288 pub limited_memory_init_val_max: Number,
289 pub limited_memory_init_val_min: Number,
290 /// `limited_memory_initialization` — which formula picks the initial
291 /// Hessian scalar σ. Matches upstream's `scalar1` (σ = sᵀy/sᵀs).
292 /// pounce shipped `scalar2` (σ = yᵀy/sᵀy) with no way to change it,
293 /// because the option was registered and never read (#677).
294 pub limited_memory_initialization: InitialApprox,
295 /// `limited_memory_init_val` — σ on the first iteration, before any
296 /// curvature pair exists, and every iteration under
297 /// `InitialApprox::Constant`. Upstream default 1.0.
298 pub limited_memory_init_val: Number,
299 /// `limited_memory_max_skipping` — consecutive skipped curvature
300 /// updates before the approximation is discarded (#686). Upstream
301 /// default 2.
302 pub limited_memory_max_skipping: Index,
303 /// Positions in the algorithm's compressed `x_var` space that enter
304 /// the problem *nonlinearly* (gh#624). `None` — the default —
305 /// approximates the Hessian over every variable, which is what the
306 /// limited-memory path has always done. When set, the quasi-Newton
307 /// update is restricted to this subspace and the Hessian is exactly
308 /// zero elsewhere. Comes from
309 /// `TNLPAdapter::quasi_newton_nonlinear_vars` (the TNLP's
310 /// `get_list_of_nonlinear_variables`, or the `num_linear_variables`
311 /// prefix fallback) and is ignored on the exact-Hessian path.
312 ///
313 /// The restoration sub-IPM must clear this: the mask indexes the
314 /// original NLP's variables, not the restoration compound primal.
315 pub limited_memory_nonlinear_vars: Option<Vec<Index>>,
316 pub line_search_method: LineSearchChoice,
317 pub warm_start_init_point: bool,
318 /// `mehrotra_algorithm` — when true, [`PdSearchDirCalc`] folds
319 /// the Mehrotra second-order complementarity term into the
320 /// search-direction RHS. Mirrors upstream's
321 /// `IpAlgBuilder.cpp:Mehrotra` flag. Requires `mu_strategy =
322 /// Adaptive` so that an affine step is computed each iteration;
323 /// [`Self::build_with_backend`] does not enforce this — the
324 /// option-parser in `application.rs` is responsible for the
325 /// cascading defaults (`mu_oracle = probing` etc.).
326 pub mehrotra_algorithm: bool,
327 /// `fast_step_computation` — when true, [`PdSearchDirCalc`] accepts
328 /// the search direction without the residual check and allows an
329 /// inexact linear solve. Mirrors upstream's flag of the same name,
330 /// default `no`. The field existed and was consumed from the day the
331 /// search-direction calculator landed, hard-coded to `false`; only
332 /// the option's read site was missing, so setting it did nothing
333 /// (gh#483 follow-up, #191 round 2).
334 pub fast_step_computation: bool,
335 /// `kappa_sigma` — factor bounding how far the bound multipliers may
336 /// deviate from their primal estimates. The clamp
337 /// (`kappa_sigma_clamp`) runs after every accepted step; `< 1`
338 /// disables the correction. Mirrors `IpIpoptAlg.cpp` (Eqn. (16)),
339 /// default `1e10`. Baked onto [`crate::ipopt_alg::IpoptAlgorithm`] by
340 /// the solve path.
341 pub kappa_sigma: Number,
342 /// `recalc_y` / `recalc_y_feas_tol` — least-square re-estimation of
343 /// the equality multipliers once feasible (#677). Registered
344 /// upstream, refused by pounce as unimplemented until now. Default
345 /// `false` matches the registry; the limited-memory path turns it on
346 /// for itself in `application.rs`, as upstream's own option text
347 /// says it does.
348 pub recalc_y: bool,
349 pub recalc_y_feas_tol: Number,
350 /// `kappa_d` — weight of the linear damping term added to the barrier
351 /// objective/gradient (and dual-infeasibility) to handle one-sided
352 /// bounds. Mirrors `IpIpoptCalculatedQuantities.cpp`, default `1e-5`.
353 /// Baked onto [`crate::ipopt_cq::IpoptCalculatedQuantities`] by the
354 /// solve path.
355 pub kappa_d: Number,
356 /// `s_max` — cap on the average multiplier magnitude used to build
357 /// the `(s_d, s_c)` scaling factors of the KKT error test
358 /// (`IpIpoptCalculatedQuantities.cpp:ComputeOptimalityErrorScaling`,
359 /// the paragraph after Eqn. (6) of the implementation paper).
360 /// Registered default `100`, which is what
361 /// [`crate::ipopt_cq::IpoptCalculatedQuantities`] already carries as
362 /// its struct default, so forwarding it is behaviour-neutral for a
363 /// run that does not set it (#551 / #677). Baked onto the cq by the
364 /// solve path, next to `kappa_d`.
365 pub s_max: Number,
366 /// `tiny_step_tol` — relative primal step size below which the full
367 /// step is accepted without line search; repeated tiny steps
368 /// terminate the solve. Mirrors `IpBacktrackingLineSearch.cpp`,
369 /// default `10·EPSILON`. Baked onto
370 /// [`crate::ipopt_alg::IpoptAlgorithm`] by the solve path.
371 pub tiny_step_tol: Number,
372 /// `tiny_step_y_tol` — dual-step threshold; when both primal and dual
373 /// steps are tiny in consecutive iterations the algorithm stops at the
374 /// best attainable accuracy. Default `1e-2`.
375 pub tiny_step_y_tol: Number,
376 /// `diverging_iterates_tol` — if `max_i |x_i|` exceeds this the solve
377 /// aborts as diverging. Default `1e20`.
378 pub diverging_iterates_tol: Number,
379 /// `dual_diverging_streak` (pounce#246) — consecutive growing-dual-
380 /// infeasibility iterations before the dual-divergence guard routes to
381 /// restoration. **Default `0` (off).**
382 ///
383 /// It defaulted to `15` when introduced, on the strength of a reported
384 /// emfl050 bad-warm-start grind. That justification did not survive being
385 /// reproduced: the measurement was caller-side JAX compilation, and the
386 /// build predating the guard solves both emfl050 instances to the same
387 /// optimum in the same time (pounce#246 / pounce#250). What remained was a
388 /// knife-edge, non-monotone effect on four of 1284 MINLPLib models — so it
389 /// is opt-in rather than imposed. See `upstream_options.rs` for the full
390 /// account.
391 pub dual_diverging_streak: Index,
392 /// `dual_divergence_retry_step_tol` (gh#884) — the scale-relative
393 /// step `max_i |d_i| / (1 + |x_i|)` at or below which the biactive
394 /// dual-divergence detector calls the primal iterate *settled*.
395 /// Default `1e-5`; `0` disables the detector without disabling the
396 /// `dual_divergence_retry` option. See `upstream_options.rs` for the
397 /// measured population behind the default.
398 pub dual_divergence_retry_step_tol: Number,
399 /// `dual_divergence_retry_du_floor` (gh#884) — the *unscaled* dual
400 /// infeasibility at or above which the same detector calls the
401 /// multipliers *diverged*. Default `1e2`. Measured in the model's own
402 /// units on purpose: the `s_d`-normalised aggregate is what hid the
403 /// defect. See `upstream_options.rs`.
404 pub dual_divergence_retry_du_floor: Number,
405 /// `resto_decline_deferrals` (gh #534) — how many times the
406 /// acceptable-point restoration decline may be deferred on a solve whose
407 /// NLP error is still contracting. Default `1`; `0` restores the pre-#534
408 /// behaviour (decline immediately, always). See `upstream_options.rs`.
409 pub resto_decline_deferrals: Index,
410 /// `resto_decline_progress_ratio` (gh #534) — required per-iteration
411 /// contraction of the NLP error before a decline is deferred. Default
412 /// `0.5`; at or above `1` the progress requirement is dropped entirely.
413 pub resto_decline_progress_ratio: Number,
414 /// `neg_curv_escapes` (gh #797) — how many times a certified stationary
415 /// point with an indefinite reduced Hessian may be left along a direction
416 /// of negative curvature instead of reported. Default `1`; `0` restores the
417 /// pre-#797 behaviour. See `upstream_options.rs`.
418 pub neg_curv_escapes: Index,
419 /// `limited_memory_ls_failure_restarts` (gh #818) — how many times a
420 /// line-search failure at an already-feasible point may re-anchor the
421 /// quasi-Newton model and retry instead of entering restoration.
422 /// Default `0`, i.e. the rung is off and a line-search failure always
423 /// hands off, which is upstream's behaviour; see
424 /// `DEFAULT_LBFGS_LS_FAILURE_RESTARTS` in `ipopt_alg.rs` for the
425 /// measurement that put it there. See `upstream_options.rs`.
426 pub limited_memory_ls_failure_restarts: Index,
427 /// `kkt_fidelity_tol` (pounce#173). Read by the algorithm as well as by the
428 /// post-solve gate, because the #200 fallback's tiebreak has to rank the two
429 /// candidate points by the status each will be *reported* under. Default
430 /// `0.0` (gate disabled).
431 pub kkt_fidelity_tol: Number,
432 pub conv_check: ConvCheckOptions,
433 pub mu: MuOptions,
434 pub line_search: LineSearchOptions,
435 pub refinement: RefinementOptions,
436 pub perturbation: PerturbationOptions,
437 pub resto: RestoOptions,
438 pub output: OutputOptions,
439 pub warm: WarmStartOptions,
440 /// SQP-specific options (consulted only when
441 /// `algorithm = ActiveSetSqp`).
442 pub sqp: crate::sqp::SqpOptions,
443 /// QP-subproblem-solver options for the active-set SQP path
444 /// (`pounce_qp::QpOptions`), threaded into the `SqpAlgorithm` via
445 /// `with_qp_options`. Consulted only when `algorithm = ActiveSetSqp`.
446 /// Populated from the `sqp_qp_*` CLI options by
447 /// `application::apply_qp_subproblem_options`.
448 pub sqp_qp: pounce_qp::QpOptions,
449 pub init: InitOptions,
450 /// Optional block-triangular / Schur KKT partition (pounce#180 item 2):
451 /// `(schur_indices, feral_cfg)`. When `Some` and the IPM path is selected
452 /// with the feral linear solver and an exact Hessian, `build_with_backend`
453 /// wraps the standard aug-system solver in a
454 /// [`crate::kkt::SchurAugSystemSolver`] over the given KKT-space indices.
455 /// The Schur solver falls back to the standard solver transparently when
456 /// the partition is unsuitable. Set via [`Self::set_kkt_schur`].
457 pub kkt_schur: Option<(Vec<usize>, pounce_feral::FeralConfig)>,
458 /// Shared tally of successful linear-solver quality escalations, handed
459 /// to the assembled
460 /// [`PdFullSpaceSolver`](crate::kkt::pd_full_space_solver::PdFullSpaceSolver)
461 /// by [`Self::build_with_backend`]. `None` leaves that solver with its
462 /// own private counter, which is what every test double and every
463 /// direct builder user gets.
464 ///
465 /// The point of sharing it is the restoration sub-solve: its inner
466 /// algorithm is assembled from a *clone* of this builder
467 /// (`resto_inner_solver::run_inner_resto`), so a `Some` here makes the
468 /// sub-solve's escalations land in the same total as the main loop's.
469 /// gh#857's exact leg escalates once in each, and counting only the
470 /// main loop would report half the trajectory change.
471 pub quality_escalation_counter: Option<Rc<std::cell::Cell<u64>>>,
472}
473
474/// Knobs read off `OptionsList` and baked into
475/// [`DefaultIterateInitializer`]. Defaults mirror
476/// `IpDefaultIterateInitializer.cpp:RegisterOptions`. The Mehrotra
477/// cascade in `application.rs` overrides `bound_push`, `bound_frac`,
478/// and `bound_mult_init_val` to upstream's more-aggressive values
479/// (`10`, `0.2`, `1.0`).
480#[derive(Debug, Clone)]
481pub struct InitOptions {
482 pub bound_push: Number,
483 pub bound_frac: Number,
484 pub slack_bound_push: Number,
485 pub slack_bound_frac: Number,
486 pub constr_mult_init_max: Number,
487 pub bound_mult_init_val: Number,
488 /// `bound_mult_init_method`: `"constant"` (default) or `"mu-based"`
489 /// (matches upstream's `IpDefaultIterateInitializer.cpp`).
490 pub bound_mult_init_method: String,
491 /// `least_square_init_primal` — replace the user's starting `x`
492 /// with the min-norm primal that satisfies the linearized
493 /// constraints. Used by the Mehrotra cascade in `application.rs`
494 /// to drop iter-0 primal infeasibility on LP-shaped problems.
495 /// Mirrors upstream `IpDefaultIterateInitializer.cpp:200-222`.
496 pub least_square_init_primal: bool,
497}
498
499impl Default for InitOptions {
500 fn default() -> Self {
501 Self {
502 bound_push: 1e-2,
503 bound_frac: 1e-2,
504 slack_bound_push: 1e-2,
505 slack_bound_frac: 1e-2,
506 constr_mult_init_max: 1e3,
507 bound_mult_init_val: 1.0,
508 bound_mult_init_method: "constant".into(),
509 least_square_init_primal: false,
510 }
511 }
512}
513
514/// Knobs read off `OptionsList` and baked into
515/// [`WarmStartIterateInitializer`]. Defaults mirror
516/// `IpWarmStartIterateInitializer.cpp:RegisterOptions`.
517///
518/// Wired today: every knob above plus the gh#606 recentering pair.
519/// `warm_start_entire_iterate` / `warm_start_same_structure` are
520/// deliberately *not* here — they name the `GetWarmStartIterate` TNLP
521/// surface pounce does not expose, and are refused by
522/// [`crate::unimplemented_options`] rather than parsed into a field
523/// nothing reads (gh#606).
524#[derive(Debug, Clone)]
525pub struct WarmStartOptions {
526 pub bound_push: Number,
527 pub bound_frac: Number,
528 pub slack_bound_push: Number,
529 pub slack_bound_frac: Number,
530 pub mult_bound_push: Number,
531 pub mult_init_max: Number,
532 pub target_mu: Number,
533 /// The value a NaN-seeded bound multiplier takes: NaN in a
534 /// user-supplied `z`/`v` seed means "unseeded, use the default".
535 /// Threaded from `builder.init.bound_mult_init_val` at build time
536 /// so the Mehrotra override and any user setting stay the single
537 /// source of truth.
538 pub bound_mult_init_val: Number,
539 /// `constr_mult_init_max`, threaded from the init options at build
540 /// time (gh#606). The warm path's reconstruction of an unseeded
541 /// equality-multiplier block runs the *cold* path's least-squares
542 /// solve, so it is capped by the cold path's cap — not by
543 /// `warm_start_mult_init_max`, which caps multipliers a caller
544 /// actually supplied and is three orders looser (1e6 vs 1e3).
545 /// Measured: on `redundant_rows`, whose duplicated equality rows
546 /// make the least-squares system singular, the looser cap let an
547 /// arbitrary estimate through and cost 7 -> 25 iterations.
548 pub constr_mult_init_max: Number,
549 /// `warm_start_recentering` (gh#606). Whether the initializer
550 /// measures the supplied point and adapts μ / the multiplier fills
551 /// to it, or keeps the pre-gh#606 universal constants.
552 pub recentering: WarmStartRecentering,
553}
554
555/// Value of `warm_start_recentering` (gh#606).
556#[derive(Debug, Clone, Copy, PartialEq, Eq)]
557pub enum WarmStartRecentering {
558 /// Pre-gh#606 behaviour: constant pushes and floors, `y` left at
559 /// zero when unseeded, μ untouched unless `warm_start_target_mu`
560 /// is set. The kill switch.
561 None,
562 /// Measure the supplied iterate's residuals and derive μ, the
563 /// bound-multiplier fills, and the equality-multiplier
564 /// reconstruction from them.
565 Residual,
566}
567
568impl WarmStartOptions {
569 /// `mult_init_max` as a usable cap: the registered `0` sentinel
570 /// means "no cap".
571 pub(crate) fn mult_init_max_or_inf(&self) -> Number {
572 if self.mult_init_max > 0.0 {
573 self.mult_init_max
574 } else {
575 Number::INFINITY
576 }
577 }
578}
579
580impl Default for WarmStartOptions {
581 fn default() -> Self {
582 Self {
583 bound_push: 1e-3,
584 bound_frac: 1e-3,
585 slack_bound_push: 1e-3,
586 slack_bound_frac: 1e-3,
587 mult_bound_push: 1e-3,
588 mult_init_max: 1e6,
589 target_mu: 0.0,
590 // seeded from the init options so the default has one
591 // home; build() re-resolves it from the live init options
592 // anyway (see `resolved_warm_options`)
593 bound_mult_init_val: InitOptions::default().bound_mult_init_val,
594 constr_mult_init_max: InitOptions::default().constr_mult_init_max,
595 recentering: WarmStartRecentering::Residual,
596 }
597 }
598}
599
600/// The warm-start options as the initializer actually receives them:
601/// `bound_mult_init_val` and `constr_mult_init_max` come from the
602/// (option-read, Mehrotra-resolved) init options, never from
603/// `WarmStartOptions`'s own copy. Split out of `build()` so the threading is testable.
604pub(crate) fn resolved_warm_options(
605 warm: &WarmStartOptions,
606 init: &InitOptions,
607) -> WarmStartOptions {
608 let mut w = warm.clone();
609 w.bound_mult_init_val = init.bound_mult_init_val;
610 w.constr_mult_init_max = init.constr_mult_init_max;
611 w
612}
613
614/// Knobs read off `OptionsList` and baked into the assembled
615/// `MonotoneMuUpdate` or `AdaptiveMuUpdate`. Defaults mirror
616/// `IpMonotoneMuUpdate.cpp` / `IpAdaptiveMuUpdate.cpp:RegisterOptions`.
617/// `mu_max` defaults to the sentinel `-1`; positive values are baked
618/// into both updaters at build time (adaptive interprets `-1` as
619/// "lazy-init from `mu_max_fact * avrg_compl`").
620#[derive(Debug, Clone)]
621pub struct MuOptions {
622 pub mu_init: Number,
623 pub mu_max: Number,
624 pub mu_max_fact: Number,
625 pub mu_min: Number,
626 pub mu_target: Number,
627 pub mu_linear_decrease_factor: Number,
628 pub mu_superlinear_decrease_power: Number,
629 pub mu_allow_fast_monotone_decrease: bool,
630 pub barrier_tol_factor: Number,
631 /// `tau_min` — floor on the fraction-to-the-boundary parameter
632 /// τ = max(tau_min, 1 − μ). Registered default 0.99, which is what
633 /// both `MonotoneMuUpdate` and `AdaptiveMuUpdate` already carry as
634 /// their struct default, so forwarding it is behaviour-neutral for
635 /// a run that does not set the option (#551 / #677). Consumed by
636 /// both updaters; the adaptive one also uses it in the monotone
637 /// mode and for the post-oracle τ = max(tau_min, 1 − NLP error).
638 pub tau_min: Number,
639 /// `sigma_max` / `sigma_min` — clamp on the centering parameter σ
640 /// chosen by `QualityFunctionMuOracle`. Only consumed when
641 /// `mu_strategy=adaptive` and `mu_oracle=quality-function`.
642 /// Defaults from `IpQualityFunctionMuOracle.cpp:RegisterOptions`.
643 pub sigma_max: Number,
644 pub sigma_min: Number,
645 /// `adaptive_mu_globalization` — globalization strategy for the
646 /// adaptive μ-selection mode. Mirrors
647 /// `IpAdaptiveMuUpdate.cpp:RegisterOptions`. Default is
648 /// `ObjConstrFilter`; the Mehrotra cascade switches to
649 /// `NeverMonotoneMode` to disable globalization entirely.
650 pub adaptive_mu_globalization: crate::mu::adaptive::AdaptiveMuGlobalization,
651 /// `quality_function_norm_type` — norm used inside the quality
652 /// function to aggregate the three KKT components. Forwarded to
653 /// `QualityFunctionMuOracle` when `mu_oracle=quality-function`.
654 pub quality_function_norm_type: crate::mu::oracle::quality_function::NormType,
655 /// `quality_function_centrality` — centrality penalty term added
656 /// to the quality function.
657 pub quality_function_centrality: crate::mu::oracle::quality_function::CentralityType,
658 /// `quality_function_balancing_term` — balancing penalty term in
659 /// the quality function (kicks in when complementarity is far
660 /// below infeasibilities).
661 pub quality_function_balancing_term: crate::mu::oracle::quality_function::BalancingTermType,
662 /// `quality_function_max_section_steps` — cap on golden-section
663 /// iterations when picking σ. Default 8.
664 pub quality_function_max_section_steps: i32,
665 /// `quality_function_section_sigma_tol` — width tolerance in
666 /// σ-space for golden section. Default 1e-2.
667 pub quality_function_section_sigma_tol: Number,
668 /// `quality_function_section_qf_tol` — relative flatness
669 /// tolerance for golden section. Default 0.0.
670 pub quality_function_section_qf_tol: Number,
671 /// `adaptive_mu_safeguard_factor` — guard for the LOQO fallback
672 /// in adaptive mode. Default 0.0.
673 pub adaptive_mu_safeguard_factor: Number,
674 /// `adaptive_mu_monotone_init_factor` — multiplier on the
675 /// average complementarity when seeding monotone mode after a
676 /// free-mode bailout. Default 0.8.
677 pub adaptive_mu_monotone_init_factor: Number,
678 /// `adaptive_mu_restore_previous_iterate` — restore the most
679 /// recent free-mode iterate when switching to fixed mode.
680 /// Default `false`.
681 pub adaptive_mu_restore_previous_iterate: bool,
682 /// `adaptive_mu_max_free_returns` (pounce#749, POUNCE
683 /// extension). Cap on how many times the adaptive strategy may
684 /// leave fixed-mu mode. `-1` = unlimited (upstream behavior).
685 pub adaptive_mu_max_free_returns: i32,
686 /// `adaptive_mu_budget_pin_fraction` (pounce#753, POUNCE
687 /// extension). Fraction of an explicitly-set `max_cpu_time` /
688 /// `max_wall_time` after which the adaptive strategy stops
689 /// returning to free-mu mode and finishes monotone. `1.0`
690 /// disables. Inert unless the caller set a time budget.
691 pub adaptive_mu_budget_pin_fraction: Number,
692 /// `adaptive_mu_kkterror_red_iters` — window length for the
693 /// `KKT_ERROR` globalization history. Default 4.
694 pub adaptive_mu_kkterror_red_iters: usize,
695 /// `adaptive_mu_kkterror_red_fact` — required relative reduction
696 /// of the KKT error over the window. Default 0.9999.
697 pub adaptive_mu_kkterror_red_fact: Number,
698 /// `adaptive_mu_kkt_norm_type` — norm used to score the iterate
699 /// in adaptive globalization decisions.
700 pub adaptive_mu_kkt_norm_type: crate::mu::adaptive::AdaptiveMuKktNorm,
701 /// `filter_margin_fact` — width factor of the margin an entry must
702 /// clear in the `obj-constr-filter` adaptive globalization test
703 /// (`margin = filter_margin_fact * min(filter_max_margin, err)`).
704 /// Only consumed when `mu_strategy=adaptive` and
705 /// `adaptive_mu_globalization=obj-constr-filter` (the default).
706 /// Default 1e-5, from `IpAdaptiveMuUpdate.cpp:RegisterOptions`.
707 pub filter_margin_fact: Number,
708 /// `filter_max_margin` — cap on the margin above. Default 1.0,
709 /// from `IpAdaptiveMuUpdate.cpp:RegisterOptions`.
710 pub filter_max_margin: Number,
711 /// `probing_iterate_quality_factor` (default 1e4, pounce-specific
712 /// — see pounce#58). When the probing (Mehrotra) μ-oracle is
713 /// about to read `curr_avrg_compl()` for its `mu_curr` input, a
714 /// single imbalanced `(s_i, z_i)` pair can inflate the average
715 /// 5+ orders above the stored `data.curr_mu`. The oracle then
716 /// returns `σ · mu_curr` ≫ previous μ, throwing the iterate out
717 /// of the convergence neighborhood. This guard short-circuits
718 /// that case by signalling restoration when the ratio
719 /// `curr_avrg_compl / curr_mu` exceeds the factor. Set to 0 or
720 /// any non-positive value to disable.
721 pub probing_iterate_quality_factor: Number,
722}
723
724impl Default for MuOptions {
725 fn default() -> Self {
726 Self {
727 mu_init: 0.1,
728 mu_max: -1.0,
729 mu_max_fact: 1e3,
730 mu_min: 1e-11,
731 mu_target: 0.0,
732 mu_linear_decrease_factor: 0.2,
733 mu_superlinear_decrease_power: 1.5,
734 mu_allow_fast_monotone_decrease: true,
735 barrier_tol_factor: 10.0,
736 tau_min: 0.99,
737 sigma_max: 1e2,
738 sigma_min: 1e-6,
739 adaptive_mu_globalization:
740 crate::mu::adaptive::AdaptiveMuGlobalization::ObjConstrFilter,
741 quality_function_norm_type:
742 crate::mu::oracle::quality_function::NormType::TwoNormSquared,
743 quality_function_centrality: crate::mu::oracle::quality_function::CentralityType::None,
744 quality_function_balancing_term:
745 crate::mu::oracle::quality_function::BalancingTermType::None,
746 quality_function_max_section_steps: 8,
747 quality_function_section_sigma_tol: 1e-2,
748 quality_function_section_qf_tol: 0.0,
749 adaptive_mu_safeguard_factor: 0.0,
750 adaptive_mu_monotone_init_factor: 0.8,
751 adaptive_mu_restore_previous_iterate: false,
752 adaptive_mu_max_free_returns: -1,
753 adaptive_mu_budget_pin_fraction: 0.75,
754 adaptive_mu_kkterror_red_iters: 4,
755 adaptive_mu_kkterror_red_fact: 0.9999,
756 adaptive_mu_kkt_norm_type: crate::mu::adaptive::AdaptiveMuKktNorm::TwoNormSquared,
757 filter_margin_fact: 1e-5,
758 filter_max_margin: 1.0,
759 probing_iterate_quality_factor: 1e4,
760 }
761 }
762}
763
764/// Knobs baked into the assembled [`BacktrackingLineSearch`]. Defaults
765/// mirror `IpBacktrackingLineSearch.cpp:RegisterOptions`.
766#[derive(Debug, Clone)]
767pub struct LineSearchOptions {
768 /// `alpha_red_factor` — fractional reduction applied to the trial
769 /// step size at every backtracking step
770 /// (`alpha *= alpha_red_factor`). Mirrors upstream's
771 /// `IpBacktrackingLineSearch::alpha_red_factor_`.
772 pub alpha_red_factor: Number,
773 /// `alpha_red_factor_min` — floor on one backtracking reduction,
774 /// which turns the fixed geometric trial sequence into a
775 /// safeguarded quadratic interpolation (gh#818). See
776 /// `BacktrackingLineSearch::next_alpha`.
777 ///
778 /// `None` — the default — means "let the Hessian mode decide", and
779 /// the two modes decide differently:
780 ///
781 /// * **limited-memory → `0.05`** (interpolation on). The
782 /// quasi-Newton model's *scale* can be wrong by orders of
783 /// magnitude in any direction its curvature pairs do not span, so
784 /// the acceptable `alpha` can be far below 1 and the fixed factor
785 /// spends `log(1/alpha)` objective evaluations walking to it.
786 /// * **exact → off** (`alpha_red_factor`, i.e. upstream's fixed
787 /// sequence). A Newton step's length is meaningful, the
788 /// acceptable `alpha` is normally within a few halvings of 1, and
789 /// the trial sequence is not what costs.
790 ///
791 /// The split is measured, not assumed. Forcing
792 /// `alpha_red_factor_min 0.05` onto the exact path moves 3 of the
793 /// 156 fixture-legs in `scripts/sweep-fixtures.sh`, and the one
794 /// that matters is a **status loss**: `eigena2` goes from
795 /// `SolveSucceeded`/27 to `SolvedToAcceptableLevel`/32. The other
796 /// two are `issue_508_infeasible_gap_1em4` 441 → 385 to the same
797 /// certificate — a gain — and an objective digit on
798 /// `hs13_bigstart`. One fixture giving up a solve is the whole
799 /// argument; a faster infeasibility certificate does not buy it
800 /// back. Before `ALPHA_INTERP_MIN_TRIALS` gated the interpolation
801 /// the same experiment moved 9 legs and cost
802 /// `infeasible_square_scaled_1em4` its infeasibility certificate
803 /// (`InfeasibleProblemDetected`/17 → `ErrorInStepComputation`/12);
804 /// the gate narrowed the damage, it did not remove the reason for
805 /// the split.
806 ///
807 /// An explicit `alpha_red_factor_min` from the user is honoured on
808 /// both paths — the mode-dependence is only in the default.
809 pub alpha_red_factor_min: Option<Number>,
810 pub watchdog_shortened_iter_trigger: Index,
811 pub watchdog_trial_iter_max: Index,
812 /// `soft_resto_pderror_reduction_factor` — required relative
813 /// reduction in the primal-dual error for a soft-resto step.
814 /// `0` disables the soft restoration phase.
815 pub soft_resto_pderror_reduction_factor: Number,
816 /// `max_soft_resto_iters` — cap on consecutive soft-resto
817 /// iterations before full restoration is forced.
818 pub max_soft_resto_iters: Index,
819 /// `accept_every_trial_step` — short-circuits the filter / alpha
820 /// loop and accepts the full fraction-to-the-boundary step every
821 /// outer iteration. Mirrors upstream's
822 /// `IpBacktrackingLineSearch::accept_every_trial_step_`. Drops
823 /// global convergence guarantees; only safe for problems where the
824 /// Newton step is already a descent step (LPs, convex QPs). The
825 /// Mehrotra cascade in `application.rs` flips this on.
826 pub accept_every_trial_step: bool,
827 /// `alpha_for_y` — policy for the equality-multiplier (y_c / y_d)
828 /// step length. Upstream default is `Primal`; the Mehrotra cascade
829 /// switches to `BoundMult`.
830 pub alpha_for_y: crate::line_search::backtracking::AlphaForY,
831 /// `accept_after_max_steps` — accept a trial point once this many
832 /// backtracking steps have been taken, even if it fails the
833 /// acceptor's tests. `-1` (the default) disables it, which is why
834 /// wiring it moves no default trajectory. Mirrors
835 /// `IpBacktrackingLineSearch.cpp:759-770`.
836 pub accept_after_max_steps: Index,
837
838 // Filter switching / Armijo / margin constants baked onto the
839 // assembled [`crate::line_search::filter_acceptor::FilterLsAcceptor`]
840 // (only when `line_search_method = Filter`). All were registered but
841 // never read (#191); defaults mirror `IpFilterLSAcceptor.cpp`.
842 /// `eta_phi` — relaxation factor in the Armijo condition (Eqn. (20)).
843 pub eta_phi: Number,
844 /// `delta` — multiplier on the constraint violation in the filter's
845 /// switching rule (Eqn. (19)); maps to
846 /// [`FilterLsAcceptor::delta_armijo`]. Default 1.0, from
847 /// `IpFilterLSAcceptor.cpp:RegisterOptions`.
848 pub delta: Number,
849 /// `theta_min_fact` — constraint-violation threshold factor in the
850 /// switching rule.
851 pub theta_min_fact: Number,
852 /// `theta_max_fact` — upper-bound factor for constraint violation in
853 /// the filter (Eqn. (21)).
854 pub theta_max_fact: Number,
855 /// `theta_max_row_scale_kappa` — multiplier on the constraint-row
856 /// count used as the floor of the `theta_max` reference.
857 /// **Opt-in**: default `0`, which is upstream's bare
858 /// `max(1, theta_0)` floor bit-for-bit. Set to `1` on a large model
859 /// that stalls from a feasible start. See
860 /// [`FilterLsAcceptor::theta_max_row_scale_kappa`].
861 pub theta_max_row_scale_kappa: Number,
862 /// `theta_max_adaptive_trigger` — consecutive line searches whose
863 /// every trial was refused at the `theta_max` gate before the
864 /// ceiling is raised. `0` disables the rule. See
865 /// [`FilterLsAcceptor::theta_max_adaptive_trigger`] (pounce#546).
866 pub theta_max_adaptive_trigger: u32,
867 /// Geometric factor applied to `theta_max` on each adaptive raise.
868 /// See [`FilterLsAcceptor::theta_max_adaptive_factor`].
869 pub theta_max_adaptive_factor: Number,
870 /// Cap on adaptive raises per solve, which is what keeps `theta_max`
871 /// finite. See [`FilterLsAcceptor::theta_max_adaptive_max_raises`].
872 pub theta_max_adaptive_max_raises: u32,
873 /// `gamma_phi` — filter margin factor for the barrier function
874 /// (Eqn. (18a)).
875 pub gamma_phi: Number,
876 /// `gamma_theta` — filter margin factor for the constraint violation
877 /// (Eqn. (18b)).
878 pub gamma_theta: Number,
879 /// `s_phi` — exponent for the linear barrier model in the switching
880 /// rule (Eqn. (19)).
881 pub s_phi: Number,
882 /// `s_theta` — exponent for the current constraint violation in the
883 /// switching rule (Eqn. (19)).
884 pub s_theta: Number,
885 /// `alpha_min_frac` — safety factor for the minimal step size before
886 /// switching to restoration (gamma_alpha, Eqn. (23)).
887 pub alpha_min_frac: Number,
888 /// `obj_max_inc` — max acceptable increase (orders of magnitude) of
889 /// the barrier objective for a trial point.
890 pub obj_max_inc: Number,
891 /// `max_filter_resets` — maximum number of filter resets allowed
892 /// (`0` disables the reset heuristic).
893 pub max_filter_resets: Index,
894 /// `filter_reset_trigger` — successive filter-rejected iterations that
895 /// trigger a filter reset.
896 pub filter_reset_trigger: Index,
897 /// `filter_theta_roundoff_retry` (gh#945) — whether a line search that
898 /// runs out of `alpha` at an iterate already feasible to round-off gets
899 /// one more pass with the filter's `theta` axis measured against
900 /// `theta`'s own evaluation noise, before the driver hands off to a
901 /// restoration phase that has nothing to minimize. **On by default**;
902 /// see `BacktrackingLineSearch::run_filter_line_search` and
903 /// `IpoptCq::theta_evaluation_noise_floor`.
904 pub filter_theta_roundoff_retry: bool,
905
906 // Penalty-acceptor constants baked onto the assembled
907 // [`crate::line_search::penalty_acceptor::PenaltyLsAcceptor`] (only
908 // when `line_search_method = penalty` / `cg-penalty`). Defaults
909 // mirror `IpPenaltyLSAcceptor.cpp:RegisterOptions`.
910 /// `nu_init` — initial value of the penalty parameter ν.
911 pub nu_init: Number,
912 /// `nu_inc` — increment added when ν is bumped.
913 pub nu_inc: Number,
914 /// `rho` — convex-combination weight in the ν update rule.
915 pub rho: Number,
916 /// `eta_penalty` — relaxation factor in the Armijo condition on the
917 /// penalty merit function.
918 pub eta_penalty: Number,
919
920 // Second-order-correction constants baked onto the assembled
921 // [`BacktrackingLineSearch`]. Registered but never read (#191);
922 // defaults mirror `IpBacktrackingLineSearch.cpp`.
923 /// `max_soc` — max second-order-correction trial steps per iteration;
924 /// `0` disables SOC.
925 pub max_soc: Index,
926 /// `kappa_soc` — sufficient-reduction factor for a SOC step to be
927 /// continued.
928 pub kappa_soc: Number,
929 /// `soc_method` — `0` (paper method) or `1` (alpha-on-rhs variant).
930 pub soc_method: Index,
931}
932
933impl Default for LineSearchOptions {
934 fn default() -> Self {
935 Self {
936 alpha_red_factor: 0.5,
937 alpha_red_factor_min: None,
938 watchdog_shortened_iter_trigger: 10,
939 watchdog_trial_iter_max: 3,
940 soft_resto_pderror_reduction_factor: 1.0 - 1e-4,
941 max_soft_resto_iters: 10,
942 accept_every_trial_step: false,
943 alpha_for_y: crate::line_search::backtracking::AlphaForY::Primal,
944 accept_after_max_steps: -1,
945 eta_phi: 1e-8,
946 delta: 1.0,
947 theta_min_fact: 1e-4,
948 theta_max_fact: 1e4,
949 theta_max_row_scale_kappa: 0.0,
950 theta_max_adaptive_trigger: 3,
951 theta_max_adaptive_factor: 100.0,
952 theta_max_adaptive_max_raises: 4,
953 gamma_phi: 1e-8,
954 gamma_theta: 1e-5,
955 s_phi: 2.3,
956 s_theta: 1.1,
957 alpha_min_frac: 0.05,
958 obj_max_inc: 5.0,
959 max_filter_resets: 5,
960 filter_reset_trigger: 5,
961 filter_theta_roundoff_retry: true,
962 nu_init: 1e-6,
963 nu_inc: 1e-4,
964 rho: 0.1,
965 eta_penalty: 1e-8,
966 max_soc: 4,
967 kappa_soc: 0.99,
968 soc_method: 0,
969 }
970 }
971}
972
973/// Inertia-correction / regularization knobs baked onto the assembled
974/// [`crate::kkt::perturbation_handler::PdPerturbationHandler`]. Field
975/// names use the option names; they map to the handler's `delta_xs_*` /
976/// `delta_cd_*` fields. Defaults mirror
977/// `IpPDPerturbationHandler.cpp:RegisterOptions`. All were registered but
978/// never read (#191).
979#[derive(Debug, Clone)]
980pub struct PerturbationOptions {
981 /// `max_hessian_perturbation` → `delta_xs_max`.
982 pub max_hessian_perturbation: Number,
983 /// `min_hessian_perturbation` → `delta_xs_min`.
984 pub min_hessian_perturbation: Number,
985 /// `perturb_inc_fact_first` → `delta_xs_first_inc_fact`.
986 pub perturb_inc_fact_first: Number,
987 /// `perturb_inc_fact` → `delta_xs_inc_fact`.
988 pub perturb_inc_fact: Number,
989 /// `perturb_dec_fact` → `delta_xs_dec_fact`.
990 pub perturb_dec_fact: Number,
991 /// `first_hessian_perturbation` → `delta_xs_init`.
992 pub first_hessian_perturbation: Number,
993 /// `jacobian_regularization_value` → `delta_cd_val`.
994 pub jacobian_regularization_value: Number,
995 /// `jacobian_regularization_exponent` → `delta_cd_exp`.
996 pub jacobian_regularization_exponent: Number,
997 /// `perturb_always_cd` — always regularize the c/d (Jacobian) block.
998 pub perturb_always_cd: bool,
999 /// `perturb_delta_c_max_rungs` → `delta_c_max_rungs` (pounce gh#592).
1000 pub perturb_delta_c_max_rungs: Index,
1001}
1002
1003impl Default for PerturbationOptions {
1004 fn default() -> Self {
1005 Self {
1006 max_hessian_perturbation: 1e20,
1007 min_hessian_perturbation: 1e-20,
1008 perturb_inc_fact_first: 100.0,
1009 perturb_inc_fact: 8.0,
1010 perturb_dec_fact: 1.0 / 3.0,
1011 first_hessian_perturbation: 1e-4,
1012 jacobian_regularization_value: 1e-8,
1013 jacobian_regularization_exponent: 0.25,
1014 perturb_always_cd: false,
1015 perturb_delta_c_max_rungs: 3,
1016 }
1017 }
1018}
1019
1020/// Restoration-phase knobs carried on the outer builder and copied into
1021/// the `RestoAlgorithmBuilder` when the restoration factory is minted
1022/// (`pounce-restoration`). The restoration builder is constructed with
1023/// defaults by each frontend and never options-configured, so these were
1024/// registered but never read (#191). Defaults mirror upstream's
1025/// restoration `RegisterOptions`.
1026#[derive(Debug, Clone)]
1027pub struct RestoOptions {
1028 /// `bound_mult_reset_threshold` — reset bound multipliers to 1 after
1029 /// restoration if the largest exceeds this.
1030 pub bound_mult_reset_threshold: Number,
1031 /// `constr_mult_reset_threshold` — ignore the least-square constraint
1032 /// multiplier estimate after restoration if its norm exceeds this
1033 /// (`0` keeps the estimate).
1034 pub constr_mult_reset_threshold: Number,
1035 /// `resto_penalty_parameter` — penalty on the slack 1-norm in the
1036 /// restoration objective (`rho`).
1037 pub resto_penalty_parameter: Number,
1038 /// `resto_proximity_weight` — proximity-term weight (`eta_factor`;
1039 /// `η = eta_factor · sqrt(μ)`).
1040 pub resto_proximity_weight: Number,
1041 /// `required_infeasibility_reduction` — the restoration sub-solve
1042 /// keeps iterating until the *original* NLP's infeasibility has been
1043 /// reduced to at most this fraction of its value at restoration entry
1044 /// (`κ_resto` in `IpRestoConvCheck.cpp:58`). `0` disables the guard,
1045 /// i.e. restoration runs until the sub-NLP itself converges.
1046 pub required_infeasibility_reduction: Number,
1047 /// `evaluate_orig_obj_at_resto_trial` — evaluate the *original*
1048 /// objective at every restoration trial point, so an iterate the
1049 /// restoration problem likes but the original cannot evaluate is
1050 /// rejected there rather than after the phase exits. Upstream default
1051 /// `yes`. `RestoAlgorithmBuilder` has consumed this since it landed;
1052 /// only the read site was missing (gh#483, #191 round 2).
1053 pub evaluate_orig_obj_at_resto_trial: bool,
1054 /// `expect_infeasible_problem` — enter restoration sooner and demand
1055 /// more infeasibility reduction before leaving it. Upstream default
1056 /// `no`. Same story: consumed, never read.
1057 pub expect_infeasible_problem: bool,
1058 /// `start_with_resto` — switch to restoration in the first iteration.
1059 /// Upstream default `no`. Same story.
1060 pub start_with_resto: bool,
1061 /// `max_resto_iter` — cap on *successive* restoration iterations
1062 /// (`IpRestoConvCheck.cpp:144`'s `maximum_resto_iters`). Consumed by
1063 /// `pounce_restoration::conv_check::RestoConvCheckAdapter`, which
1064 /// returns `MaxIterExceeded` once the count is reached; the value
1065 /// used to be the hard-coded `RESTO_MAX_SUCCESSIVE_ITERS` in
1066 /// `resto_inner_solver.rs`, so setting the option did nothing
1067 /// (#551 / #677). The field is named after the option here, but the
1068 /// consumer's field is `maximum_resto_iters` — which is why grepping
1069 /// for the option name found nothing (#551 caution 2).
1070 ///
1071 /// **This default deliberately differs from the registered one.**
1072 /// `upstream_options.rs` registers Ipopt's `3000000`; pounce has
1073 /// enforced `3000` since the cap landed. Wiring the option must not
1074 /// change what an unset option does, so the effective cap stays
1075 /// `3000` and only an explicit `max_resto_iter` moves it. Raising
1076 /// the default to upstream's number is a trajectory change (it
1077 /// lets a restoration that pounce currently cuts off at 3000 keep
1078 /// going) and belongs to a change that measures it.
1079 pub max_resto_iter: i32,
1080}
1081
1082impl Default for RestoOptions {
1083 fn default() -> Self {
1084 Self {
1085 bound_mult_reset_threshold: 1e3,
1086 constr_mult_reset_threshold: 0.0,
1087 resto_penalty_parameter: 1e3,
1088 resto_proximity_weight: 1.0,
1089 required_infeasibility_reduction: 0.9,
1090 evaluate_orig_obj_at_resto_trial: true,
1091 expect_infeasible_problem: false,
1092 start_with_resto: false,
1093 // NOT the registered default (3000000) — see the field docs.
1094 max_resto_iter: 3000,
1095 }
1096 }
1097}
1098
1099/// Iterative-refinement knobs baked onto the assembled
1100/// [`crate::kkt::pd_full_space_solver::PdFullSpaceSolver`]. Defaults
1101/// mirror `IpPDFullSpaceSolver.cpp:RegisterOptions`. All were registered
1102/// but never read (#191).
1103#[derive(Debug, Clone)]
1104pub struct RefinementOptions {
1105 /// `min_refinement_steps` — minimum iterative-refinement steps per
1106 /// linear solve.
1107 pub min_refinement_steps: Index,
1108 /// `max_refinement_steps` — maximum iterative-refinement steps.
1109 pub max_refinement_steps: Index,
1110 /// `residual_ratio_max` — refine until the residual test ratio drops
1111 /// below this (or `max_refinement_steps` is reached).
1112 pub residual_ratio_max: Number,
1113 /// `residual_ratio_singular` — above this ratio after failed
1114 /// refinement, the system is declared singular.
1115 pub residual_ratio_singular: Number,
1116 /// `residual_improvement_factor` — minimum per-step reduction of the
1117 /// residual test ratio before refinement is aborted.
1118 pub residual_improvement_factor: Number,
1119 /// `neg_curv_test_tol` — tolerance α_n of the inertia-free curvature
1120 /// test of Zavala & Chiang (2014). Zero (the registered default)
1121 /// disables the heuristic and keeps the inertia check; positive
1122 /// turns the inertia check off and accepts the factorization only
1123 /// when the computed direction passes the curvature test in
1124 /// `PdFullSpaceSolver::solve_once`.
1125 pub neg_curv_test_tol: Number,
1126 /// `neg_curv_test_reg` — whether the curvature test includes the
1127 /// primal regularization δ_x‖dx‖² + δ_s‖ds‖². Registered default
1128 /// `yes`; `no` reproduces the original Ipopt form that ignores it.
1129 /// Only consulted when `neg_curv_test_tol > 0`.
1130 pub neg_curv_test_reg: bool,
1131}
1132
1133impl Default for RefinementOptions {
1134 fn default() -> Self {
1135 Self {
1136 min_refinement_steps: 1,
1137 max_refinement_steps: 10,
1138 residual_ratio_max: 1e-10,
1139 residual_ratio_singular: 1e-5,
1140 residual_improvement_factor: 0.999_999_999,
1141 neg_curv_test_tol: 0.0,
1142 neg_curv_test_reg: true,
1143 }
1144 }
1145}
1146
1147/// Knobs baked into the assembled [`OrigIterationOutput`]. Defaults
1148/// mirror `IpOrigIterationOutput.cpp:RegisterOptions` /
1149/// `IpAlgorithmRegOp.cpp`.
1150#[derive(Debug, Clone)]
1151pub struct OutputOptions {
1152 pub print_frequency_iter: Index,
1153 pub print_frequency_time: Number,
1154 /// `print_info_string` (default `false`). When on, the iter row
1155 /// ends with the contents of `IpoptData::info_string` so users
1156 /// can read the per-iteration diagnostic tags.
1157 pub print_info_string: bool,
1158 /// `inf_pr_output` — `"original"` (default) prints the unscaled
1159 /// NLP primal infeasibility; `"internal"` prints the internal
1160 /// reformulated violation. Only meaningful once NLP-side scaling
1161 /// is in play; until then both modes produce the same number.
1162 pub inf_pr_output_internal: bool,
1163}
1164
1165impl Default for OutputOptions {
1166 fn default() -> Self {
1167 Self {
1168 print_frequency_iter: 1,
1169 print_frequency_time: 0.0,
1170 print_info_string: false,
1171 inf_pr_output_internal: false,
1172 }
1173 }
1174}
1175
1176impl Default for AlgorithmBuilder {
1177 fn default() -> Self {
1178 Self {
1179 algorithm: AlgorithmChoice::default(),
1180 linear_solver: LinearSolverChoice::Feral,
1181 linear_system_scaling: LinearSystemScalingChoice::None,
1182 linear_scaling_on_demand: true,
1183 mu_strategy: MuStrategyChoice::Monotone,
1184 mu_oracle: MuOracleKind::QualityFunction,
1185 hessian_approximation: HessianApproxChoice::Exact,
1186 partitioned_update_type: UpdateType::Sr1,
1187 partitioned_update_type_was_set: false,
1188 partitioned_max_element: 64,
1189 objective_nonlinear_vars: None,
1190 partitioned_curvature_cap: Number::INFINITY,
1191 partitioned_elements: crate::hess::partitioned_quasi_newton::ElementMode::PerConstraint,
1192 partitioned_block_size: 64,
1193 fd_hessian_pattern: crate::hess::fd_hessian::FdPatternSource::Declared,
1194 fd_hessian_coloring: crate::hess::fd_hessian::FdColoring::Cpr,
1195 fd_hessian_reuse_tol: 0.0,
1196 limited_memory_update_type: UpdateType::Bfgs,
1197 limited_memory_max_history: 6,
1198 limited_memory_init_val_max: 1e8,
1199 limited_memory_init_val_min: 1e-8,
1200 limited_memory_initialization: InitialApprox::Scalar1,
1201 limited_memory_init_val: 1.0,
1202 limited_memory_max_skipping: 2,
1203 limited_memory_nonlinear_vars: None,
1204 line_search_method: LineSearchChoice::Filter,
1205 warm_start_init_point: false,
1206 mehrotra_algorithm: false,
1207 fast_step_computation: false,
1208 kappa_sigma: 1e10,
1209 recalc_y: false,
1210 recalc_y_feas_tol: 1e-6,
1211 kappa_d: 1e-5,
1212 s_max: 100.0,
1213 tiny_step_tol: 10.0 * Number::EPSILON,
1214 tiny_step_y_tol: 1e-2,
1215 diverging_iterates_tol: 1e20,
1216 dual_divergence_retry_step_tol: 1e-5,
1217 dual_divergence_retry_du_floor: 1e2,
1218 dual_diverging_streak: 0,
1219 resto_decline_deferrals: 1,
1220 resto_decline_progress_ratio: 0.5,
1221 neg_curv_escapes: 1,
1222 limited_memory_ls_failure_restarts: 0,
1223 kkt_fidelity_tol: 0.0,
1224 conv_check: ConvCheckOptions::default(),
1225 mu: MuOptions::default(),
1226 line_search: LineSearchOptions::default(),
1227 refinement: RefinementOptions::default(),
1228 perturbation: PerturbationOptions::default(),
1229 resto: RestoOptions::default(),
1230 output: OutputOptions::default(),
1231 warm: WarmStartOptions::default(),
1232 sqp: crate::sqp::SqpOptions::default(),
1233 sqp_qp: pounce_qp::QpOptions::sqp_subproblem(),
1234 init: InitOptions::default(),
1235 kkt_schur: None,
1236 quality_escalation_counter: None,
1237 }
1238 }
1239}
1240
1241impl AlgorithmBuilder {
1242 pub fn new() -> Self {
1243 Self::default()
1244 }
1245
1246 /// Install a Schur KKT partition (pounce#180 item 2). `schur_indices` are
1247 /// KKT-space indices (`0..dim`, the `x,s,c,d` block order the aug-system
1248 /// solver assembles); `cfg` configures the per-block feral solvers. Only
1249 /// honored on the IPM + feral + exact-Hessian path by
1250 /// [`Self::build_with_backend`]; ignored otherwise.
1251 pub fn set_kkt_schur(&mut self, schur_indices: Vec<usize>, cfg: pounce_feral::FeralConfig) {
1252 self.kkt_schur = Some((schur_indices, cfg));
1253 }
1254
1255 /// Assemble the strategy bundle without a search-direction
1256 /// calculator. Used by structural unit tests that don't want to
1257 /// pull in a linear-solver backend.
1258 pub fn build(&self) -> AlgorithmBundle {
1259 self.build_inner(None)
1260 }
1261
1262 /// Same as [`Self::build`] but also constructs the
1263 /// `SymLinearSolver → AugSystemSolver → PdFullSpaceSolver →
1264 /// PdSearchDirCalc` chain via the supplied `factory`.
1265 pub fn build_with_backend(&self, mut factory: LinearBackendFactory) -> AlgorithmBundle {
1266 let backend = factory(self.linear_solver);
1267 let make_scaling = || -> Option<Box<dyn pounce_linsol::TSymScalingMethod>> {
1268 match self.linear_system_scaling {
1269 LinearSystemScalingChoice::None => None,
1270 LinearSystemScalingChoice::Ruiz => {
1271 Some(Box::new(pounce_linsol::RuizTSymScalingMethod::new()))
1272 }
1273 LinearSystemScalingChoice::Mc19 => {
1274 tracing::warn!(target: "pounce::algorithm",
1275 "pounce: linear_system_scaling=mc19 not yet implemented; using no scaling"
1276 );
1277 None
1278 }
1279 LinearSystemScalingChoice::SlackBased => {
1280 Some(Box::new(pounce_linsol::SlackBasedTSymScalingMethod::new()))
1281 }
1282 }
1283 };
1284 let linsol = TSymLinearSolver::new(backend, make_scaling(), self.linear_scaling_on_demand);
1285 let inner_aug = StdAugSystemSolver::new(linsol);
1286 // Limited-memory mode publishes the Hessian as a
1287 // `LowRankUpdateSymMatrix`; wrap the standard solver in the
1288 // Sherman-Morrison-Woodbury low-rank solver so the augmented
1289 // system factorizes only the diagonal `B0` and the quasi-Newton
1290 // update is applied as a rank-`m` correction (`O(n·m)` memory).
1291 let is_lbfgs = matches!(
1292 self.hessian_approximation,
1293 HessianApproxChoice::LimitedMemory
1294 );
1295 let aug_solver: Box<dyn AugSystemSolver> = if is_lbfgs {
1296 // A second, independent inner solver for the Hessian-free
1297 // solves. Both see one W shape each for their whole life, so
1298 // neither re-runs the backend's symbolic factorization when
1299 // the other's shape comes round — see
1300 // `LowRankAugSystemSolver::with_bypass_solver` (gh#730).
1301 let bypass_linsol = TSymLinearSolver::new(
1302 factory(self.linear_solver),
1303 make_scaling(),
1304 self.linear_scaling_on_demand,
1305 );
1306 Box::new(LowRankAugSystemSolver::with_bypass_solver(
1307 Box::new(inner_aug),
1308 Box::new(StdAugSystemSolver::new(bypass_linsol)),
1309 ))
1310 } else if let Some((indices, cfg)) = self.kkt_schur.clone() {
1311 // Block-triangular / Schur KKT path (pounce#180 item 2). Only on the
1312 // exact-Hessian feral path — the Schur backend is feral-specific,
1313 // and the L-BFGS low-rank Woodbury wrapper owns the (2,2) block.
1314 // The Schur solver falls back to `StdAugSystemSolver` transparently
1315 // when the partition is unsuitable, so a stray hook never breaks a
1316 // solve; we gate on `linear_solver == Feral` here to avoid silently
1317 // ignoring a user's explicit MA57 selection.
1318 if matches!(self.linear_solver, LinearSolverChoice::Feral) {
1319 Box::new(crate::kkt::SchurAugSystemSolver::new(
1320 inner_aug, indices, cfg,
1321 ))
1322 } else {
1323 Box::new(inner_aug)
1324 }
1325 } else {
1326 Box::new(inner_aug)
1327 };
1328 // Inertia-correction / Jacobian-regularization constants (#191):
1329 // registered but previously never read. Defaults equal the
1330 // registered defaults. `perturb_always_cd` goes through the setter
1331 // because it also rebuilds the initial jac-degeneracy state.
1332 let mut ph = PdPerturbationHandler::new();
1333 ph.delta_xs_max = self.perturbation.max_hessian_perturbation;
1334 ph.delta_xs_min = self.perturbation.min_hessian_perturbation;
1335 ph.delta_xs_first_inc_fact = self.perturbation.perturb_inc_fact_first;
1336 ph.delta_xs_inc_fact = self.perturbation.perturb_inc_fact;
1337 ph.delta_xs_dec_fact = self.perturbation.perturb_dec_fact;
1338 ph.delta_xs_init = self.perturbation.first_hessian_perturbation;
1339 ph.delta_cd_val = self.perturbation.jacobian_regularization_value;
1340 ph.delta_cd_exp = self.perturbation.jacobian_regularization_exponent;
1341 ph.set_perturb_always_cd(self.perturbation.perturb_always_cd);
1342 ph.delta_c_max_rungs = self.perturbation.perturb_delta_c_max_rungs;
1343 let perturb = Rc::new(RefCell::new(ph));
1344 let mut pd_solver = PdFullSpaceSolver::new(aug_solver, perturb);
1345 // Iterative-refinement constants (#191): registered but previously
1346 // never read, so overrides were silently dropped. Defaults equal
1347 // the registered defaults.
1348 pd_solver.min_refinement_steps = self.refinement.min_refinement_steps;
1349 pd_solver.max_refinement_steps = self.refinement.max_refinement_steps;
1350 pd_solver.residual_ratio_max = self.refinement.residual_ratio_max;
1351 pd_solver.residual_ratio_singular = self.refinement.residual_ratio_singular;
1352 pd_solver.residual_improvement_factor = self.refinement.residual_improvement_factor;
1353 // Inertia-free curvature test (#551 / #677). Both were registered
1354 // and never read; `neg_curv_test_tol` defaults to 0, which leaves
1355 // the heuristic off and the inertia check on, so this changes
1356 // nothing for a run that does not set it.
1357 pd_solver.neg_curv_test_tol = self.refinement.neg_curv_test_tol;
1358 pd_solver.neg_curv_test_reg = self.refinement.neg_curv_test_reg;
1359 // gh#857: share the escalation tally with the caller, so the
1360 // restoration sub-solve built from a clone of this builder counts
1361 // into the same total.
1362 if let Some(counter) = self.quality_escalation_counter.as_ref() {
1363 pd_solver.set_quality_escalation_counter(Rc::clone(counter));
1364 }
1365 let mut search_dir = PdSearchDirCalc::new(pd_solver);
1366 search_dir.mehrotra_algorithm = self.mehrotra_algorithm;
1367 search_dir.fast_step_computation = self.fast_step_computation;
1368 self.build_inner(Some(search_dir))
1369 }
1370
1371 /// Phase 5b assembly path for the SQP algorithm. Consults
1372 /// `self.algorithm`: when `ActiveSetSqp`, constructs an
1373 /// `SqpAlgorithm` using the supplied backend factory for the
1374 /// QP subproblem solver; otherwise returns `None` so the
1375 /// caller can fall back to the IPM `build_with_backend`.
1376 ///
1377 /// Sister to `build_with_backend`: the SQP algorithm doesn't
1378 /// share `AlgorithmBundle`'s shape (no mu_update / no IPM
1379 /// line search), so the two paths return different types.
1380 pub fn build_sqp_with_backend(
1381 &self,
1382 mut factory: LinearBackendFactory,
1383 ) -> Option<crate::sqp::SqpAlgorithm> {
1384 if !matches!(self.algorithm, AlgorithmChoice::ActiveSetSqp) {
1385 return None;
1386 }
1387 let backend = factory(self.linear_solver);
1388 let qp_solver = pounce_qp::ParametricActiveSetSolver::new(backend);
1389 Some(
1390 crate::sqp::SqpAlgorithm::new(qp_solver, self.sqp.clone())
1391 .with_qp_options(self.sqp_qp.clone()),
1392 )
1393 }
1394
1395 fn build_inner(&self, search_dir: Option<PdSearchDirCalc>) -> AlgorithmBundle {
1396 let mu_update: Box<dyn crate::mu::r#trait::MuUpdate> = match self.mu_strategy {
1397 MuStrategyChoice::Monotone => {
1398 let mut m = MonotoneMuUpdate::new();
1399 m.mu_init = self.mu.mu_init;
1400 // `mu_max` sentinel `-1` keeps the monotone default
1401 // (1e5); only override on a user-supplied positive.
1402 if self.mu.mu_max > 0.0 {
1403 m.mu_max = self.mu.mu_max;
1404 }
1405 m.mu_min = self.mu.mu_min;
1406 m.mu_target = self.mu.mu_target;
1407 m.mu_linear_decrease_factor = self.mu.mu_linear_decrease_factor;
1408 m.mu_superlinear_decrease_power = self.mu.mu_superlinear_decrease_power;
1409 m.mu_allow_fast_monotone_decrease = self.mu.mu_allow_fast_monotone_decrease;
1410 m.barrier_tol_factor = self.mu.barrier_tol_factor;
1411 m.tau_min = self.mu.tau_min;
1412 m.compl_inf_tol = self.conv_check.compl_inf_tol;
1413 Box::new(m)
1414 }
1415 MuStrategyChoice::Adaptive => {
1416 let mut adaptive = AdaptiveMuUpdate::new();
1417 adaptive.mu_oracle = self.mu_oracle;
1418 adaptive.mu_init = self.mu.mu_init;
1419 // Adaptive treats `mu_max == -1` as "lazy init from
1420 // `mu_max_fact * curr_avrg_compl`" — forward the
1421 // sentinel as-is.
1422 adaptive.mu_max = self.mu.mu_max;
1423 adaptive.mu_max_fact = self.mu.mu_max_fact;
1424 adaptive.mu_min = self.mu.mu_min;
1425 adaptive.compl_inf_tol = self.conv_check.compl_inf_tol;
1426 adaptive.mu_linear_decrease_factor = self.mu.mu_linear_decrease_factor;
1427 adaptive.mu_superlinear_decrease_power = self.mu.mu_superlinear_decrease_power;
1428 adaptive.barrier_tol_factor = self.mu.barrier_tol_factor;
1429 adaptive.tau_min = self.mu.tau_min;
1430 adaptive.sigma_min = self.mu.sigma_min;
1431 adaptive.sigma_max = self.mu.sigma_max;
1432 adaptive.adaptive_mu_globalization = self.mu.adaptive_mu_globalization;
1433 adaptive.qf_norm_type = self.mu.quality_function_norm_type;
1434 adaptive.qf_centrality_type = self.mu.quality_function_centrality;
1435 adaptive.qf_balancing_term = self.mu.quality_function_balancing_term;
1436 adaptive.qf_max_section_steps = self.mu.quality_function_max_section_steps;
1437 adaptive.qf_section_sigma_tol = self.mu.quality_function_section_sigma_tol;
1438 adaptive.qf_section_qf_tol = self.mu.quality_function_section_qf_tol;
1439 adaptive.probing_iterate_quality_factor = self.mu.probing_iterate_quality_factor;
1440 adaptive.adaptive_mu_safeguard_factor = self.mu.adaptive_mu_safeguard_factor;
1441 adaptive.adaptive_mu_monotone_init_factor =
1442 self.mu.adaptive_mu_monotone_init_factor;
1443 adaptive.restore_accepted_iterate = self.mu.adaptive_mu_restore_previous_iterate;
1444 adaptive.max_free_returns = self.mu.adaptive_mu_max_free_returns;
1445 adaptive.budget_pin_fraction = self.mu.adaptive_mu_budget_pin_fraction;
1446 // The pin measures against the same budget the
1447 // convergence check enforces (pounce#753); the
1448 // `conv_check` copies are the ones the application
1449 // also hands to `Deadline::new`.
1450 adaptive.max_cpu_time = self.conv_check.max_cpu_time;
1451 adaptive.max_wall_time = self.conv_check.max_wall_time;
1452 adaptive.adaptive_mu_kkterror_red_iters = self.mu.adaptive_mu_kkterror_red_iters;
1453 adaptive.adaptive_mu_kkterror_red_fact = self.mu.adaptive_mu_kkterror_red_fact;
1454 adaptive.adaptive_mu_kkt_norm = self.mu.adaptive_mu_kkt_norm_type;
1455 adaptive.filter_margin_fact = self.mu.filter_margin_fact;
1456 adaptive.filter_max_margin = self.mu.filter_max_margin;
1457 Box::new(adaptive)
1458 }
1459 };
1460
1461 let acceptor: Box<dyn BacktrackingLsAcceptor> = match self.line_search_method {
1462 LineSearchChoice::Filter => {
1463 // Filter switching / Armijo / margin constants (#191):
1464 // registered but previously never read. Set them on the
1465 // concrete acceptor before boxing; defaults equal the
1466 // registered defaults, so a run that doesn't set them is
1467 // unchanged.
1468 let mut f = FilterLsAcceptor::default();
1469 f.eta_phi = self.line_search.eta_phi;
1470 f.delta_armijo = self.line_search.delta;
1471 f.theta_min_fact = self.line_search.theta_min_fact;
1472 f.theta_max_fact = self.line_search.theta_max_fact;
1473 f.theta_max_row_scale_kappa = self.line_search.theta_max_row_scale_kappa;
1474 f.theta_max_adaptive_trigger = self.line_search.theta_max_adaptive_trigger;
1475 f.theta_max_adaptive_factor = self.line_search.theta_max_adaptive_factor;
1476 f.theta_max_adaptive_max_raises = self.line_search.theta_max_adaptive_max_raises;
1477 f.gamma_phi = self.line_search.gamma_phi;
1478 f.gamma_theta = self.line_search.gamma_theta;
1479 f.s_phi = self.line_search.s_phi;
1480 f.s_theta = self.line_search.s_theta;
1481 f.alpha_min_frac = self.line_search.alpha_min_frac;
1482 f.obj_max_inc = self.line_search.obj_max_inc;
1483 f.max_filter_resets = self.line_search.max_filter_resets;
1484 f.filter_reset_trigger = self.line_search.filter_reset_trigger;
1485 Box::new(f)
1486 }
1487 // Penalty-acceptor constants: same direct-field pattern as
1488 // the filter arm above. `reset()` re-seeds ν (and `last_nu`)
1489 // from the freshly-set `nu_init`, which `default()` had
1490 // seeded from the registered default.
1491 LineSearchChoice::Penalty | LineSearchChoice::CgPenalty => {
1492 // CG-penalty acceptor lands with the rest of the
1493 // CG-penalty path; fall back to the penalty acceptor's
1494 // surface for now.
1495 let mut p = PenaltyLsAcceptor::default();
1496 p.nu_init = self.line_search.nu_init;
1497 p.nu_inc = self.line_search.nu_inc;
1498 p.rho = self.line_search.rho;
1499 p.eta_penalty = self.line_search.eta_penalty;
1500 p.reset();
1501 Box::new(p)
1502 }
1503 };
1504 let mut line_search = BacktrackingLineSearch::new(acceptor);
1505 line_search.alpha_red_factor = self.line_search.alpha_red_factor;
1506 line_search.filter_theta_roundoff_retry = self.line_search.filter_theta_roundoff_retry;
1507 // Resolve `None` against the Hessian mode; see the field's doc
1508 // for the measurement behind the split (gh#818).
1509 line_search.alpha_red_factor_min =
1510 self.line_search
1511 .alpha_red_factor_min
1512 .unwrap_or(match self.hessian_approximation {
1513 // The criterion in the field's doc is whether the
1514 // model's *scale* is trustworthy, not whether it is
1515 // the exact Hessian. A quasi-Newton `B` can be wrong
1516 // by orders of magnitude in any direction its
1517 // curvature pairs do not span, which is as true of
1518 // the partitioned elements as of the limited-memory
1519 // ones — they are the same update on a finer
1520 // decomposition.
1521 HessianApproxChoice::LimitedMemory | HessianApproxChoice::Partitioned => 0.05,
1522 // Equal to `alpha_red_factor`, so the clamp in
1523 // `next_alpha` collapses and the sequence is upstream's.
1524 //
1525 // `FiniteDifference` belongs here rather than above:
1526 // it recovers the Lagrangian Hessian itself by
1527 // probing the analytic Jacobian, so it carries no
1528 // curvature history and its step is a Newton step
1529 // whose length is meaningful. It is the same
1530 // distinction `hessian_at_current` draws — a pure
1531 // function of `(x, y)` on one side, a history-carrying
1532 // `B` on the other. And the measurement behind the
1533 // split cuts this way too: forcing `0.05` onto a path
1534 // with a meaningful step length cost `eigena2` a
1535 // solve.
1536 HessianApproxChoice::Exact | HessianApproxChoice::FiniteDifference => {
1537 self.line_search.alpha_red_factor
1538 }
1539 });
1540 line_search.watchdog_shortened_iter_trigger =
1541 self.line_search.watchdog_shortened_iter_trigger;
1542 line_search.watchdog_trial_iter_max = self.line_search.watchdog_trial_iter_max;
1543 line_search.soft_resto_pderror_reduction_factor =
1544 self.line_search.soft_resto_pderror_reduction_factor;
1545 line_search.max_soft_resto_iters = self.line_search.max_soft_resto_iters;
1546 line_search.accept_every_trial_step = self.line_search.accept_every_trial_step;
1547 line_search.alpha_for_y = self.line_search.alpha_for_y;
1548 line_search.accept_after_max_steps = self.line_search.accept_after_max_steps;
1549 // Second-order-correction constants (#191): registered but
1550 // previously never read. Same direct-field pattern as the
1551 // watchdog knobs above.
1552 line_search.max_soc = self.line_search.max_soc;
1553 line_search.kappa_soc = self.line_search.kappa_soc;
1554 line_search.soc_method = self.line_search.soc_method;
1555
1556 let conv_check: Box<dyn crate::conv_check::r#trait::ConvCheck> =
1557 Box::new(OptErrorConvCheck {
1558 tol: self.conv_check.tol,
1559 dual_inf_tol: self.conv_check.dual_inf_tol,
1560 constr_viol_tol: self.conv_check.constr_viol_tol,
1561 compl_inf_tol: self.conv_check.compl_inf_tol,
1562 acceptable_tol: self.conv_check.acceptable_tol,
1563 acceptable_dual_inf_tol: self.conv_check.acceptable_dual_inf_tol,
1564 acceptable_constr_viol_tol: self.conv_check.acceptable_constr_viol_tol,
1565 acceptable_compl_inf_tol: self.conv_check.acceptable_compl_inf_tol,
1566 acceptable_obj_change_tol: self.conv_check.acceptable_obj_change_tol,
1567 acceptable_iter: self.conv_check.acceptable_iter,
1568 max_iter: self.conv_check.max_iter,
1569 max_cpu_time: self.conv_check.max_cpu_time,
1570 max_wall_time: self.conv_check.max_wall_time,
1571 acceptable_count: 0,
1572 last_acceptable_obj: None,
1573 infeas_stationarity_tol: self.conv_check.infeas_stationarity_tol,
1574 infeas_viol_kappa: self.conv_check.infeas_viol_kappa,
1575 infeas_max_streak: self.conv_check.infeas_max_streak,
1576 infeas_streak: 0,
1577 obj_scale_certificate_threshold: self.conv_check.obj_scale_certificate_threshold,
1578 primal_noise_floor_kappa: self.conv_check.primal_noise_floor_kappa,
1579 acceptable_progress_kappa: self.conv_check.acceptable_progress_kappa,
1580 acceptable_window: std::collections::VecDeque::new(),
1581 acceptable_progress_refusals: 0,
1582 dual_inf_scale_kappa: self.conv_check.dual_inf_scale_kappa,
1583 dual_floor_reported: false,
1584 veto_fired: false,
1585 acceptable_veto_fired: false,
1586 masked_acceptable_veto_fired: false,
1587 veto_extra_iters: 0,
1588 rel_infeas_extra_iters: 0,
1589 prev_rel_viol: f64::NAN,
1590 });
1591
1592 let init: Box<dyn crate::init::r#trait::IterateInitializer> = if self.warm_start_init_point
1593 {
1594 Box::new(WarmStartIterateInitializer::with_options(
1595 resolved_warm_options(&self.warm, &self.init),
1596 ))
1597 } else {
1598 let mut d = DefaultIterateInitializer::with_eq_mult_calculator(Box::new(
1599 LeastSquareMults::new(),
1600 ));
1601 d.bound_push = self.init.bound_push;
1602 d.bound_frac = self.init.bound_frac;
1603 d.slack_bound_push = self.init.slack_bound_push;
1604 d.slack_bound_frac = self.init.slack_bound_frac;
1605 d.constr_mult_init_max = self.init.constr_mult_init_max;
1606 d.bound_mult_init_val = self.init.bound_mult_init_val;
1607 d.bound_mult_init_method = self.init.bound_mult_init_method.clone();
1608 d.least_square_init_primal = self.init.least_square_init_primal;
1609 Box::new(d)
1610 };
1611
1612 let eq_mult: Box<dyn crate::eq_mult::r#trait::EqMultCalculator> =
1613 Box::new(LeastSquareMults::new());
1614
1615 let hess: Box<dyn crate::hess::r#trait::HessianUpdater> = match self.hessian_approximation {
1616 HessianApproxChoice::Exact => Box::new(ExactHessianUpdater::new()),
1617 HessianApproxChoice::LimitedMemory => Box::new(LimMemQuasiNewtonUpdater {
1618 update_type: self.limited_memory_update_type,
1619 max_history: self.limited_memory_max_history,
1620 init_val_max: self.limited_memory_init_val_max,
1621 init_val_min: self.limited_memory_init_val_min,
1622 initial_approx: self.limited_memory_initialization,
1623 init_val: self.limited_memory_init_val,
1624 max_skipping: self.limited_memory_max_skipping,
1625 nonlinear_vars: self.limited_memory_nonlinear_vars.clone(),
1626 ..LimMemQuasiNewtonUpdater::default()
1627 }),
1628 HessianApproxChoice::Partitioned => {
1629 let mut u =
1630 crate::hess::partitioned_quasi_newton::PartitionedQuasiNewtonUpdater::new(
1631 self.partitioned_update_type,
1632 );
1633 u.max_element = self.partitioned_max_element;
1634 u.objective_vars = self.objective_nonlinear_vars.clone();
1635 u.curvature_cap = self.partitioned_curvature_cap;
1636 u.mode = self.partitioned_elements;
1637 u.block_size = self.partitioned_block_size;
1638 // Damped BFGS is the right pairing for a Lagrangian
1639 // block: unlike a single constraint's Hessian, the
1640 // Lagrangian's is the object the IPM wants a positive
1641 // definite model of, and the sign problem that makes
1642 // damping wrong per constraint does not arise. Only when
1643 // the caller has not named a formula.
1644 if self.partitioned_elements
1645 == crate::hess::partitioned_quasi_newton::ElementMode::PrimalBlock
1646 && !self.partitioned_update_type_was_set
1647 {
1648 u.update_type = UpdateType::Bfgs;
1649 }
1650 u.init_val_min = self.limited_memory_init_val_min;
1651 u.init_val_max = self.limited_memory_init_val_max;
1652 Box::new(u)
1653 }
1654 HessianApproxChoice::FiniteDifference => {
1655 let mut u = crate::hess::fd_hessian::FdHessianUpdater::new(self.fd_hessian_pattern);
1656 u.coloring = self.fd_hessian_coloring;
1657 u.reuse_tol = self.fd_hessian_reuse_tol;
1658 u.objective_vars = self.objective_nonlinear_vars.clone();
1659 u.nonlinear_vars = self.limited_memory_nonlinear_vars.clone();
1660 Box::new(u)
1661 }
1662 };
1663
1664 let iter_output: Box<dyn crate::output::r#trait::IterationOutput> = {
1665 use crate::output::orig::{InfPrTag, PrintInfoString};
1666 let mut o = OrigIterationOutput::new();
1667 o.print_frequency_iter = self.output.print_frequency_iter;
1668 o.print_frequency_time = self.output.print_frequency_time;
1669 o.print_info_string = if self.output.print_info_string {
1670 PrintInfoString::Yes
1671 } else {
1672 PrintInfoString::No
1673 };
1674 o.inf_pr_output = if self.output.inf_pr_output_internal {
1675 InfPrTag::Internal
1676 } else {
1677 InfPrTag::Original
1678 };
1679 Box::new(o)
1680 };
1681
1682 AlgorithmBundle {
1683 mu_update,
1684 conv_check,
1685 init,
1686 eq_mult,
1687 hess,
1688 line_search,
1689 iter_output,
1690 search_dir,
1691 }
1692 }
1693}
1694
1695#[cfg(test)]
1696mod tests {
1697 use super::*;
1698
1699 #[test]
1700 fn warm_options_take_the_init_default_not_their_own() {
1701 let mut init = InitOptions::default();
1702 init.bound_mult_init_val = 10.0; // the Mehrotra override value
1703 let mut warm = WarmStartOptions::default();
1704 warm.bound_mult_init_val = 123.0; // stale copy must lose
1705 let resolved = resolved_warm_options(&warm, &init);
1706 assert_eq!(resolved.bound_mult_init_val, 10.0);
1707 // everything else passes through untouched
1708 assert_eq!(resolved.mult_bound_push, warm.mult_bound_push);
1709 assert_eq!(resolved.target_mu, warm.target_mu);
1710 }
1711
1712 #[test]
1713 fn default_builder_assembles() {
1714 let bundle = AlgorithmBuilder::new().build();
1715 // Sanity: the placeholder traits compile and the boxed
1716 // strategies don't panic on construction.
1717 let _ = bundle.line_search.acceptor();
1718 assert!(bundle.search_dir.is_none());
1719 }
1720
1721 #[test]
1722 fn build_with_backend_assembles_search_dir_chain() {
1723 // Drive the builder with the FERAL backend factory; the
1724 // resulting bundle should expose a populated `PdSearchDirCalc`.
1725 let factory: LinearBackendFactory = Box::new(|_| {
1726 Box::new(pounce_feral::FeralSolverInterface::new())
1727 as Box<dyn SparseSymLinearSolverInterface>
1728 });
1729 let bundle = AlgorithmBuilder::new().build_with_backend(factory);
1730 assert!(bundle.search_dir.is_some());
1731 }
1732
1733 #[test]
1734 fn limited_memory_sr1_propagates() {
1735 let b = AlgorithmBuilder {
1736 hessian_approximation: HessianApproxChoice::LimitedMemory,
1737 limited_memory_update_type: UpdateType::Sr1,
1738 ..AlgorithmBuilder::default()
1739 };
1740 let _bundle = b.build();
1741 }
1742
1743 #[test]
1744 fn every_strategy_combination_assembles_without_panic() {
1745 let solvers = [LinearSolverChoice::Ma57, LinearSolverChoice::Feral];
1746 let mu = [MuStrategyChoice::Monotone, MuStrategyChoice::Adaptive];
1747 let hess = [
1748 HessianApproxChoice::Exact,
1749 HessianApproxChoice::LimitedMemory,
1750 ];
1751 let ls = [
1752 LineSearchChoice::Filter,
1753 LineSearchChoice::CgPenalty,
1754 LineSearchChoice::Penalty,
1755 ];
1756 for &linear_solver in &solvers {
1757 for &mu_strategy in &mu {
1758 for &hessian_approximation in &hess {
1759 for &line_search_method in &ls {
1760 let _ = AlgorithmBuilder {
1761 algorithm: AlgorithmChoice::default(),
1762 linear_solver,
1763 linear_system_scaling: LinearSystemScalingChoice::None,
1764 linear_scaling_on_demand: true,
1765 mu_strategy,
1766 mu_oracle: MuOracleKind::QualityFunction,
1767 hessian_approximation,
1768 partitioned_update_type: UpdateType::Sr1,
1769 partitioned_update_type_was_set: false,
1770 partitioned_max_element: 64,
1771 objective_nonlinear_vars: None,
1772 partitioned_curvature_cap: Number::INFINITY,
1773 partitioned_elements:
1774 crate::hess::partitioned_quasi_newton::ElementMode::PerConstraint,
1775 partitioned_block_size: 64,
1776 fd_hessian_pattern: crate::hess::fd_hessian::FdPatternSource::Declared,
1777 fd_hessian_coloring: crate::hess::fd_hessian::FdColoring::Cpr,
1778 fd_hessian_reuse_tol: 0.0,
1779 limited_memory_update_type: UpdateType::Bfgs,
1780 limited_memory_max_history: 6,
1781 limited_memory_init_val_max: 1e8,
1782 limited_memory_init_val_min: 1e-8,
1783 limited_memory_initialization: InitialApprox::Scalar1,
1784 limited_memory_init_val: 1.0,
1785 limited_memory_max_skipping: 2,
1786 limited_memory_nonlinear_vars: None,
1787 line_search_method,
1788 warm_start_init_point: false,
1789 mehrotra_algorithm: false,
1790 fast_step_computation: false,
1791 kappa_sigma: 1e10,
1792 recalc_y: false,
1793 recalc_y_feas_tol: 1e-6,
1794 kappa_d: 1e-5,
1795 s_max: 100.0,
1796 tiny_step_tol: 10.0 * Number::EPSILON,
1797 tiny_step_y_tol: 1e-2,
1798 diverging_iterates_tol: 1e20,
1799 dual_diverging_streak: 0,
1800 dual_divergence_retry_step_tol: 1e-5,
1801 dual_divergence_retry_du_floor: 1e2,
1802 resto_decline_deferrals: 1,
1803 resto_decline_progress_ratio: 0.5,
1804 neg_curv_escapes: 1,
1805 limited_memory_ls_failure_restarts: 0,
1806 kkt_fidelity_tol: 0.0,
1807 conv_check: ConvCheckOptions::default(),
1808 mu: MuOptions::default(),
1809 line_search: LineSearchOptions::default(),
1810 refinement: RefinementOptions::default(),
1811 perturbation: PerturbationOptions::default(),
1812 resto: RestoOptions::default(),
1813 output: OutputOptions::default(),
1814 warm: WarmStartOptions::default(),
1815 sqp: crate::sqp::SqpOptions::default(),
1816 sqp_qp: pounce_qp::QpOptions::sqp_subproblem(),
1817 init: InitOptions::default(),
1818 kkt_schur: None,
1819 quality_escalation_counter: None,
1820 }
1821 .build();
1822 }
1823 }
1824 }
1825 }
1826 }
1827}