frame_roundtrip/
frame_roundtrip.rs

1//! Worked example for PR 7 of issue #53 — reduction frame +
2//! postsolve multiplier recovery.
3//!
4//! Full problem (3 variables `b1, b2, y`, 3 constraints):
5//!
6//! ```text
7//!   min   f(b1, b2, y) = 10 b1 + 4 b2 + y²
8//!   s.t.  c0(b):  2 b1 +   b2       - 3 = 0   ← dropped (equality)
9//!         c1(b):    b1 -   b2       + 1 = 0   ← dropped (equality)
10//!         c2(b):    b1 +   b2 +  y  - 5 = 0   ← kept    (equality)
11//! ```
12//!
13//! Block solve fixes `(b1, b2) = (2/3, 5/3)` by satisfying `c0` and
14//! `c1`. The reduced problem has one variable (`y`) and one row
15//! (`c2`), which gives `y* = 8/3`. The IPM then converges to that.
16//!
17//! Postsolve must:
18//! 1. Lift `x` back to full space `(b1*, b2*, y*)`.
19//! 2. Recover λ for the dropped rows by solving the full-space KKT
20//!    stationarity at the fixed variables.
21//! 3. Verify the full-space KKT residual is below tolerance.
22//!
23//! Run with:
24//! ```bash
25//! cargo run -p pounce-presolve --example frame_roundtrip
26//! ```
27
28#![allow(clippy::expect_used, clippy::unwrap_used)]
29
30use pounce_presolve::reduction_frame::ReductionFrame;
31
32fn main() {
33    // Full-space problem shape.
34    let n_vars = 3;
35    let n_rows = 3;
36
37    // Block solve fixed b1, b2.
38    let frame = ReductionFrame::new(
39        n_vars,
40        n_rows,
41        vec![0, 1],                 // fixed_vars: b1, b2
42        vec![2.0 / 3.0, 5.0 / 3.0], // their values
43        vec![0, 1],                 // dropped_rows: c0, c1
44    );
45
46    println!("frame_roundtrip example — PR 7 of pounce#53");
47    println!(
48        "full vars: {}, reduced vars: {}",
49        frame.n_full_vars(),
50        frame.n_reduced_vars()
51    );
52    println!(
53        "full rows: {}, reduced rows: {}",
54        frame.n_full_rows(),
55        frame.n_reduced_rows()
56    );
57
58    // Simulated IPM result on the reduced problem.
59    let y_star = 8.0 / 3.0;
60    let x_reduced = [y_star];
61    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0
62
63    // Lift x to full space.
64    let x_full = frame.lift_x(&x_reduced);
65    println!(
66        "lifted x:      [{:.6}, {:.6}, {:.6}]",
67        x_full[0], x_full[1], x_full[2]
68    );
69
70    // Lift λ to full space — zeros at dropped row indices, to be
71    // filled in by multiplier recovery.
72    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
73    println!("lifted λ (pre-recovery): {:?}", lambda_full);
74
75    // Build the full Jacobian at x_full (rows × vars, row-major).
76    let jac_full = [
77        2.0, 1.0, 0.0, // c0: 2 b1 + b2
78        1.0, -1.0, 0.0, // c1: b1 - b2
79        1.0, 1.0, 1.0, // c2: b1 + b2 + y
80    ];
81    // Objective gradient at x_full.
82    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];
83
84    // Recover the dropped multipliers.
85    let lam_dropped = frame
86        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
87        .expect("non-singular");
88
89    // Fill them into the full λ vector.
90    for (k, &r) in frame.dropped_rows.iter().enumerate() {
91        lambda_full[r] = lam_dropped[k];
92    }
93    println!("lifted λ (post-recovery): {:?}", lambda_full);
94
95    // Verify full-space stationarity at ALL variables.
96    println!("KKT stationarity residual (should be ≈ 0):");
97    let mut max_res = 0.0_f64;
98    for i in 0..n_vars {
99        let mut s = grad_f[i];
100        for r in 0..n_rows {
101            s -= jac_full[r * n_vars + i] * lambda_full[r];
102        }
103        println!("  var {i}: {:.3e}", s);
104        max_res = max_res.max(s.abs());
105    }
106    println!("max |residual|: {:.3e}", max_res);
107    assert!(max_res < 1e-12);
108}
frame_roundtrip/frame_roundtrip.rs

frame_roundtrip/
frame_roundtrip.rs
