Skip to main content

pounce_cli/
verify.rs

1//! `pounce verify <problem.nl> <claim.sol>` — independent solution checker.
2//!
3//! # Why this exists
4//!
5//! When pounce is a *tool an agent calls*, the agent should never be the
6//! thing you trust for "the solution satisfies the constraints." Trust
7//! belongs to a small, deterministic checker that re-derives the answer
8//! from the **canonical** problem — not from the agent's narration and not
9//! even from the solver's own exit string. Optimization is the rare setting
10//! where this is cheap: a claimed `x*` is just numbers, and feasibility is
11//! one constraint evaluation (`g_l ≤ g(x*) ≤ g_u`, `x_l ≤ x* ≤ x_u`),
12//! `O(nnz)` work with no resolve.
13//!
14//! `pounce verify` loads the canonical `.nl`, reads a claimed `.sol`, and
15//! reports the worst constraint/bound violation (and, when the `.sol`
16//! carries constraint duals, a first-order/KKT stationarity residual). It
17//! defends the three agent-workflow failure modes:
18//!
19//! * **fabrication** ("here's a solution that looks like pounce ran") —
20//!   invented numbers fail the residual check against the real model;
21//! * **ignoring the solver** — a downstream consumer gates on the receipt's
22//!   `verified: true` plus the problem hash, not on prose;
23//! * **solving the wrong problem** (dropping/relaxing a constraint to dodge
24//!   infeasibility) — the check runs against the *canonical* constraints
25//!   and bounds, so a point that is only feasible for a relaxed model is
26//!   caught here.
27//!
28//! The JSON receipt content-addresses both inputs by SHA-256 so a consumer
29//! can confirm *which* problem was verified. When the `POUNCE_VERIFY_KEY`
30//! environment variable holds a secret the agent does not have, the receipt
31//! is additionally signed with HMAC-SHA256 over a float-free preimage (see
32//! [`signing_preimage`]) — so an agent cannot mint a receipt that a consumer
33//! holding the key will accept. The consumer recomputes the HMAC over the
34//! same preimage and compares.
35//!
36//! Verdict / exit code: `0` when every violation is within tolerance
37//! (`FEASIBLE`); `20` when a violation exceeds tolerance (`INFEASIBLE`);
38//! `2` on a usage or I/O error. Optimality is reported but, by default,
39//! does not gate — feasibility is the rigorous, sign-convention-independent
40//! guarantee; pass `--require-optimal` to also gate on the stationarity
41//! residual.
42//!
43//! # Two different complementarity quantities (gh #516)
44//!
45//! "Complementarity" names two distinct residuals, and printing either one
46//! under the bare label invites a comparison against the other:
47//!
48//! * **constraint** complementarity — `max_i |λ_i| · dist(g_i, nearest
49//!   finite side)` over **rows**, from the `.sol`'s constraint duals. This
50//!   is the one `verify` has always computed.
51//! * **bound** complementarity — `max_j max(|z_L·(x−x_L)|, |z_U·(x_U−x)|)`
52//!   over **variables**, from the bound multipliers. This is what Ipopt
53//!   prints as `Complementarity`, and it needs the `ipopt_zL_out` /
54//!   `ipopt_zU_out` `.sol` suffixes.
55//!
56//! They can differ by many orders of magnitude at the same point and
57//! neither is wrong for what it measures. `verify` therefore names both
58//! explicitly, reads the bound multipliers when the `.sol` carries them,
59//! and says "not checked" — rather than nothing — when it does not.
60//!
61//! Those same suffixes also sharpen stationarity: without them the residual
62//! is bound-*projected* and cannot see a bound multiplier that is missing or
63//! wrong (gh #495); with them the exact residual is available, and
64//! `--require-optimal` gates on it.
65
66use crate::nl_reader::{self, DiscreteCensus};
67use pounce_common::tolerance::is_negligible;
68use pounce_common::types::{Number, lower_bound_present, upper_bound_present};
69use pounce_nlp::tnlp::{BoundsInfo, IndexStyle, SparsityRequest, TNLP};
70use std::path::PathBuf;
71use std::process::ExitCode;
72
73/// Parsed `verify` subcommand arguments.
74#[derive(Debug, Clone)]
75pub struct VerifyArgs {
76    pub nl: PathBuf,
77    pub sol: PathBuf,
78    /// Max `|violation|` of any constraint or bound still called feasible.
79    ///
80    /// Applied **per row, relative to the row's magnitude** — the accepting
81    /// test is `|viol| <= feas_tol * max(1, |row|)` (see [`row_is_violated`]),
82    /// so on a row of size `2e6` a violation of `1.0` passes at `1e-6`. Use
83    /// [`abs_feas_tol`](VerifyArgs::abs_feas_tol) for a plain absolute test.
84    pub feas_tol: Number,
85    /// `--abs-feas-tol T` — judge every row by the plain absolute test
86    /// `|viol| <= T`, replacing the relative test above (gh#987).
87    pub abs_feas_tol: Option<Number>,
88    /// Max stationarity residual still called first-order optimal.
89    pub opt_tol: Number,
90    /// `--json-output PATH` — write the machine-readable receipt to PATH.
91    pub json_output: Option<PathBuf>,
92    /// `--require-optimal` — also gate the exit code on the stationarity
93    /// residual (needs duals in the `.sol`).
94    pub require_optimal: bool,
95}
96
97impl Default for VerifyArgs {
98    fn default() -> Self {
99        VerifyArgs {
100            nl: PathBuf::new(),
101            sol: PathBuf::new(),
102            feas_tol: 1e-6,
103            abs_feas_tol: None,
104            opt_tol: 1e-6,
105            json_output: None,
106            require_optimal: false,
107        }
108    }
109}
110
111const USAGE: &str = "\
112Usage: pounce verify <problem.nl> <claim.sol> [OPTIONS]
113
114Independently check that the solution in <claim.sol> satisfies the
115constraints and bounds of the canonical problem <problem.nl>. Re-derives
116feasibility from the model itself — it does not trust the .sol's status
117line or rerun the solver.
118
119Arguments:
120  <problem.nl>            canonical AMPL .nl problem (the source of truth)
121  <claim.sol>            claimed AMPL .sol solution to check
122
123Options:
124  --feas-tol <t>         feasibility tolerance (default 1e-6). RELATIVE per
125                         row: a row passes when |violation| <= t*max(1,|row|),
126                         so a violation of 1.0 on a row of size 2e6 passes at
127                         1e-6. Use --abs-feas-tol for an absolute test.
128  --abs-feas-tol <t>     judge every row by the absolute test |violation| <= t
129                         (replaces the relative --feas-tol test)
130  --opt-tol <t>          stationarity tolerance (default 1e-6)
131  --require-optimal      also fail if the KKT stationarity residual
132                         exceeds --opt-tol (needs duals in the .sol)
133  --json-output <path>   write a JSON verification receipt to <path>
134  -h, --help             print this message
135
136Complementarity: two different residuals carry that name, and they can
137differ by many orders of magnitude at the same point.
138  * constraint complementarity (rows, |lambda|*slack) is computed from the
139    .sol's constraint duals and is always reported alongside stationarity.
140  * bound complementarity (vars, |z|*slack) is the quantity Ipopt prints as
141    `Complementarity`. It needs the bound multipliers, which reach a .sol
142    only as the `ipopt_zL_out` / `ipopt_zU_out` suffixes; without them it is
143    reported as `not checked`, never as a number.
144Do not compare the row quantity against a solver's `Complementarity` line.
145
146Exit code: 0 = verified feasible, 20 = violation exceeds tolerance,
1472 = usage/IO error.";
148
149/// Entry point dispatched from `main` when argv[1] == "verify".
150pub fn run_from_argv(rest: &[String]) -> ExitCode {
151    let args = match parse_verify_argv(rest) {
152        Ok(Some(a)) => a,
153        Ok(None) => {
154            // help was requested
155            println!("{USAGE}");
156            return ExitCode::SUCCESS;
157        }
158        Err(msg) => {
159            eprintln!("pounce verify: {msg}");
160            eprintln!("{USAGE}");
161            return ExitCode::from(2);
162        }
163    };
164    run(&args)
165}
166
167fn parse_verify_argv(rest: &[String]) -> Result<Option<VerifyArgs>, String> {
168    let mut a = VerifyArgs::default();
169    let mut positionals: Vec<PathBuf> = Vec::new();
170    let mut it = rest.iter();
171    while let Some(arg) = it.next() {
172        match arg.as_str() {
173            "-h" | "--help" => return Ok(None),
174            "--feas-tol" => {
175                let v = it.next().ok_or("--feas-tol requires a value")?;
176                a.feas_tol = v.parse().map_err(|e| format!("--feas-tol: {e}"))?;
177            }
178            "--abs-feas-tol" => {
179                let v = it.next().ok_or("--abs-feas-tol requires a value")?;
180                a.abs_feas_tol = Some(v.parse().map_err(|e| format!("--abs-feas-tol: {e}"))?);
181            }
182            "--opt-tol" => {
183                let v = it.next().ok_or("--opt-tol requires a value")?;
184                a.opt_tol = v.parse().map_err(|e| format!("--opt-tol: {e}"))?;
185            }
186            "--require-optimal" => a.require_optimal = true,
187            "--json-output" => {
188                let v = it.next().ok_or("--json-output requires a value")?;
189                a.json_output = Some(PathBuf::from(v));
190            }
191            other if other.starts_with('-') => {
192                return Err(format!("unknown flag `{other}`"));
193            }
194            _ => positionals.push(PathBuf::from(arg)),
195        }
196    }
197    match positionals.len() {
198        0 | 1 => Err("expected two positional arguments: <problem.nl> <claim.sol>".to_string()),
199        2 => {
200            a.nl = positionals[0].clone();
201            a.sol = positionals[1].clone();
202            Ok(Some(a))
203        }
204        n => Err(format!("expected 2 positional arguments, got {n}")),
205    }
206}
207
208/// The fully-evaluated verification result. Serialized to the JSON
209/// receipt and rendered to the console.
210#[derive(Debug)]
211pub struct VerifyOutcome {
212    pub n_vars: usize,
213    pub n_cons: usize,
214    pub nl_sha256: String,
215    pub sol_sha256: String,
216    pub solve_result_num: Option<i32>,
217    pub feas_tol: Number,
218    pub opt_tol: Number,
219    // feasibility
220    pub max_con_violation: Number,
221    pub worst_con: Option<RowReport>,
222    pub max_bound_violation: Number,
223    pub worst_bound: Option<RowReport>,
224    pub feasible: bool,
225    // optimality (only when duals supplied)
226    pub objective: Option<Number>,
227    pub duals_present: bool,
228    pub stationarity: Option<Number>,
229    pub dual_sign: Option<i32>,
230    /// `max_i |λ_i| · dist(g_i, active side)` over **rows**. NOT the
231    /// quantity a solver reports as `Complementarity` — see
232    /// [`bound_complementarity`](VerifyOutcome::bound_complementarity).
233    pub constraint_complementarity: Option<Number>,
234    /// Whether the `.sol` carried `ipopt_zL_out` / `ipopt_zU_out`.
235    pub bound_multipliers_present: bool,
236    /// `max_j max(|z_L·(x−x_L)|, |z_U·(x_U−x)|)` over **variables** — the
237    /// quantity Ipopt prints as `Complementarity`. `None` when the `.sol`
238    /// carried no bound multipliers, in which case it is *not checked*
239    /// rather than zero.
240    pub bound_complementarity: Option<Number>,
241    /// Exact (non-projected) dual infeasibility
242    /// `‖∇f + sign·Jᵀλ − (z_L^suffix + z_U^suffix)‖∞`, available only when
243    /// both duals and bound multipliers are present.
244    pub stationarity_with_bound_multipliers: Option<Number>,
245    pub optimal: Option<bool>,
246    /// The `.nl` header's declared binary / integer census (gh#987).
247    pub discrete: Option<DiscreteCensus>,
248    /// Max distance to the nearest integer over the declared integer
249    /// columns, and the worst column. `None` when integrality was not
250    /// checked (no discrete variables, or their columns cannot be told apart
251    /// from the header alone).
252    pub integrality: Option<(Number, usize)>,
253    /// The tolerance used for the absolute test, when `--abs-feas-tol` was set.
254    pub abs_feas_tol: Option<Number>,
255    // final
256    pub verified: bool,
257}
258
259#[derive(Debug, Clone)]
260pub struct RowReport {
261    pub index: usize,
262    pub name: String,
263    pub value: Number,
264    pub lo: Number,
265    pub hi: Number,
266    pub violation: Number,
267}
268
269// `is_finite_bound(b) = b > NLP_LOWER_BOUND_INF && b < NLP_UPPER_BOUND_INF`
270// used to live here — a *band* membership test applied to lower and upper
271// bounds alike (gh #403). A real upper bound of `-5e20` failed it, so
272// `box_violation` scored `0.0` against it and `verify` reported ACCEPTED for a
273// `.sol` that violates a declared bound. Presence is directional; use
274// `lower_bound_present` / `upper_bound_present` from `pounce_common::types`,
275// picking the one that matches the side you hold.
276
277/// `g_l ≤ v ≤ g_u` violation: how far `v` is outside the box, 0 if inside.
278///
279/// A non-finite `v` (NaN or ±∞) is treated as an infinite violation, never
280/// as feasible: `NaN`-laden arithmetic would otherwise collapse to `0.0`
281/// through `f64::max` (which drops NaN operands) and let a fabricated `.sol`
282/// slip past the feasibility gate — the exact threat this checker defends
283/// against. An unbounded variable pinned at ±∞ is likewise not a real point.
284/// The natural magnitude of a row, for a scale-relative feasibility test.
285///
286/// `verify` reads a `.nl` and a `.sol`; no solver scaling has been applied, so
287/// the magnitude has to come from the row's own numbers — the evaluated value
288/// and whichever bounds are finite. Infinite bounds carry no magnitude
289/// information and are skipped.
290pub(crate) fn row_magnitude(value: Number, lo: Number, hi: Number) -> Number {
291    let mut m = if value.is_finite() { value.abs() } else { 0.0 };
292    if lower_bound_present(lo) {
293        m = m.max(lo.abs());
294    }
295    if upper_bound_present(hi) {
296        m = m.max(hi.abs());
297    }
298    m
299}
300
301/// Whether a row's violation is real, judged relative to the row's own
302/// magnitude.
303///
304/// An absolute tolerance is meaningless against a row evaluating near `1e13`:
305/// `--feas-tol 1e-6` is unreachable there, so a solution correct to eleven
306/// relative digits was reported REJECTED. Scaling the tolerance by the row
307/// magnitude makes the verdict independent of how the model happens to be
308/// written.
309///
310/// Uses the **accepting** direction (`is_negligible`), which is never stricter
311/// than the plain absolute `tol`. A pure relative test was tried first and
312/// rejected genuine solutions: the solver converges to *absolute* residuals, so
313/// on a row of magnitude `1e-3` a residual of `1e-8` is converged, while a
314/// relative test at `tol = 1e-6` would demand `1e-9`.
315///
316/// The non-finite case is handled here rather than inside the primitive, which
317/// reports an unjudgeable value as not-negligible-and-not-significant. A `.sol`
318/// carrying `NaN` or `±inf` is not a point at all and must be rejected — which
319/// is what `box_violation` returning infinity encodes.
320pub(crate) fn row_is_violated(viol: Number, magnitude: Number, feas_tol: Number) -> bool {
321    if !viol.is_finite() {
322        return true;
323    }
324    !is_negligible(viol, magnitude, feas_tol)
325}
326
327/// Row verdict under the active tolerance mode: the relative accepting test by
328/// default, the plain absolute test when `--abs-feas-tol` was given.
329fn judge_row(args: &VerifyArgs, viol: Number, magnitude: Number) -> bool {
330    match args.abs_feas_tol {
331        Some(t) => !viol.is_finite() || viol > t,
332        None => row_is_violated(viol, magnitude, args.feas_tol),
333    }
334}
335
336/// The one-line warning for a `.nl` that declares integer variables, shared by
337/// the solve path and `verify` (gh#987). `None` when there are none.
338pub fn integrality_notice(c: Option<DiscreteCensus>) -> Option<String> {
339    let c = c.filter(|c| c.total() > 0)?;
340    Some(format!(
341        "the .nl declares {} discrete variable(s) ({} binary and {} integer appearing \
342         linearly, {} appearing nonlinearly); pounce is a continuous solver and solved \
343         the RELAXATION, so the result is a bound on the MIP optimum, not an \
344         integer-feasible solution",
345        c.total(),
346        c.binary,
347        c.integer,
348        c.nonlinear_discrete()
349    ))
350}
351
352pub(crate) fn box_violation(v: Number, lo: Number, hi: Number) -> Number {
353    if !v.is_finite() {
354        return Number::INFINITY;
355    }
356    let below = if lower_bound_present(lo) {
357        lo - v
358    } else {
359        Number::NEG_INFINITY
360    };
361    let above = if upper_bound_present(hi) {
362        v - hi
363    } else {
364        Number::NEG_INFINITY
365    };
366    below.max(above).max(0.0)
367}
368
369pub fn run(args: &VerifyArgs) -> ExitCode {
370    let outcome = match evaluate(args) {
371        Ok(o) => o,
372        Err(msg) => {
373            eprintln!("pounce verify: {msg}");
374            return ExitCode::from(2);
375        }
376    };
377    print_report(args, &outcome);
378
379    if let Some(path) = &args.json_output {
380        let json = receipt_json(args, &outcome);
381        if let Err(e) = std::fs::write(path, json.as_bytes()) {
382            eprintln!(
383                "pounce verify: failed to write receipt {}: {e}",
384                path.display()
385            );
386            return ExitCode::from(2);
387        }
388        let signed = std::env::var(KEY_ENV)
389            .map(|k| !k.is_empty())
390            .unwrap_or(false);
391        println!(
392            "  receipt: {}{}",
393            path.display(),
394            if signed {
395                "  (signed: HMAC-SHA256)"
396            } else {
397                ""
398            }
399        );
400    }
401
402    if outcome.verified {
403        ExitCode::SUCCESS
404    } else {
405        ExitCode::from(20)
406    }
407}
408
409fn evaluate(args: &VerifyArgs) -> Result<VerifyOutcome, String> {
410    // --- read + hash the two inputs (content-address the receipt) ---
411    let nl_bytes =
412        std::fs::read(&args.nl).map_err(|e| format!("cannot read {}: {e}", args.nl.display()))?;
413    let sol_bytes =
414        std::fs::read(&args.sol).map_err(|e| format!("cannot read {}: {e}", args.sol.display()))?;
415    let nl_sha256 = sha256::hex(&nl_bytes);
416    let sol_sha256 = sha256::hex(&sol_bytes);
417
418    // --- canonical problem ---
419    let prob = nl_reader::read_nl_file(&args.nl)?;
420    let n = prob.n;
421    let m = prob.m;
422    let con_names = prob.con_names.clone();
423    let var_names = prob.var_names.clone();
424    let prob_discrete = prob.n_discrete;
425    let prob_nl_counts = prob.nl_counts;
426    // `NlTnlp` evaluates a maximize objective negated, while the `.sol`
427    // bound blocks are `sigma*z_L` / `-sigma*z_U` (gh#959); the exact
428    // dual-infeasibility check has to bring them into the same frame.
429    let obj_sigma: Number = if prob.minimize { 1.0 } else { -1.0 };
430    let mut tnlp = nl_reader::NlTnlp::new(prob);
431
432    let info = tnlp
433        .get_nlp_info()
434        .ok_or("get_nlp_info failed on the .nl")?;
435    let nnz = info.nnz_jac_g.max(0) as usize;
436    let fortran = matches!(info.index_style, IndexStyle::Fortran);
437
438    // --- claimed solution ---
439    let sol_text = String::from_utf8_lossy(&sol_bytes);
440    let parsed = parse_sol(&sol_text)?;
441    if parsed.x.len() != n {
442        return Err(format!(
443            "solution has {} primal values but the problem has {n} variables \
444             (is this the right .sol for this .nl?)",
445            parsed.x.len()
446        ));
447    }
448    let x = parsed.x;
449    let duals_present = !parsed.lambda.is_empty();
450    if duals_present && parsed.lambda.len() != m {
451        return Err(format!(
452            "solution carries {} dual values but the problem has {m} constraints",
453            parsed.lambda.len()
454        ));
455    }
456
457    // --- bounds ---
458    let mut x_l = vec![0.0; n];
459    let mut x_u = vec![0.0; n];
460    let mut g_l = vec![0.0; m];
461    let mut g_u = vec![0.0; m];
462    if !tnlp.get_bounds_info(BoundsInfo {
463        x_l: &mut x_l,
464        x_u: &mut x_u,
465        g_l: &mut g_l,
466        g_u: &mut g_u,
467    }) {
468        return Err("get_bounds_info failed".to_string());
469    }
470
471    // --- bound feasibility ---
472    let mut max_bound_violation = 0.0_f64;
473    let mut worst_bound: Option<RowReport> = None;
474    let mut any_bound_violated = false;
475    for j in 0..n {
476        let viol = box_violation(x[j], x_l[j], x_u[j]);
477        if judge_row(args, viol, row_magnitude(x[j], x_l[j], x_u[j])) {
478            any_bound_violated = true;
479        }
480        if viol > max_bound_violation {
481            max_bound_violation = viol;
482            worst_bound = Some(RowReport {
483                index: j,
484                name: name_at(&var_names, j, 'x'),
485                value: x[j],
486                lo: x_l[j],
487                hi: x_u[j],
488                violation: viol,
489            });
490        }
491    }
492
493    // --- constraint feasibility ---
494    let mut g = vec![0.0; m];
495    if !tnlp.eval_g(&x, true, &mut g) {
496        return Err("eval_g failed at the claimed solution".to_string());
497    }
498    let mut max_con_violation = 0.0_f64;
499    let mut worst_con: Option<RowReport> = None;
500    let mut any_con_violated = false;
501    for i in 0..m {
502        let viol = box_violation(g[i], g_l[i], g_u[i]);
503        if judge_row(args, viol, row_magnitude(g[i], g_l[i], g_u[i])) {
504            any_con_violated = true;
505        }
506        if viol > max_con_violation {
507            max_con_violation = viol;
508            worst_con = Some(RowReport {
509                index: i,
510                name: name_at(&con_names, i, 'c'),
511                value: g[i],
512                lo: g_l[i],
513                hi: g_u[i],
514                violation: viol,
515            });
516        }
517    }
518
519    // Per-row and scale-relative: a single absolute threshold across rows of
520    // wildly different magnitude answers a different question for each of them.
521    // Integrality (gh#987). pounce solves the continuous relaxation, so a
522    // `.sol` for a MIP `.nl` can be feasible for every row and still be
523    // fractional. Header line 7 says how many columns are discrete — `nbv`
524    // and `niv` for those appearing only linearly, `nlvbi nlvci nlvoi` for
525    // those appearing nonlinearly — and Gay's variable ordering places them
526    // (see [`DiscreteCensus::integer_columns`]). When the columns cannot be
527    // identified the check is skipped and the report says NOT CHECKED rather
528    // than implying it passed.
529    let discrete = prob_discrete;
530    let integrality = discrete
531        .filter(|d| d.total() > 0)
532        .and_then(|d| d.integer_columns(n, prob_nl_counts))
533        .filter(|cols| !cols.is_empty())
534        .map(|cols| {
535            let mut worst = (0.0_f64, cols[0]);
536            for j in cols {
537                let xj = x[j];
538                let dist = if xj.is_finite() {
539                    (xj - xj.round()).abs()
540                } else {
541                    Number::INFINITY
542                };
543                if dist > worst.0 {
544                    worst = (dist, j);
545                }
546            }
547            worst
548        });
549    let integrality_ok =
550        integrality.is_none_or(|(d, _)| d <= args.abs_feas_tol.unwrap_or(args.feas_tol));
551    let feasible = !any_con_violated && !any_bound_violated && integrality_ok;
552
553    // --- objective ---
554    let objective = tnlp.eval_f(&x, true);
555
556    // --- bound multipliers, when the `.sol` exported them (gh #516) ---
557    //
558    // The bound complementarity Ipopt prints as `Complementarity` cannot be
559    // computed from the primal and the constraint duals alone; it needs
560    // `z_L` / `z_U`, which reach a `.sol` only as the `ipopt_zL_out` /
561    // `ipopt_zU_out` variable suffixes. Absent them the quantity is *not
562    // checked* — never silently reported as the row quantity.
563    let bound_multipliers_present = parsed.z_l.is_some() || parsed.z_u.is_some();
564    let z_l_suf = parsed.z_l.clone().unwrap_or_else(|| vec![0.0; n]);
565    let z_u_suf = parsed.z_u.clone().unwrap_or_else(|| vec![0.0; n]);
566    let bound_complementarity = if bound_multipliers_present {
567        Some(bound_complementarity(&x, &x_l, &x_u, &z_l_suf, &z_u_suf))
568    } else {
569        None
570    };
571
572    // --- first-order / KKT stationarity (only when duals are supplied) ---
573    let mut stationarity = None;
574    let mut dual_sign = None;
575    let mut constraint_complementarity = None;
576    let mut stationarity_with_bound_multipliers = None;
577    let mut optimal = None;
578    // A problem with no rows has no constraint duals to carry, so `∇f` alone
579    // is the Lagrangian gradient and the residual is available from an empty
580    // dual block — which is what a `.sol` for a bounds-only model has.
581    if duals_present || m == 0 {
582        let lambda = &parsed.lambda;
583
584        // ∇f(x*)
585        let mut grad_f = vec![0.0; n];
586        tnlp.eval_grad_f(&x, true, &mut grad_f);
587
588        // Jacobian triplets (structure then values).
589        let mut irow = vec![0i32; nnz];
590        let mut jcol = vec![0i32; nnz];
591        tnlp.eval_jac_g(
592            Some(&x),
593            true,
594            SparsityRequest::Structure {
595                irow: &mut irow,
596                jcol: &mut jcol,
597            },
598        );
599        let mut jval = vec![0.0; nnz];
600        tnlp.eval_jac_g(
601            Some(&x),
602            true,
603            SparsityRequest::Values { values: &mut jval },
604        );
605
606        // AMPL's dual sign convention can flip relative to ours; rather
607        // than guess, compute the bound-projected stationarity residual
608        // for both signs and keep the better one. A genuine KKT point is
609        // stationary for exactly one of them; we report which.
610        let s_pos = lagrangian_gradient(1.0, &grad_f, &irow, &jcol, &jval, fortran, lambda);
611        let s_neg = lagrangian_gradient(-1.0, &grad_f, &irow, &jcol, &jval, fortran, lambda);
612        let resid_pos = bound_projected_residual(&s_pos, &x, &x_l, &x_u, args.opt_tol);
613        let resid_neg = bound_projected_residual(&s_neg, &x, &x_l, &x_u, args.opt_tol);
614        let (best_resid, sign, s) = if resid_pos <= resid_neg {
615            (resid_pos, 1, &s_pos)
616        } else {
617            (resid_neg, -1, &s_neg)
618        };
619        stationarity = Some(best_resid);
620        dual_sign = Some(sign);
621        constraint_complementarity = Some(row_complementarity(lambda, &g, &g_l, &g_u));
622
623        // With the bound multipliers in hand the residual no longer has to
624        // be projected: the exact dual infeasibility is available, and it
625        // is what a solver reports. It is also the strictly sharper check —
626        // the projection can only *remove* residual — so `--require-optimal`
627        // gates on it whenever it exists.
628        if bound_multipliers_present {
629            stationarity_with_bound_multipliers =
630                Some(exact_dual_infeasibility(s, &z_l_suf, &z_u_suf, obj_sigma));
631        }
632        let gate = stationarity_with_bound_multipliers.unwrap_or(best_resid);
633        optimal = Some(gate <= args.opt_tol);
634    }
635
636    // Verified = feasible (always required) AND, if --require-optimal,
637    // also first-order optimal.
638    let verified = feasible && (!args.require_optimal || optimal.unwrap_or(false));
639
640    Ok(VerifyOutcome {
641        n_vars: n,
642        n_cons: m,
643        nl_sha256,
644        sol_sha256,
645        solve_result_num: parsed.solve_result_num,
646        feas_tol: args.feas_tol,
647        opt_tol: args.opt_tol,
648        max_con_violation,
649        worst_con,
650        max_bound_violation,
651        worst_bound,
652        feasible,
653        objective,
654        duals_present,
655        stationarity,
656        dual_sign,
657        constraint_complementarity,
658        bound_multipliers_present,
659        bound_complementarity,
660        stationarity_with_bound_multipliers,
661        optimal,
662        discrete,
663        integrality,
664        abs_feas_tol: args.abs_feas_tol,
665        verified,
666    })
667}
668
669/// `s = ∇f + sign·Jᵀλ` — the part of the Lagrangian gradient the constraint
670/// duals can account for, before any bound multiplier enters.
671fn lagrangian_gradient(
672    sign: Number,
673    grad_f: &[Number],
674    irow: &[i32],
675    jcol: &[i32],
676    jval: &[Number],
677    fortran: bool,
678    lambda: &[Number],
679) -> Vec<Number> {
680    let n = grad_f.len();
681    let off = if fortran { 1 } else { 0 };
682    let mut s = grad_f.to_vec();
683    for k in 0..jval.len() {
684        let row = (irow[k] as usize).wrapping_sub(off);
685        let col = (jcol[k] as usize).wrapping_sub(off);
686        if row < lambda.len() && col < n {
687            s[col] += sign * jval[k] * lambda[row];
688        }
689    }
690    s
691}
692
693/// Bound-**projected** stationarity (a.k.a. "dual infeasibility"): for each
694/// variable, the part of `s` that a valid sign-constrained bound multiplier
695/// `z_L, z_U ≥ 0` cannot absorb. Returns `‖projected s‖∞`.
696///
697/// Stationarity is `z_L − z_U = s`, so a positive `s_j` is absorbable by the
698/// lower bound and a negative one by the upper. A bound may absorb it when
699/// the multiplier that would take is complementary with the variable's gap
700/// to that bound: `|s_j|·gap ≤ compl_tol`, or the gap is within `1e-8`
701/// relative. The product is the test an interior-point answer passes. It
702/// stops a barrier's distance from the bound, not on it, and that distance
703/// is `μ/z`, so an activity test on the gap alone misreads its active bounds
704/// as interior variables and reports their whole gradient as residual. On
705/// PGLib case6468_rte, variables 1.9e-8 inside a bound of 0.2189 (under the
706/// relative test's 1.2e-8) read as a residual of 0.44 against a solver dual
707/// infeasibility of 2.3e-8, with their `ipopt_zU_out` matching the gradient
708/// to every printed digit.
709///
710/// This is a *relaxation*: it projects out exactly the component a bound
711/// multiplier would carry, so it cannot see a missing or wrong `z` (gh #495).
712/// When the `.sol` exports the multipliers, prefer
713/// [`exact_dual_infeasibility`].
714fn bound_projected_residual(
715    s: &[Number],
716    x: &[Number],
717    x_l: &[Number],
718    x_u: &[Number],
719    compl_tol: Number,
720) -> Number {
721    let n = s.len();
722    let mut dual_inf = 0.0_f64;
723    for j in 0..n {
724        let fixed = lower_bound_present(x_l[j])
725            && upper_bound_present(x_u[j])
726            && (x_u[j] - x_l[j]).abs() <= 1e-12;
727        if fixed {
728            continue;
729        }
730        let absorbable = |present: bool, bound: Number, gap: Number, z: Number| {
731            present && (gap.abs() <= 1e-8 * (1.0 + bound.abs()) || z * gap.max(0.0) <= compl_tol)
732        };
733        let r = if s[j] > 0.0 {
734            if absorbable(lower_bound_present(x_l[j]), x_l[j], x[j] - x_l[j], s[j]) {
735                0.0
736            } else {
737                s[j]
738            }
739        } else if absorbable(upper_bound_present(x_u[j]), x_u[j], x_u[j] - x[j], -s[j]) {
740            0.0
741        } else {
742            -s[j]
743        };
744        dual_inf = dual_inf.max(r);
745    }
746    dual_inf
747}
748
749/// Exact dual infeasibility `‖s − (z_L^suffix + z_U^suffix)‖∞`, with `s` the
750/// [`lagrangian_gradient`] at the sign matching the `.sol`'s dual convention.
751///
752/// Stationarity in pounce's internal convention is
753/// `∇f + Jᵀλ − z_L + z_U = 0` with `z_L, z_U ≥ 0`, and the `.sol` suffixes
754/// carry `ipopt_zL_out = +z_L`, `ipopt_zU_out = −z_U` — both equal to the
755/// objective-gradient component at the bound, matching Ipopt 3.14 (gh #296).
756/// So `−z_L + z_U` is exactly `−(zL_out + zU_out)`, and no sign has to be
757/// guessed here beyond the one already chosen for `λ`.
758///
759/// Unlike [`bound_projected_residual`] this sees a bound multiplier that is
760/// missing or wrong, because nothing is projected away.
761fn exact_dual_infeasibility(
762    s: &[Number],
763    z_l_suf: &[Number],
764    z_u_suf: &[Number],
765    obj_sigma: Number,
766) -> Number {
767    let mut dual_inf = 0.0_f64;
768    for (j, &s_j) in s.iter().enumerate() {
769        let z = obj_sigma
770            * (z_l_suf.get(j).copied().unwrap_or(0.0) + z_u_suf.get(j).copied().unwrap_or(0.0));
771        dual_inf = dual_inf.max((s_j - z).abs());
772    }
773    dual_inf
774}
775
776/// Bound complementarity over **variables**:
777/// `max_j max(|z_L·(x−x_L)|, |z_U·(x_U−x)|)` — the quantity Ipopt prints as
778/// `Complementarity` (gh #516). Only variables with a finite bound on the
779/// side in question contribute.
780///
781/// Magnitudes throughout, so the result does not depend on which sign
782/// convention the writer used for the multipliers, nor on which side of a
783/// bound the point sits.
784fn bound_complementarity(
785    x: &[Number],
786    x_l: &[Number],
787    x_u: &[Number],
788    z_l_suf: &[Number],
789    z_u_suf: &[Number],
790) -> Number {
791    let mut comp = 0.0_f64;
792    for j in 0..x.len() {
793        if lower_bound_present(x_l[j]) {
794            let z = z_l_suf.get(j).copied().unwrap_or(0.0);
795            comp = comp.max((z * (x[j] - x_l[j])).abs());
796        }
797        if upper_bound_present(x_u[j]) {
798            let z = z_u_suf.get(j).copied().unwrap_or(0.0);
799            comp = comp.max((z * (x_u[j] - x[j])).abs());
800        }
801    }
802    comp
803}
804
805/// `max_i |λ_i| · dist(g_i, active side)` over constraints with a finite
806/// range — a constraint with a nonzero multiplier should be active.
807/// Equalities (`g_l == g_u`) contribute 0. Best-effort, informational.
808///
809/// This is **constraint** complementarity, over rows, and is not the
810/// quantity a solver reports as `Complementarity` — that one is
811/// [`bound_complementarity`], over variables. The two are unrelated in
812/// magnitude; see the module docs (gh #516).
813fn row_complementarity(lambda: &[Number], g: &[Number], g_l: &[Number], g_u: &[Number]) -> Number {
814    let mut comp = 0.0_f64;
815    for i in 0..lambda.len() {
816        // An equality needs *both* bounds present (gh #403): `g_l = g_u = -5e20`
817        // is the one-sided `g <= -5e20`, not an equality at `-5e20`, and
818        // skipping it here would drop a real complementarity term.
819        if lower_bound_present(g_l[i])
820            && upper_bound_present(g_u[i])
821            && (g_u[i] - g_l[i]).abs() <= 1e-12
822        {
823            continue; // equality: multiplier is free, no complementarity
824        }
825        let dl = if lower_bound_present(g_l[i]) {
826            (g[i] - g_l[i]).abs()
827        } else {
828            Number::INFINITY
829        };
830        let du = if upper_bound_present(g_u[i]) {
831            (g_u[i] - g[i]).abs()
832        } else {
833            Number::INFINITY
834        };
835        let dist = dl.min(du);
836        if dist.is_finite() {
837            comp = comp.max(lambda[i].abs() * dist);
838        }
839    }
840    comp
841}
842
843pub(crate) fn name_at(names: &[String], i: usize, kind: char) -> String {
844    match names.get(i) {
845        Some(s) if !s.is_empty() => s.clone(),
846        _ => format!("{kind}[{i}]"),
847    }
848}
849
850// ---------------------------------------------------------------------------
851// AMPL .sol parser (the inverse of `crate::nl_writer`).
852// ---------------------------------------------------------------------------
853
854#[derive(Debug)]
855struct ParsedSol {
856    x: Vec<Number>,
857    lambda: Vec<Number>,
858    solve_result_num: Option<i32>,
859    /// `ipopt_zL_out` variable suffix, densified to `n`, when present.
860    z_l: Option<Vec<Number>>,
861    /// `ipopt_zU_out` variable suffix, densified to `n`, when present.
862    z_u: Option<Vec<Number>>,
863}
864
865/// Parse the ASCII AMPL `.sol` form pounce writes: a free-text banner, a
866/// blank line, `Options`, an option count + that many option words, the
867/// four-integer count block `<n_dual> <m> <n_primal> <n>`, then the dual
868/// block followed by the primal block, then an optional `objno` line and any
869/// number of suffix blocks.
870fn parse_sol(text: &str) -> Result<ParsedSol, String> {
871    // Find the "Options" delimiter line, then tokenize everything after it.
872    let mut after_options = None;
873    for (i, line) in text.lines().enumerate() {
874        if line.trim() == "Options" {
875            after_options = Some(i);
876            break;
877        }
878    }
879    let start = after_options.ok_or("malformed .sol: no `Options` section found")?;
880    let tail: String = text.lines().skip(start + 1).collect::<Vec<_>>().join(" ");
881    let mut toks = tail.split_whitespace();
882
883    let nopts: usize = toks
884        .next()
885        .ok_or("malformed .sol: missing option count")?
886        .parse()
887        .map_err(|e| format!("malformed .sol: bad option count: {e}"))?;
888    for _ in 0..nopts {
889        toks.next()
890            .ok_or("malformed .sol: truncated option words")?;
891    }
892
893    let next_usize = |toks: &mut std::str::SplitWhitespace, what: &str| -> Result<usize, String> {
894        toks.next()
895            .ok_or_else(|| format!("malformed .sol: missing {what}"))?
896            .parse::<usize>()
897            .map_err(|e| format!("malformed .sol: bad {what}: {e}"))
898    };
899    let n_dual = next_usize(&mut toks, "dual count")?;
900    let _m = next_usize(&mut toks, "constraint count")?;
901    let n_primal = next_usize(&mut toks, "primal count")?;
902    let _n = next_usize(&mut toks, "variable count")?;
903
904    let mut lambda = Vec::with_capacity(n_dual);
905    for k in 0..n_dual {
906        let t = toks
907            .next()
908            .ok_or_else(|| format!("malformed .sol: truncated dual block at {k}"))?;
909        lambda.push(
910            t.parse::<Number>()
911                .map_err(|e| format!("malformed .sol: bad dual {k}: {e}"))?,
912        );
913    }
914    let mut x = Vec::with_capacity(n_primal);
915    for k in 0..n_primal {
916        let t = toks
917            .next()
918            .ok_or_else(|| format!("malformed .sol: truncated primal block at {k}"))?;
919        x.push(
920            t.parse::<Number>()
921                .map_err(|e| format!("malformed .sol: bad primal {k}: {e}"))?,
922        );
923    }
924
925    // Trailing section: an optional `objno <objno> <solve_result_num>` and
926    // any number of suffix blocks.
927    let rest: Vec<&str> = toks.collect();
928    let (solve_result_num, var_suffixes) = parse_sol_tail(&rest, n_primal);
929    let suffix = |name: &str| -> Option<Vec<Number>> {
930        var_suffixes
931            .iter()
932            .find(|(n, _)| n == name)
933            .map(|(_, v)| v.clone())
934    };
935
936    Ok(ParsedSol {
937        x,
938        lambda,
939        solve_result_num,
940        z_l: suffix("ipopt_zL_out"),
941        z_u: suffix("ipopt_zU_out"),
942    })
943}
944
945/// Walk the tokens after the primal block: an optional
946/// `objno <objno> <solve_result_num>` and any number of suffix blocks, each
947/// `suffix <kind> <nvalues> <namelen> <tablen> <tabline>`, the name on its
948/// own line, then `<idx> <value>` pairs (see `pounce_nl::sol_writer` for the
949/// shape pounce writes and Ipopt's AMPL interface writes back).
950///
951/// Returns the `solve_result_num` and every **variable-indexed real**
952/// suffix, densified to `n` — a `.sol` sparse-trims zero entries, so an
953/// absent index means zero, not missing.
954///
955/// A malformed or unsupported block stops the walk and keeps what was read
956/// so far: a `.sol` is still perfectly usable for the feasibility check that
957/// is this tool's actual gate, and a parse error there must not turn a
958/// checkable solution into an I/O failure.
959fn parse_sol_tail(rest: &[&str], n: usize) -> (Option<i32>, Vec<(String, Vec<Number>)>) {
960    let mut solve_result_num = None;
961    let mut out: Vec<(String, Vec<Number>)> = Vec::new();
962    let mut i = 0;
963    while i < rest.len() {
964        match rest[i] {
965            "objno" => {
966                solve_result_num = rest.get(i + 2).and_then(|t| t.parse::<i32>().ok());
967                i += 3;
968            }
969            "suffix" => {
970                let int_at = |k: usize| rest.get(i + k).and_then(|t| t.parse::<i64>().ok());
971                let (Some(kind), Some(nvalues), Some(tablen)) = (int_at(1), int_at(2), int_at(4))
972                else {
973                    break;
974                };
975                let (Some(name), true) = (rest.get(i + 6), nvalues >= 0) else {
976                    break;
977                };
978                let name = (*name).to_string();
979                i += 7;
980                // A suffix value table follows the name as free text we
981                // cannot delimit by whitespace, so its tokens would be
982                // mis-read as values. Neither pounce nor Ipopt writes one.
983                if tablen != 0 {
984                    break;
985                }
986                // Low two bits pick the target (0 = var), 0x4 flags a real
987                // payload — ASL's `ASL_Sufkind_*` bits.
988                let want = (kind & 0x3) == 0 && (kind & 0x4) != 0;
989                let mut dense = vec![0.0; n];
990                let mut complete = true;
991                for _ in 0..nvalues as usize {
992                    let (Some(it), Some(vt)) = (rest.get(i), rest.get(i + 1)) else {
993                        complete = false;
994                        break;
995                    };
996                    if let (true, Ok(idx), Ok(v)) =
997                        (want, it.parse::<usize>(), vt.parse::<Number>())
998                        && idx < n
999                    {
1000                        dense[idx] = v;
1001                    }
1002                    i += 2;
1003                }
1004                if !complete {
1005                    break;
1006                }
1007                if want {
1008                    out.push((name, dense));
1009                }
1010            }
1011            _ => i += 1,
1012        }
1013    }
1014    (solve_result_num, out)
1015}
1016
1017// ---------------------------------------------------------------------------
1018// Console + JSON rendering.
1019// ---------------------------------------------------------------------------
1020
1021fn print_report(args: &VerifyArgs, o: &VerifyOutcome) {
1022    println!("pounce verify — independent solution check");
1023    println!(
1024        "  problem : {}  ({} vars, {} cons)",
1025        args.nl.display(),
1026        o.n_vars,
1027        o.n_cons
1028    );
1029    println!("            sha256:{}", o.nl_sha256);
1030    println!("  solution: {}", args.sol.display());
1031    println!("            sha256:{}", o.sol_sha256);
1032    if let Some(srn) = o.solve_result_num {
1033        println!("  claimed solve_result_num: {srn}");
1034    }
1035    println!();
1036    match o.abs_feas_tol {
1037        Some(t) => println!("  feasibility (absolute tol {t:.1e}):"),
1038        None => println!(
1039            "  feasibility (tol {:.1e}, relative per row: |viol| <= tol*max(1,|row|)):",
1040            o.feas_tol
1041        ),
1042    }
1043    print_row(
1044        "max constraint violation",
1045        o.max_con_violation,
1046        &o.worst_con,
1047    );
1048    print_row(
1049        "max bound violation     ",
1050        o.max_bound_violation,
1051        &o.worst_bound,
1052    );
1053    match (o.discrete.filter(|d| d.total() > 0), o.integrality) {
1054        (Some(d), Some((v, j))) => println!(
1055            "    integrality ({} discrete columns: {} binary + {} integer linear, {} nonlinear): \
1056             max distance to an integer {v:.3e} at x[{j}]",
1057            d.total(),
1058            d.binary,
1059            d.integer,
1060            d.nonlinear_discrete()
1061        ),
1062        (Some(d), None) => println!(
1063            "    integrality: NOT CHECKED — the .nl declares {} discrete variable(s) ({} appearing \
1064             nonlinearly) whose columns cannot be identified from the header",
1065            d.total(),
1066            d.nonlinear_discrete()
1067        ),
1068        (None, _) => {}
1069    }
1070    if let Some(obj) = o.objective {
1071        println!("  objective at x*: {obj:.10e}");
1072    }
1073    if o.stationarity.is_some() || o.bound_multipliers_present {
1074        let source = match (o.duals_present, o.bound_multipliers_present) {
1075            (true, true) => "duals + bound multipliers supplied",
1076            (true, false) => "duals supplied",
1077            (false, true) => "bound multipliers supplied",
1078            (false, false) => "no rows, so no duals to supply",
1079        };
1080        println!();
1081        println!("  optimality (tol {:.1e}, {source}):", o.opt_tol);
1082        if let Some(s) = o.stationarity {
1083            let sign = o.dual_sign.unwrap_or(1);
1084            println!(
1085                "    KKT stationarity residual (bound-projected)  : {s:.3e}  (dual sign {sign:+})"
1086            );
1087        }
1088        if let Some(s) = o.stationarity_with_bound_multipliers {
1089            println!("    dual infeasibility (with z_L/z_U suffixes)   : {s:.3e}");
1090        }
1091        // Two different residuals answer to "complementarity", and the row
1092        // one is NOT what a solver prints as `Complementarity` — label both
1093        // by what they range over so the numbers cannot be crossed (gh #516).
1094        if let Some(c) = o.constraint_complementarity {
1095            println!("    constraint complementarity (rows, |λ|·slack) : {c:.3e}");
1096        }
1097        match o.bound_complementarity {
1098            Some(c) => println!("    bound complementarity (vars, |z|·slack)      : {c:.3e}"),
1099            None => {
1100                println!(
1101                    "    bound complementarity (vars, |z|·slack)      : not checked \
1102                     — the .sol carries no"
1103                );
1104                println!(
1105                    "      ipopt_zL_out/ipopt_zU_out suffixes. This, not the row line \
1106                     above, is the"
1107                );
1108                println!("      quantity a solver reports as `Complementarity`.");
1109            }
1110        }
1111    } else {
1112        println!();
1113        println!("  optimality: not checked (.sol carried no duals)");
1114    }
1115    println!();
1116    if o.verified && !args.require_optimal {
1117        if let Some(s) = o.stationarity_with_bound_multipliers {
1118            if s > o.opt_tol {
1119                println!(
1120                    "  CAVEAT: dual infeasibility {s:.3e} exceeds --opt-tol {:.1e}; VERIFIED here \
1121                     means FEASIBLE only, not optimal (pass --require-optimal to gate on it).",
1122                    o.opt_tol
1123                );
1124            }
1125        } else if let Some(s) = o.stationarity {
1126            if s > o.opt_tol {
1127                println!(
1128                    "  CAVEAT: KKT stationarity residual {s:.3e} exceeds --opt-tol {:.1e}; \
1129                     VERIFIED here means FEASIBLE only, not optimal (pass --require-optimal to gate on it).",
1130                    o.opt_tol
1131                );
1132            }
1133        }
1134    }
1135    let verdict = if o.verified {
1136        "VERIFIED — solution is feasible for the canonical problem".to_string()
1137    } else if !o.feasible
1138        && o.integrality
1139            .is_some_and(|(d, _)| d > o.abs_feas_tol.unwrap_or(o.feas_tol))
1140        && o.max_con_violation <= o.abs_feas_tol.unwrap_or(o.feas_tol)
1141        && o.max_bound_violation <= o.abs_feas_tol.unwrap_or(o.feas_tol)
1142    {
1143        "REJECTED — solution is not integer-feasible (declared integer columns are fractional)"
1144            .to_string()
1145    } else if !o.feasible {
1146        "REJECTED — solution VIOLATES the canonical constraints".to_string()
1147    } else if o.optimal.is_none() {
1148        // Feasible, --require-optimal was asked for, but optimality could not
1149        // be checked at all because the .sol carried no duals — say so rather
1150        // than implying we found it non-optimal.
1151        "REJECTED — feasible, but --require-optimal needs duals and the .sol \
1152         carried none"
1153            .to_string()
1154    } else {
1155        "REJECTED — feasible but not first-order optimal (--require-optimal)".to_string()
1156    };
1157    println!("  VERDICT: {verdict}");
1158}
1159
1160fn print_row(label: &str, v: Number, worst: &Option<RowReport>) {
1161    match worst {
1162        Some(r) => println!(
1163            "    {label}: {v:.3e}  at {} (value {:.6e}, bounds [{:.6e}, {:.6e}])",
1164            r.name, r.value, r.lo, r.hi
1165        ),
1166        None => println!("    {label}: {v:.3e}"),
1167    }
1168}
1169
1170/// Environment variable holding the HMAC key. When set (non-empty) and a
1171/// `--json-output` receipt is requested, the receipt is signed.
1172pub const KEY_ENV: &str = "POUNCE_VERIFY_KEY";
1173
1174/// The exact byte string that gets HMAC-signed. Deliberately **float-free**
1175/// — only hex hashes, integer counts, and the verdict — so any language
1176/// reproduces it byte-for-byte (no float-formatting parity problems between
1177/// Rust and a Python/JS consumer). One `key=value` per line, fixed order,
1178/// trailing newline. The consumer re-derives this from the receipt fields,
1179/// recomputes `HMAC-SHA256(key, preimage)`, and compares to `signature`.
1180/// Documented in `docs/src/verify.md`.
1181///
1182/// The signed fields are exactly the security-critical bindings: *which*
1183/// problem (`nl_sha256`), *which* solution (`sol_sha256`), the problem
1184/// dimensions, and the verdict. The numeric violations in the receipt are
1185/// supporting evidence; trust flows from the hashes + `verified` flag.
1186pub fn signing_preimage(o: &VerifyOutcome) -> String {
1187    format!(
1188        "pounce-verify-receipt/v1\n\
1189         verify_version=1\n\
1190         nl_sha256={}\n\
1191         sol_sha256={}\n\
1192         n_vars={}\n\
1193         n_cons={}\n\
1194         feasible={}\n\
1195         verified={}\n\
1196         verdict={}\n",
1197        o.nl_sha256,
1198        o.sol_sha256,
1199        o.n_vars,
1200        o.n_cons,
1201        o.feasible,
1202        o.verified,
1203        if o.verified { "VERIFIED" } else { "REJECTED" },
1204    )
1205}
1206
1207fn receipt_json(args: &VerifyArgs, o: &VerifyOutcome) -> String {
1208    use serde_json::json;
1209    let worst_con = o.worst_con.as_ref().map(row_json);
1210    let worst_bound = o.worst_bound.as_ref().map(row_json);
1211    let optimality = if o.duals_present || o.bound_multipliers_present {
1212        // Optimality is a property of a FEASIBLE point, so this must not report
1213        // `true` for one that violates the constraints. The stationarity
1214        // residual of an infeasible point can be legitimately zero, which
1215        // previously surfaced as `optimality.optimal: true` inside a receipt
1216        // whose verdict was REJECTED — the top-level fields were correct, but a
1217        // consumer reading this nested field alone was told the opposite.
1218        // The raw residuals are still reported unconditioned: they are useful
1219        // for diagnosing *why* a point failed.
1220        let optimal = o.optimal.map(|opt| opt && o.feasible);
1221        json!({
1222            "available": true,
1223            "objective": o.objective,
1224            "stationarity_residual": o.stationarity,
1225            "dual_sign": o.dual_sign,
1226            "stationarity_residual_with_bound_multipliers":
1227                o.stationarity_with_bound_multipliers,
1228            "constraint_complementarity_residual": o.constraint_complementarity,
1229            "bound_complementarity_residual": o.bound_complementarity,
1230            "bound_multipliers_present": o.bound_multipliers_present,
1231            // Deprecated alias, kept so a v1 consumer does not break. Its bare
1232            // name is the trap gh #516 is about: read
1233            // `constraint_complementarity_residual` instead.
1234            "complementarity_residual": o.constraint_complementarity,
1235            "optimal": optimal,
1236            "note": "`stationarity_residual` is the BOUND-PROJECTED dual infeasibility from \
1237                     the .sol's constraint duals, with bound multipliers inferred from \
1238                     activity; the sign is chosen to match the supplied dual convention. \
1239                     `constraint_complementarity_residual` is max_i |lambda_i| * dist(g_i, \
1240                     nearest finite side) over ROWS — it is NOT what a solver reports as \
1241                     `Complementarity`. That is `bound_complementarity_residual`, \
1242                     max_j max(|z_L*(x-x_L)|, |z_U*(x_U-x)|) over VARIABLES, available only \
1243                     when the .sol carries the ipopt_zL_out/ipopt_zU_out suffixes (null \
1244                     otherwise, meaning not checked — not zero). When those suffixes are \
1245                     present, `stationarity_residual_with_bound_multipliers` is the exact, \
1246                     unprojected residual and is what `--require-optimal` gates on. \
1247                     `complementarity_residual` is a deprecated alias of \
1248                     `constraint_complementarity_residual`. Feasibility is the rigorous \
1249                     gate, and `optimal` is reported false for an infeasible point \
1250                     regardless of its stationarity residual."
1251        })
1252    } else {
1253        json!({ "available": false })
1254    };
1255    let mut receipt = json!({
1256        "pounce_verify_version": 1,
1257        "solver": format!("pounce {}", env!("CARGO_PKG_VERSION")),
1258        "problem": {
1259            "path": args.nl.display().to_string(),
1260            "sha256": o.nl_sha256,
1261            "n_vars": o.n_vars,
1262            "n_cons": o.n_cons,
1263        },
1264        "solution": {
1265            "path": args.sol.display().to_string(),
1266            "sha256": o.sol_sha256,
1267            "claimed_solve_result_num": o.solve_result_num,
1268            "duals_present": o.duals_present,
1269        },
1270        "tolerances": {
1271            "feasibility": o.feas_tol,
1272            "absolute_feasibility": o.abs_feas_tol,
1273            "feasibility_mode": if o.abs_feas_tol.is_some() { "absolute" } else { "relative_per_row" },
1274            "optimality": o.opt_tol,
1275        },
1276        "integrality": match (o.discrete, o.integrality) {
1277            (Some(d), Some((v, j))) if d.total() > 0 => json!({
1278                "checked": true, "binary": d.binary, "integer": d.integer,
1279                "nonlinear": d.nonlinear_discrete(),
1280                "max_distance_to_integer": v, "worst_column": j,
1281            }),
1282            (Some(d), None) if d.total() > 0 => json!({
1283                "checked": false, "binary": d.binary, "integer": d.integer,
1284                "nonlinear": d.nonlinear_discrete(),
1285            }),
1286            _ => json!({ "checked": false, "binary": 0, "integer": 0, "nonlinear": 0 }),
1287        },
1288        "feasibility": {
1289            "max_constraint_violation": o.max_con_violation,
1290            "worst_constraint": worst_con,
1291            "max_bound_violation": o.max_bound_violation,
1292            "worst_bound": worst_bound,
1293            "feasible": o.feasible,
1294        },
1295        "optimality": optimality,
1296        "verdict": if o.verified { "VERIFIED" } else { "REJECTED" },
1297        "verified": o.verified,
1298    });
1299
1300    // Sign the receipt when a key is present. The signature covers the
1301    // float-free `signing_preimage`, NOT the pretty JSON, so a consumer in
1302    // any language can recompute it without matching float formatting.
1303    if let Ok(key) = std::env::var(KEY_ENV) {
1304        if !key.is_empty() {
1305            if let Some(obj) = receipt.as_object_mut() {
1306                let sig = sha256::hmac_hex(key.as_bytes(), signing_preimage(o).as_bytes());
1307                obj.insert("signature_alg".into(), json!("HMAC-SHA256"));
1308                obj.insert(
1309                    "signed_fields".into(),
1310                    json!([
1311                        "verify_version",
1312                        "nl_sha256",
1313                        "sol_sha256",
1314                        "n_vars",
1315                        "n_cons",
1316                        "feasible",
1317                        "verified",
1318                        "verdict"
1319                    ]),
1320                );
1321                obj.insert("signature".into(), json!(sig));
1322            }
1323        }
1324    }
1325
1326    serde_json::to_string_pretty(&receipt).unwrap_or_else(|_| "{}".to_string())
1327}
1328
1329fn row_json(r: &RowReport) -> serde_json::Value {
1330    serde_json::json!({
1331        "index": r.index,
1332        "name": r.name,
1333        "value": r.value,
1334        "lower": r.lo,
1335        "upper": r.hi,
1336        "violation": r.violation,
1337    })
1338}
1339
1340// ---------------------------------------------------------------------------
1341// Self-contained SHA-256 (FIPS 180-4) — content-addresses the receipt's
1342// inputs with zero new dependencies, matching the crate's hand-rolled,
1343// dependency-light style. Known-answer tested below.
1344// ---------------------------------------------------------------------------
1345
1346pub mod sha256 {
1347    const K: [u32; 64] = [
1348        0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1, 0x923f82a4,
1349        0xab1c5ed5, 0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3, 0x72be5d74, 0x80deb1fe,
1350        0x9bdc06a7, 0xc19bf174, 0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc, 0x2de92c6f,
1351        0x4a7484aa, 0x5cb0a9dc, 0x76f988da, 0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7,
1352        0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967, 0x27b70a85, 0x2e1b2138, 0x4d2c6dfc,
1353        0x53380d13, 0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85, 0xa2bfe8a1, 0xa81a664b,
1354        0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070, 0x19a4c116,
1355        0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
1356        0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208, 0x90befffa, 0xa4506ceb, 0xbef9a3f7,
1357        0xc67178f2,
1358    ];
1359
1360    /// Raw 32-byte SHA-256 digest.
1361    pub fn digest(data: &[u8]) -> [u8; 32] {
1362        let mut h: [u32; 8] = [
1363            0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a, 0x510e527f, 0x9b05688c, 0x1f83d9ab,
1364            0x5be0cd19,
1365        ];
1366
1367        // Pad: message || 0x80 || 0x00... || 64-bit big-endian bit length.
1368        let bit_len = (data.len() as u64).wrapping_mul(8);
1369        let mut msg = data.to_vec();
1370        msg.push(0x80);
1371        while msg.len() % 64 != 56 {
1372            msg.push(0);
1373        }
1374        msg.extend_from_slice(&bit_len.to_be_bytes());
1375
1376        let mut w = [0u32; 64];
1377        for chunk in msg.chunks_exact(64) {
1378            for i in 0..16 {
1379                w[i] = u32::from_be_bytes([
1380                    chunk[4 * i],
1381                    chunk[4 * i + 1],
1382                    chunk[4 * i + 2],
1383                    chunk[4 * i + 3],
1384                ]);
1385            }
1386            for i in 16..64 {
1387                let s0 = w[i - 15].rotate_right(7) ^ w[i - 15].rotate_right(18) ^ (w[i - 15] >> 3);
1388                let s1 = w[i - 2].rotate_right(17) ^ w[i - 2].rotate_right(19) ^ (w[i - 2] >> 10);
1389                w[i] = w[i - 16]
1390                    .wrapping_add(s0)
1391                    .wrapping_add(w[i - 7])
1392                    .wrapping_add(s1);
1393            }
1394
1395            let (mut a, mut b, mut c, mut d, mut e, mut f, mut g, mut hh) =
1396                (h[0], h[1], h[2], h[3], h[4], h[5], h[6], h[7]);
1397            for i in 0..64 {
1398                let s1 = e.rotate_right(6) ^ e.rotate_right(11) ^ e.rotate_right(25);
1399                let ch = (e & f) ^ ((!e) & g);
1400                let t1 = hh
1401                    .wrapping_add(s1)
1402                    .wrapping_add(ch)
1403                    .wrapping_add(K[i])
1404                    .wrapping_add(w[i]);
1405                let s0 = a.rotate_right(2) ^ a.rotate_right(13) ^ a.rotate_right(22);
1406                let maj = (a & b) ^ (a & c) ^ (b & c);
1407                let t2 = s0.wrapping_add(maj);
1408                hh = g;
1409                g = f;
1410                f = e;
1411                e = d.wrapping_add(t1);
1412                d = c;
1413                c = b;
1414                b = a;
1415                a = t1.wrapping_add(t2);
1416            }
1417            h[0] = h[0].wrapping_add(a);
1418            h[1] = h[1].wrapping_add(b);
1419            h[2] = h[2].wrapping_add(c);
1420            h[3] = h[3].wrapping_add(d);
1421            h[4] = h[4].wrapping_add(e);
1422            h[5] = h[5].wrapping_add(f);
1423            h[6] = h[6].wrapping_add(g);
1424            h[7] = h[7].wrapping_add(hh);
1425        }
1426
1427        let mut out = [0u8; 32];
1428        for (i, word) in h.iter().enumerate() {
1429            out[4 * i..4 * i + 4].copy_from_slice(&word.to_be_bytes());
1430        }
1431        out
1432    }
1433
1434    fn to_hex(bytes: &[u8]) -> String {
1435        let mut out = String::with_capacity(bytes.len() * 2);
1436        for b in bytes {
1437            out.push_str(&format!("{b:02x}"));
1438        }
1439        out
1440    }
1441
1442    /// Lowercase-hex SHA-256 of `data`.
1443    pub fn hex(data: &[u8]) -> String {
1444        to_hex(&digest(data))
1445    }
1446
1447    /// HMAC-SHA256(key, msg) per RFC 2104, raw 32 bytes.
1448    pub fn hmac(key: &[u8], msg: &[u8]) -> [u8; 32] {
1449        const BLOCK: usize = 64;
1450        let mut k = [0u8; BLOCK];
1451        if key.len() > BLOCK {
1452            k[..32].copy_from_slice(&digest(key));
1453        } else {
1454            k[..key.len()].copy_from_slice(key);
1455        }
1456        let mut ipad = [0x36u8; BLOCK];
1457        let mut opad = [0x5cu8; BLOCK];
1458        for i in 0..BLOCK {
1459            ipad[i] ^= k[i];
1460            opad[i] ^= k[i];
1461        }
1462        let mut inner = Vec::with_capacity(BLOCK + msg.len());
1463        inner.extend_from_slice(&ipad);
1464        inner.extend_from_slice(msg);
1465        let inner_digest = digest(&inner);
1466        let mut outer = Vec::with_capacity(BLOCK + 32);
1467        outer.extend_from_slice(&opad);
1468        outer.extend_from_slice(&inner_digest);
1469        digest(&outer)
1470    }
1471
1472    /// HMAC-SHA256 as lowercase hex.
1473    pub fn hmac_hex(key: &[u8], msg: &[u8]) -> String {
1474        to_hex(&hmac(key, msg))
1475    }
1476}
1477
1478#[cfg(test)]
1479mod tests {
1480    use super::*;
1481    use crate::nl_writer::{SolutionFile, format_sol};
1482    use pounce_common::types::{NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF};
1483
1484    #[test]
1485    fn sha256_known_answers() {
1486        // FIPS 180-4 test vectors.
1487        assert_eq!(
1488            sha256::hex(b""),
1489            "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"
1490        );
1491        assert_eq!(
1492            sha256::hex(b"abc"),
1493            "ba7816bf8f01cfea414140de5dae2223b00361a396177a9cb410ff61f20015ad"
1494        );
1495        assert_eq!(
1496            sha256::hex(b"The quick brown fox jumps over the lazy dog"),
1497            "d7a8fbb307d7809469ca9abcb0082e4f8d5651e46d3cdb762d02d0bf37c9e592"
1498        );
1499    }
1500
1501    #[test]
1502    fn hmac_sha256_known_answers() {
1503        // RFC 4231 test case 2.
1504        assert_eq!(
1505            sha256::hmac_hex(b"Jefe", b"what do ya want for nothing?"),
1506            "5bdcc146bf60754e6a042426089575c75a003f089d2739839dec58b964ec3843"
1507        );
1508        // RFC 4231 test case 1: key = 0x0b * 20, data = "Hi There".
1509        assert_eq!(
1510            sha256::hmac_hex(&[0x0b; 20], b"Hi There"),
1511            "b0344c61d8db38535ca8afceaf0bf12b881dc200c9833da726e9376c2e32cff7"
1512        );
1513    }
1514
1515    #[test]
1516    fn parse_sol_round_trips_writer() {
1517        // Writer is the inverse we must match exactly. Derive the banner
1518        // from the crate version so this fixture never goes stale on a
1519        // version bump (the round-trip is agnostic to the exact string).
1520        let message = format!(
1521            "POUNCE {}: Optimal Solution Found",
1522            env!("CARGO_PKG_VERSION")
1523        );
1524        let payload = SolutionFile {
1525            message: &message,
1526            x: &[1.0, 2.5, -0.5, 100.0],
1527            mult_g: &[0.1, -0.2],
1528            solve_result_num: 0,
1529            suffixes: &[],
1530        };
1531        let text = format_sol(&payload);
1532        let parsed = parse_sol(&text).expect("parse");
1533        assert_eq!(parsed.x.len(), 4);
1534        assert_eq!(parsed.lambda.len(), 2);
1535        assert!((parsed.x[1] - 2.5).abs() < 1e-15);
1536        assert!((parsed.x[3] - 100.0).abs() < 1e-12);
1537        // The primal round-trips as an identity, but the dual block does
1538        // NOT: `format_sol` negates pounce's internal multipliers into the
1539        // AMPL marginal convention (gh #271), and `parse_sol` reads the
1540        // file back verbatim. So a `mult_g` of +0.1 must come back as a
1541        // parsed dual of -0.1. Asserting identity here is what previously
1542        // let the sign defect pass unnoticed.
1543        assert!((parsed.lambda[0] + 0.1).abs() < 1e-15);
1544        assert!((parsed.lambda[1] - 0.2).abs() < 1e-15);
1545        assert_eq!(parsed.solve_result_num, Some(0));
1546    }
1547
1548    #[test]
1549    fn parse_sol_handles_no_duals() {
1550        let payload = SolutionFile {
1551            message: "msg",
1552            x: &[3.0, 4.0],
1553            mult_g: &[],
1554            solve_result_num: 200,
1555            suffixes: &[],
1556        };
1557        let text = format_sol(&payload);
1558        let parsed = parse_sol(&text).expect("parse");
1559        assert_eq!(parsed.x, vec![3.0, 4.0]);
1560        assert!(parsed.lambda.is_empty());
1561        assert_eq!(parsed.solve_result_num, Some(200));
1562    }
1563
1564    #[test]
1565    fn box_violation_basic() {
1566        // inside
1567        assert_eq!(box_violation(5.0, 0.0, 10.0), 0.0);
1568        // below lower
1569        assert!((box_violation(-2.0, 0.0, 10.0) - 2.0).abs() < 1e-15);
1570        // above upper
1571        assert!((box_violation(13.0, 0.0, 10.0) - 3.0).abs() < 1e-15);
1572        // one-sided (no upper)
1573        assert_eq!(box_violation(1e30, 0.0, NLP_UPPER_BOUND_INF), 0.0);
1574    }
1575
1576    #[test]
1577    fn box_violation_rejects_non_finite() {
1578        // Regression: a fabricated `.sol` carrying NaN must register an
1579        // infinite violation, not slip through as feasible. Before the
1580        // `is_finite` guard, `NaN.max(_).max(0.0)` collapsed to `0.0`
1581        // (f64::max drops NaN operands) and the checker reported VERIFIED.
1582        assert_eq!(box_violation(Number::NAN, 0.0, 10.0), Number::INFINITY);
1583        // ±∞ pinned at an unbounded variable is not a real point either.
1584        assert_eq!(
1585            box_violation(Number::INFINITY, 0.0, NLP_UPPER_BOUND_INF),
1586            Number::INFINITY
1587        );
1588        assert_eq!(
1589            box_violation(Number::NEG_INFINITY, NLP_LOWER_BOUND_INF, 10.0),
1590            Number::INFINITY
1591        );
1592    }
1593
1594    /// **gh #403.** `verify` exists to be the independent check on a `.sol`.
1595    /// A checker that under-reports is worse than its blast radius suggests.
1596    ///
1597    /// `is_finite_bound` was a *band* membership test —
1598    /// `b > NLP_LOWER_BOUND_INF && b < NLP_UPPER_BOUND_INF` — applied to `lo`
1599    /// and `hi` alike. A real upper bound of `-5e20` failed it, so `above`
1600    /// became `-inf` and the violation read `0.0`: **ACCEPTED for a `.sol` that
1601    /// violates a declared bound.**
1602    #[test]
1603    fn a_bound_past_the_opposite_sentinel_still_scores_a_violation() {
1604        // x <= -5e20, no lower bound. The point 0.0 violates it by 5e20.
1605        let v = box_violation(0.0, NLP_LOWER_BOUND_INF, -5.0e20);
1606        assert_eq!(
1607            v, 5.0e20,
1608            "0 is 5e20 above an upper bound of -5e20; scoring it 0.0 lets a \
1609             fabricated .sol past the feasibility gate"
1610        );
1611        // Mirror: x >= 5e20, no upper bound.
1612        assert_eq!(box_violation(0.0, 5.0e20, NLP_UPPER_BOUND_INF), 5.0e20);
1613        // A point that does satisfy the same bound still scores zero.
1614        assert_eq!(box_violation(-6.0e20, NLP_LOWER_BOUND_INF, -5.0e20), 0.0);
1615        // And the sentinels themselves still mean "no bound".
1616        assert_eq!(
1617            box_violation(1e30, NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF),
1618            0.0
1619        );
1620    }
1621
1622    // -----------------------------------------------------------------
1623    // gh #516 — the two complementarity quantities.
1624    // -----------------------------------------------------------------
1625
1626    /// The bound multipliers reach a `.sol` only as suffix blocks, so the
1627    /// parser has to pick them out of the trailing section — past `objno`
1628    /// and past whatever other suffixes the writer emitted.
1629    #[test]
1630    fn parse_sol_reads_the_bound_multiplier_suffixes() {
1631        use crate::nl_writer::{SolSuffix, SolSuffixTarget, SolSuffixValues};
1632        let payload = SolutionFile {
1633            message: "msg",
1634            x: &[1.0, -1.0, 0.0],
1635            mult_g: &[0.5],
1636            solve_result_num: 0,
1637            suffixes: &[
1638                // An unrelated block first: the walk must step over it.
1639                SolSuffix {
1640                    name: "sens_sol_state_1".to_string(),
1641                    target: SolSuffixTarget::Var,
1642                    values: SolSuffixValues::Real(vec![9.0, 9.0, 9.0]),
1643                },
1644                SolSuffix {
1645                    name: "ipopt_zL_out".to_string(),
1646                    target: SolSuffixTarget::Var,
1647                    values: SolSuffixValues::Real(vec![0.0, 2.0, 0.0]),
1648                },
1649                SolSuffix {
1650                    name: "ipopt_zU_out".to_string(),
1651                    target: SolSuffixTarget::Var,
1652                    values: SolSuffixValues::Real(vec![-4.0, 0.0, 0.0]),
1653                },
1654            ],
1655        };
1656        let parsed = parse_sol(&format_sol(&payload)).expect("parse");
1657        assert_eq!(parsed.solve_result_num, Some(0), "objno still parses");
1658        // Densified back to `n`: the writer sparse-trims zeros, so an absent
1659        // index means zero — not a short vector, and not "missing".
1660        assert_eq!(parsed.z_l, Some(vec![0.0, 2.0, 0.0]));
1661        assert_eq!(parsed.z_u, Some(vec![-4.0, 0.0, 0.0]));
1662    }
1663
1664    /// No suffixes → bound complementarity is *not checked*, and must stay
1665    /// `None` rather than collapse to a comfortable `0.0`.
1666    #[test]
1667    fn parse_sol_reports_absent_bound_multipliers_as_absent() {
1668        let payload = SolutionFile {
1669            message: "msg",
1670            x: &[1.0],
1671            mult_g: &[0.5],
1672            solve_result_num: 0,
1673            suffixes: &[],
1674        };
1675        let parsed = parse_sol(&format_sol(&payload)).expect("parse");
1676        assert!(parsed.z_l.is_none() && parsed.z_u.is_none());
1677    }
1678
1679    /// `min (x−3)² + (y+2)²  s.t.  x ≤ 1, y ≥ −1` — the model whose export
1680    /// convention is pinned in `main.rs` (gh #296): `ipopt_zL_out = +z_L`,
1681    /// `ipopt_zU_out = −z_U`, both equal to `∂f/∂x` at the bound.
1682    ///
1683    /// At the exact optimum every slack is zero, so bound complementarity is
1684    /// zero whichever sign convention the writer used — the check is on
1685    /// magnitudes. Off the optimum it is `|z| · slack`.
1686    #[test]
1687    fn bound_complementarity_is_z_times_slack_over_variables() {
1688        let x_l = [NLP_LOWER_BOUND_INF, -1.0];
1689        let x_u = [1.0, NLP_UPPER_BOUND_INF];
1690        // Exactly on both bounds: no slack anywhere.
1691        assert_eq!(
1692            bound_complementarity(&[1.0, -1.0], &x_l, &x_u, &[0.0, 2.0], &[-4.0, 0.0]),
1693            0.0
1694        );
1695        // Pull x off its upper bound by 1e-3 while keeping z_U: the product
1696        // is the residual, and the sign of the multiplier does not enter.
1697        let c = bound_complementarity(&[0.999, -1.0], &x_l, &x_u, &[0.0, 2.0], &[-4.0, 0.0]);
1698        assert!((c - 4.0e-3).abs() < 1e-12, "got {c}");
1699        let flipped = bound_complementarity(&[0.999, -1.0], &x_l, &x_u, &[0.0, 2.0], &[4.0, 0.0]);
1700        assert_eq!(c, flipped, "magnitudes only — no sign convention assumed");
1701        // A variable with no bound on the side in question contributes
1702        // nothing, however large its (meaningless) multiplier.
1703        assert_eq!(
1704            bound_complementarity(
1705                &[0.0],
1706                &[NLP_LOWER_BOUND_INF],
1707                &[NLP_UPPER_BOUND_INF],
1708                &[1e6],
1709                &[1e6]
1710            ),
1711            0.0
1712        );
1713    }
1714
1715    /// The exact residual uses the multipliers the `.sol` actually carries,
1716    /// so it sees what the bound-projected one projects away — the gh #495
1717    /// blind spot: a bound multiplier that is missing or wrong leaves the
1718    /// projected residual at `0.0`.
1719    #[test]
1720    fn exact_dual_infeasibility_sees_what_the_projection_hides() {
1721        // `min (x−3)² s.t. x ≤ 1`: x* = 1, ∇f = −4, so z_U = 4 and the
1722        // exported suffix is `ipopt_zU_out = −4`.
1723        let s = [-4.0];
1724        let x = [1.0];
1725        let x_l = [NLP_LOWER_BOUND_INF];
1726        let x_u = [1.0];
1727        assert_eq!(exact_dual_infeasibility(&s, &[0.0], &[-4.0], 1.0), 0.0);
1728
1729        // Projection: x sits on its upper bound, so a valid z_U absorbs the
1730        // whole negative gradient and the residual reads zero — with *no*
1731        // multiplier supplied at all.
1732        assert_eq!(bound_projected_residual(&s, &x, &x_l, &x_u, 1e-6), 0.0);
1733        // The exact check does not get to assume one exists.
1734        assert_eq!(exact_dual_infeasibility(&s, &[0.0], &[0.0], 1.0), 4.0);
1735        // Nor that it has the right sign.
1736        assert_eq!(exact_dual_infeasibility(&s, &[0.0], &[4.0], 1.0), 8.0);
1737    }
1738
1739    /// An interior-point answer stops `μ/z` short of an active bound, not on
1740    /// it. Numbers from PGLib case6468_rte: `x` 1.9e-8 inside an upper
1741    /// bound of 0.2189 with gradient −0.4399 is an active bound carrying
1742    /// `z_U = 0.44` at complementarity 8.4e-9 — not an interior variable
1743    /// with a residual of 0.44, which is what a gap-only activity test
1744    /// (1.2e-8 here) reported. The same gradient on a variable well inside
1745    /// its box is still a residual.
1746    #[test]
1747    fn a_barrier_stand_off_from_a_bound_is_still_active() {
1748        let x_l = [-0.1276];
1749        let x_u = [0.2189];
1750        let s = [-0.4399];
1751        let near = [0.2189 - 1.9e-8];
1752        assert_eq!(bound_projected_residual(&s, &near, &x_l, &x_u, 1e-6), 0.0);
1753        let inside = [0.05];
1754        assert_eq!(
1755            bound_projected_residual(&s, &inside, &x_l, &x_u, 1e-6),
1756            0.4399
1757        );
1758        // Only the bound on the matching side absorbs: a positive gradient
1759        // at the upper bound needs a negative z_U and stays a residual.
1760        assert_eq!(
1761            bound_projected_residual(&[0.4399], &near, &x_l, &x_u, 1e-6),
1762            0.4399
1763        );
1764    }
1765
1766    /// **gh #516.** Constraint complementarity (rows) and bound
1767    /// complementarity (variables) are different quantities at the same
1768    /// point, and can disagree by orders of magnitude. Printing either under
1769    /// a bare `complementarity residual` label invites the comparison that
1770    /// cost two people an afternoon in #505; this test pins the fact that
1771    /// makes the label matter.
1772    #[test]
1773    fn row_and_bound_complementarity_are_different_quantities() {
1774        // One inequality row `g ≥ 0`, slack 4.5e-2, multiplier 1 — a real
1775        // row-complementarity residual.
1776        let rows = row_complementarity(&[1.0], &[4.5e-2], &[0.0], &[NLP_UPPER_BOUND_INF]);
1777        assert!((rows - 4.5e-2).abs() < 1e-15);
1778        // The same point's variables sit hard on their bounds: bound
1779        // complementarity is eleven orders of magnitude smaller.
1780        let bounds = bound_complementarity(
1781            &[1.0],
1782            &[NLP_LOWER_BOUND_INF],
1783            &[1.0 + 1e-11],
1784            &[0.0],
1785            &[-1.0],
1786        );
1787        assert!(bounds < 1e-10, "got {bounds}");
1788        assert!(
1789            rows / bounds > 1e8,
1790            "the two must not be read as one number"
1791        );
1792    }
1793
1794    /// The same predicate sizes a row's magnitude for the scale-relative
1795    /// feasibility test. A row written at `5e20` must report that magnitude,
1796    /// not fall back to its evaluated value alone.
1797    #[test]
1798    fn row_magnitude_counts_a_bound_past_the_opposite_sentinel() {
1799        assert_eq!(
1800            row_magnitude(1.0, NLP_LOWER_BOUND_INF, -5.0e20),
1801            5.0e20,
1802            "the row's own upper bound is its magnitude"
1803        );
1804        assert_eq!(row_magnitude(1.0, 5.0e20, NLP_UPPER_BOUND_INF), 5.0e20);
1805        // Absent on both sides: only the evaluated value carries magnitude.
1806        assert_eq!(
1807            row_magnitude(3.0, NLP_LOWER_BOUND_INF, NLP_UPPER_BOUND_INF),
1808            3.0
1809        );
1810    }
1811}