Struct ReductionFrame 

pub struct ReductionFrame {
    pub fixed_vars: Vec<usize>,
    pub fixed_values: Vec<Number>,
    pub dropped_rows: Vec<usize>,
    pub var_map: Vec<Option<usize>>,
    pub row_map: Vec<Option<usize>>,
}

One layer of the postsolve stack. Built once per accepted block elimination by PR 8’s orchestrator.

Fields

fixed_vars: Vec<usize>

Inner-variable indices fixed by this layer, in ascending order.

fixed_values: Vec<Number>

Their values at the block-solve fixed point.

dropped_rows: Vec<usize>

Inner equality-row indices dropped by this layer, in ascending order. dropped_rows.len() == fixed_vars.len().

var_map: Vec<Option<usize>>

var_map[i] = Some(reduced_idx) if i survives this layer, None if i is in fixed_vars.

row_map: Vec<Option<usize>>

Same for rows.

Implementations

impl ReductionFrame

pub fn new( n_vars: usize, n_rows: usize, fixed_vars: Vec<usize>, fixed_values: Vec<Number>, dropped_rows: Vec<usize>, ) -> Self

Build a frame from the (sorted) lists of fixed variables / values / dropped rows and the full-space problem shape.

Examples found in repository

examples/frame_roundtrip.rs (lines 38-44)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (lines 220-226)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

pub fn n_full_vars(&self) -> usize

Examples found in repository

examples/frame_roundtrip.rs (line 49)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

pub fn n_full_rows(&self) -> usize

Examples found in repository

examples/frame_roundtrip.rs (line 54)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

pub fn n_reduced_vars(&self) -> usize

Examples found in repository

examples/frame_roundtrip.rs (line 50)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (line 233)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

pub fn n_reduced_rows(&self) -> usize

Examples found in repository

examples/frame_roundtrip.rs (line 55)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (line 234)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

pub fn project_x(&self, x_full: &[Number]) -> Vec<Number>

Project a full-space x vector into reduced space (drop the fixed entries).

pub fn lift_x(&self, x_reduced: &[Number]) -> Vec<Number>

Lift a reduced x back to full space, splicing the fixed values back into their original positions.

Examples found in repository

examples/frame_roundtrip.rs (line 64)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (line 242)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

pub fn project_lambda(&self, lambda_full: &[Number]) -> Vec<Number>

Project a full-space λ vector into reduced space.

pub fn lift_lambda(&self, lambda_reduced: &[Number]) -> Vec<Number>

Lift a reduced λ back to full space, with zeros at dropped row indices. (Real values for dropped rows come from Self::recover_dropped_multipliers.)

Examples found in repository

examples/frame_roundtrip.rs (line 72)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (line 268)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

pub fn recover_dropped_multipliers( &self, grad_f: &[Number], jac_full_row_major: &[Number], lambda_full: &[Number], ) -> Result<Vec<Number>, BlockSolveError>

Recover the k = fixed_vars.len() dropped-row multipliers via dense LU on the full-space KKT stationarity equations at the fixed variables. Returns one entry per self.dropped_rows (in the same order).

Assumption: fixed variables are interior to their original bounds at the optimum (so z_l = z_u = 0 for them).

Inputs

• grad_f — objective gradient at the full-space optimum (length n_full_vars).
• jac_full_row_major — dense full-space Jacobian (n_full_rows × n_full_vars) at the optimum.
• lambda_full — multipliers for kept rows; entries at dropped-row positions are ignored.

Example

use pounce_presolve::reduction_frame::ReductionFrame;

// 1 var, 1 row, dropped:  c(x) = x - 3 = 0, obj f = 4 x.
// Stationarity:  4 - 1 * λ = 0  →  λ = 4.
let frame = ReductionFrame::new(1, 1, vec![0], vec![3.0], vec![0]);
let grad_f = [4.0];
let jac = [1.0];
let lambda_full = [0.0]; // dropped, ignored
let lam = frame
    .recover_dropped_multipliers(&grad_f, &jac, &lambda_full)
    .unwrap();
assert!((lam[0] - 4.0).abs() < 1e-12);

Examples found in repository

examples/frame_roundtrip.rs (line 86)

fn main() {
    // Full-space problem shape.
    let n_vars = 3;
    let n_rows = 3;

    // Block solve fixed b1, b2.
    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        vec![0, 1],                 // fixed_vars: b1, b2
        vec![2.0 / 3.0, 5.0 / 3.0], // their values
        vec![0, 1],                 // dropped_rows: c0, c1
    );

    println!("frame_roundtrip example — PR 7 of pounce#53");
    println!(
        "full vars: {}, reduced vars: {}",
        frame.n_full_vars(),
        frame.n_reduced_vars()
    );
    println!(
        "full rows: {}, reduced rows: {}",
        frame.n_full_rows(),
        frame.n_reduced_rows()
    );

    // Simulated IPM result on the reduced problem.
    let y_star = 8.0 / 3.0;
    let x_reduced = [y_star];
    let lambda_reduced = [2.0 * y_star]; // stationarity at y: 2y - λ_kept = 0

    // Lift x to full space.
    let x_full = frame.lift_x(&x_reduced);
    println!(
        "lifted x:      [{:.6}, {:.6}, {:.6}]",
        x_full[0], x_full[1], x_full[2]
    );

    // Lift λ to full space — zeros at dropped row indices, to be
    // filled in by multiplier recovery.
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("lifted λ (pre-recovery): {:?}", lambda_full);

    // Build the full Jacobian at x_full (rows × vars, row-major).
    let jac_full = [
        2.0, 1.0, 0.0, // c0: 2 b1 + b2
        1.0, -1.0, 0.0, // c1: b1 - b2
        1.0, 1.0, 1.0, // c2: b1 + b2 + y
    ];
    // Objective gradient at x_full.
    let grad_f = [10.0, 4.0, 2.0 * x_full[2]];

    // Recover the dropped multipliers.
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .expect("non-singular");

    // Fill them into the full λ vector.
    for (k, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[k];
    }
    println!("lifted λ (post-recovery): {:?}", lambda_full);

    // Verify full-space stationarity at ALL variables.
    println!("KKT stationarity residual (should be ≈ 0):");
    let mut max_res = 0.0_f64;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("  var {i}: {:.3e}", s);
        max_res = max_res.max(s.abs());
    }
    println!("max |residual|: {:.3e}", max_res);
    assert!(max_res < 1e-12);
}

examples/pipeline_demo.rs (line 280)

fn main() {
    let n_vars = 4; // a, b, c, d
    let n_rows = 4;
    // Jacobian sparsity (row-major).
    let jac_pattern = [
        (0, 0), // c0: a
        (1, 1),
        (1, 2), // c1: b + c
        (2, 1),
        (2, 2), // c2: b - c
        (3, 0),
        (3, 1),
        (3, 3), // c3: a + b + d
    ];
    let jac_irow: Vec<i32> = jac_pattern.iter().map(|(r, _)| *r as i32).collect();
    let jac_jcol: Vec<i32> = jac_pattern.iter().map(|(_, c)| *c as i32).collect();
    let g_l = [1.0, 5.0, -1.0, 4.0];
    let g_u = [1.0, 5.0, -1.0, 4.0]; // all equalities

    let probe = ProbeView {
        n_vars,
        m_rows: n_rows,
        jac_irow: &jac_irow,
        jac_jcol: &jac_jcol,
        jac_values: None,
        g_l: &g_l,
        g_u: &g_u,
        linearity: None,
        one_based: false,
        eq_tol: 1e-12,
        excluded_vars: None,
        excluded_rows: None,
    };

    println!("pipeline_demo — PRs 2-7 of pounce#53");
    println!("======================================");

    // ----- Stage 1: incidence + matching -----
    let inc = EqualityIncidence::from_probe(&probe);
    println!(
        "\n[Stage 1] EqualityIncidence: {} eq rows, {} vars",
        inc.n_eq_rows(),
        inc.n_vars
    );
    let m = hopcroft_karp(&inc);
    println!("[Stage 1] Matching: size = {}", m.size);

    // ----- Stage 2: DM partition -----
    let dm = DulmageMendelsohnPartition::from_matching(&inc, &m);
    println!("\n[Stage 2] DM partition:");
    println!("   over rows={:?}, cols={:?}", dm.over_rows, dm.over_cols);
    println!(
        "   square rows={:?}, cols={:?}",
        dm.square_rows, dm.square_cols
    );
    println!(
        "   under rows={:?}, cols={:?}",
        dm.under_rows, dm.under_cols
    );

    // ----- Stage 3: components -----
    let comps = SquareComponents::of_square_part(&inc, &m, &dm);
    println!(
        "\n[Stage 3] Components: {} square block(s)",
        comps.components.len()
    );
    for (i, c) in comps.components.iter().enumerate() {
        println!("   component {i}: rows={:?}, cols={:?}", c.eq_rows, c.cols);
    }

    // ----- Stage 4: BTF -----
    println!("\n[Stage 4] BTF (per component):");
    let ineq = InequalityIncidence::from_probe(&probe);

    // Walk each component's BTF, solve each block, build the frame.
    let mut all_fixed_vars: Vec<usize> = Vec::new();
    let mut all_fixed_values: Vec<Number> = Vec::new();
    let mut all_dropped_rows: Vec<usize> = Vec::new();
    let mut x_running = vec![0.0; n_vars];
    // Full A and b for this problem, used both for block extraction
    // and the final KKT check.
    let a_full = [
        [1.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 1.0, 0.0],
        [0.0, 1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0, 1.0],
    ];
    let b_full = [1.0, 5.0, -1.0, 4.0];

    for comp in &comps.components {
        let btf = BlockTriangularForm::of_component(&inc, &m, comp);
        for (bi, block) in btf.blocks.iter().enumerate() {
            let class = classify_block(block, &ineq, &Default::default());
            println!(
                "   component cols={:?}, block {bi}: rows={:?}, cols={:?}, class={class:?}",
                comp.cols, block.eq_rows, block.cols
            );
            assert_eq!(class, AuxiliaryCouplingClass::PureEquality);

            // Build a small linear system for this block, splicing in
            // the values of already-solved variables.
            let k = block.eq_rows.len();
            let mut a_block = vec![0.0; k * k];
            let mut b_block = vec![0.0; k];
            for (ii, &r) in block.eq_rows.iter().enumerate() {
                let mut rhs = b_full[r];
                for j in 0..n_vars {
                    if let Some(jj) = block.cols.iter().position(|&c| c == j) {
                        a_block[ii * k + jj] = a_full[r][j];
                    } else {
                        rhs -= a_full[r][j] * x_running[j];
                    }
                }
                b_block[ii] = rhs;
            }

            // Solve.
            let mut eqs = LinearBlockEqs {
                a: a_block,
                b: b_block,
                n: k,
            };
            let mut solver = DampedNewtonSolver;
            let opt = BlockSolveOptions::default();
            let out = solver.solve(&vec![0.0; k], &mut eqs, &opt).unwrap();

            // Record the block's solution.
            for (ii, &c) in block.cols.iter().enumerate() {
                x_running[c] = out.x[ii];
                all_fixed_vars.push(c);
                all_fixed_values.push(out.x[ii]);
            }
            for &r in &block.eq_rows {
                all_dropped_rows.push(r);
            }
        }
    }

    // ----- Stage 5: lift via the reduction frame -----
    println!("\n[Stage 5] After block solves:");
    println!("   x_running (a, b, c, d) = {:?}", x_running);
    // For this demo we declare ALL fixed vars and rows as one
    // composite frame — in the real PR 8 orchestrator, each
    // elimination layer becomes its own frame on the stack.
    // Sort by variable index for the frame's contract.
    let mut order: Vec<usize> = (0..all_fixed_vars.len()).collect();
    order.sort_by_key(|&i| all_fixed_vars[i]);
    let fixed_vars_sorted: Vec<usize> = order.iter().map(|&i| all_fixed_vars[i]).collect();
    let fixed_values_sorted: Vec<Number> = order.iter().map(|&i| all_fixed_values[i]).collect();
    let mut dropped_rows_sorted = all_dropped_rows.clone();
    dropped_rows_sorted.sort_unstable();

    let frame = ReductionFrame::new(
        n_vars,
        n_rows,
        fixed_vars_sorted.clone(),
        fixed_values_sorted,
        dropped_rows_sorted.clone(),
    );
    println!(
        "   ReductionFrame: fixed_vars={:?}, dropped_rows={:?}",
        frame.fixed_vars, frame.dropped_rows
    );
    println!(
        "   reduced shape: {} vars × {} rows",
        frame.n_reduced_vars(),
        frame.n_reduced_rows()
    );

    // For this problem BTF decomposed every row into a block — the
    // reduced shape is 0×0 and no IPM step is needed. The lift is
    // trivial: x_full = the values the block solves produced.
    let x_reduced: Vec<Number> = Vec::new();
    println!("\n[Stage 6] Reduced problem is empty (presolve solved it all).");
    let x_full = frame.lift_x(&x_reduced);
    println!("   lifted x_full = {:?}", x_full);

    // Sanity: should match (1, 2, 3, 1).
    let expected = [1.0, 2.0, 3.0, 1.0];
    let max_x_err = x_full
        .iter()
        .zip(expected.iter())
        .map(|(a, b)| (a - b).abs())
        .fold(0.0_f64, f64::max);
    println!("   max |x - expected|: {:.3e}", max_x_err);
    assert!(max_x_err < 1e-9);

    // ----- Stage 7: multiplier recovery -----
    // Objective: f = (a-1)² + (b-2)² + (c-3)² + d²
    // ∇f at the optimum (1, 2, 3, 1): (0, 0, 0, 2).
    let grad_f = [
        2.0 * (x_full[0] - 1.0),
        2.0 * (x_full[1] - 2.0),
        2.0 * (x_full[2] - 3.0),
        2.0 * x_full[3],
    ];
    println!("\n[Stage 7] ∇f at the optimum: {:?}", grad_f);

    // No kept rows → reduced λ is empty.
    let lambda_reduced: Vec<Number> = Vec::new();
    let mut lambda_full = frame.lift_lambda(&lambda_reduced);
    println!("   reduced λ = {:?}", lambda_reduced);
    println!("   lifted λ (pre-recovery) = {:?}", lambda_full);

    // Build full Jacobian for the recovery LU.
    let mut jac_full = vec![0.0; n_rows * n_vars];
    for (r, row) in a_full.iter().enumerate() {
        for (c, &v) in row.iter().enumerate() {
            jac_full[r * n_vars + c] = v;
        }
    }
    let lam_dropped = frame
        .recover_dropped_multipliers(&grad_f, &jac_full, &lambda_full)
        .unwrap();
    for (idx, &r) in frame.dropped_rows.iter().enumerate() {
        lambda_full[r] = lam_dropped[idx];
    }
    println!("   recovered λ (full) = {:?}", lambda_full);

    // Final KKT residual at every variable.
    println!("\n[Stage 8] Full-space KKT residual:");
    let mut max_kkt: Number = 0.0;
    for i in 0..n_vars {
        let mut s = grad_f[i];
        for r in 0..n_rows {
            s -= jac_full[r * n_vars + i] * lambda_full[r];
        }
        println!("   ∂L/∂x[{i}] = {:.3e}", s);
        max_kkt = max_kkt.max(s.abs());
    }
    println!("   max |KKT residual| = {:.3e}", max_kkt);
    assert!(max_kkt < 1e-10);
    println!("\n✓ Pipeline end-to-end correct.");
}

Trait Implementations

impl Clone for ReductionFrame

fn clone(&self) -> ReductionFrame

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for ReductionFrame

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Default for ReductionFrame

fn default() -> ReductionFrame

Returns the “default value” for a type.