gam-sae 0.3.155

//! Streaming block updates must use one coordinated gamma/frame state, descend
//! in their conditional objective, and certify the model actually returned.

use super::BlockSparseStreamState;
use crate::sparse_dict::BlockSparseConfig;
use ndarray::{Array2, array};

fn coupled_fixture() -> (Array2<f32>, Array2<f32>, BlockSparseConfig) {
    let x = array![[1.0_f32, 2.0], [-2.0, 1.0], [0.5, 0.2], [-0.7, -0.4]];
    let decoder = array![[1.0_f32, 0.0], [0.6, 0.8]];
    let config = BlockSparseConfig {
        n_blocks: 2,
        block_size: 1,
        block_topk: 2,
        max_epochs: 64,
        minibatch: 4,
        block_tile: 2,
        frame_ridge: 0.0,
        aux_k: 0,
        matryoshka_prefix: false,
        tolerance: 1e-6,
    };
    (x, decoder, config)
}

#[test]
fn frame_proposals_and_certificates_are_identical_across_worker_counts() {
    let mut seed = 2826_u64;
    let mut sample = || {
        seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1);
        (seed >> 33) as f32 / 2147483648.0 - 0.5
    };
    let x = Array2::from_shape_fn((37, 8), |_| sample());
    let mut decoder = Array2::from_shape_fn((12, 8), |_| sample());
    for block in 0..5 {
        let mut frame = decoder
            .slice(ndarray::s![block * 2..(block + 1) * 2, ..])
            .to_owned();
        crate::sparse_dict::block::gram_schmidt_rows(&mut frame);
        decoder
            .slice_mut(ndarray::s![block * 2..(block + 1) * 2, ..])
            .assign(&frame);
    }
    decoder.slice_mut(ndarray::s![10.., ..]).fill(0.0);
    let mut config = BlockSparseConfig::new(6, 2);
    config.block_topk = 3;
    config.minibatch = 13;
    config.aux_k = 0;
    let mut reference = None;
    for workers in [1, 4] {
        let pool = rayon::ThreadPoolBuilder::new()
            .num_threads(workers)
            .build()
            .unwrap();
        let observed = pool.install(|| {
            let mut state =
                BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
            let mut epochs = Vec::new();
            for _ in 0..3 {
                state.partial_fit(x.view()).unwrap();
                let stats = state.end_epoch().unwrap();
                assert!(
                    state
                        .decoder
                        .slice(ndarray::s![10.., ..])
                        .iter()
                        .all(|&v| v == 0.0)
                );
                assert!(stats.frame_residual.is_finite());
                epochs.push((
                    state.decoder.clone(),
                    stats.gamma,
                    stats.explained_variance,
                    stats.gamma_residual,
                    stats.frame_residual,
                    stats.converged,
                    state.last_second.clone(),
                ));
            }
            epochs
        });
        if let Some(expected) = reference.as_ref() {
            assert_eq!(
                &observed, expected,
                "worker count changed the proposed model or certificate"
            );
        } else {
            reference = Some(observed);
        }
    }
}

#[test]
fn parallel_stream_moments_match_dense_reference_across_batches_and_shards() {
    let x = Array2::from_shape_fn((17, 4), |(row, feature)| {
        if row == 0 {
            0.0
        } else {
            ((row * 7 + feature * 13) as f32 * 0.17).sin()
        }
    });
    let decoder = array![
        [1.0_f32, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0],
        [0.6, 0.0, 0.0, 0.8],
        [0.0, 0.0, 1.0, 0.0],
    ];
    // Independent dense projector algebra. All three blocks are selected on
    // nonzero rows; the zero row exercises padded slots without phantom usage.
    let weights = x.mapv(f64::from).dot(&decoder.mapv(f64::from).t());
    let total = weights.dot(&decoder.mapv(f64::from));
    let gamma = 0.37_f32;
    let baseline_gamma = 0.61_f32;
    let residual = x.mapv(f64::from) - &total * gamma as f64;
    let baseline_residual = x.mapv(f64::from) - &total * baseline_gamma as f64;
    let expected_rss = residual.iter().map(|v| v * v).sum::<f64>();
    let expected_baseline_rss = baseline_residual.iter().map(|v| v * v).sum::<f64>();
    for threads in [1, 4] {
        let pool = rayon::ThreadPoolBuilder::new()
            .num_threads(threads)
            .build()
            .unwrap();
        pool.install(|| {
            for batch in [1, 5, 17] {
                for shard_rows in [3, 17] {
                    let mut config = BlockSparseConfig::new(3, 2);
                    config.block_topk = 3;
                    config.minibatch = batch;
                    config.aux_k = 2;
                    let mut state =
                        BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
                    state.gamma = gamma;
                    state.pending_birth = Some(super::PendingBlockBirth {
                        block: 0,
                        baseline_decoder: decoder.clone(),
                        baseline_gamma,
                        baseline_rss: 0.0,
                        baseline_rows: 0,
                        baseline_usage: vec![0; 3],
                        baseline_second: (0..3).map(|_| Array2::zeros((2, 2))).collect(),
                    });
                    for shard in x.axis_chunks_iter(ndarray::Axis(0), shard_rows) {
                        state.partial_fit(shard).unwrap();
                    }
                    assert_eq!(state.row_count, x.nrows());
                    assert_eq!(state.usage, vec![16; 3]);
                    assert_eq!(state.alive_count, 3);
                    assert!((state.rss - expected_rss).abs() < 1e-11);
                    let num = x
                        .iter()
                        .zip(total.iter())
                        .map(|(&x, &v)| x as f64 * v)
                        .sum::<f64>();
                    let den = total.iter().map(|v| v * v).sum::<f64>();
                    assert!((state.gamma_num - num).abs() < 1e-11);
                    assert!((state.gamma_den - den).abs() < 1e-11);
                    let pending = state.pending_birth.as_ref().unwrap();
                    assert_eq!(pending.baseline_rows, x.nrows());
                    assert_eq!(pending.baseline_usage, state.usage);
                    assert!((pending.baseline_rss - expected_baseline_rss).abs() < 1e-11);
                    for block in 0..3 {
                        let w = weights.slice(ndarray::s![.., block * 2..(block + 1) * 2]);
                        let own = w.dot(
                            &decoder
                                .slice(ndarray::s![block * 2..(block + 1) * 2, ..])
                                .mapv(f64::from),
                        );
                        let v = &own * 3.0 - &total;
                        let x64 = x.mapv(f64::from);
                        let frame = decoder
                            .slice(ndarray::s![block * 2..(block + 1) * 2, ..])
                            .mapv(f64::from);
                        let expected_coupling = (x64.t().dot(&v) + v.t().dot(&x64)).dot(&frame.t());
                        let expected_energy = x64.iter().map(|v| v * v).sum::<f64>();
                        let expected_bound = x64
                            .outer_iter()
                            .zip(v.outer_iter())
                            .map(|(x, v)| {
                                (x.dot(&x).sqrt() * v.dot(&v).sqrt() - x.dot(&v)).max(0.0)
                            })
                            .sum::<f64>();
                        assert!((state.data_energy[block] - expected_energy).abs() < 1e-11);
                        assert!((state.negative_bound[block] - expected_bound).abs() < 1e-11);
                        let expected_data = x.mapv(f64::from).t().dot(&w);
                        let expected_second = w.t().dot(&w);
                        for (got, expected) in [
                            (&state.coupling[block], expected_coupling),
                            (&state.data_cross[block], expected_data),
                            (&state.second[block], expected_second.clone()),
                            (
                                &pending.baseline_second[block],
                                expected_second * (baseline_gamma as f64).powi(2),
                            ),
                        ] {
                            assert!(
                                got.iter()
                                    .zip(expected.iter())
                                    .all(|(a, b)| (a - b).abs() < 1e-11)
                            );
                        }
                    }
                    for feature in 0..4 {
                        let column = x.column(feature);
                        assert!(
                            (state.col_sum[feature]
                                - column.iter().map(|&v| v as f64).sum::<f64>())
                            .abs()
                                < 1e-11
                        );
                        assert!(
                            (state.col_sumsq[feature]
                                - column.iter().map(|&v| (v as f64).powi(2)).sum::<f64>())
                            .abs()
                                < 1e-11
                        );
                    }
                    let mut worst: Vec<usize> = (0..x.nrows()).collect();
                    worst.sort_by(|&a, &b| {
                        let norm = |row| residual.row(row).iter().map(|v| v * v).sum::<f64>();
                        norm(b).total_cmp(&norm(a)).then(a.cmp(&b))
                    });
                    let ranked = state.reservoir.ranked();
                    assert_eq!(ranked.len(), 4);
                    for (entry, &row) in ranked.iter().zip(&worst) {
                        assert_eq!(entry.global_index, row as u64);
                        assert!(
                            entry
                                .residual
                                .iter()
                                .zip(residual.row(row).iter())
                                .all(|(&a, &b)| (a as f64 - b).abs() < 1e-6)
                        );
                    }
                }
            }
        });
    }
}

#[test]
fn frame_refresh_uses_the_new_gamma_and_descends_from_any_initial_scale_2825() {
    let (x, decoder, config) = coupled_fixture();
    for orientation in [-1.0, 1.0] {
        let decoder = decoder.mapv(|value| orientation * value);
        for gamma in [0.2, 1.0, 2.0] {
            let mut state =
                BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
            state.gamma = gamma;
            for row in x.outer_iter() {
                state
                    .partial_fit(row.insert_axis(ndarray::Axis(0)))
                    .unwrap();
            }
            let stats = state.end_epoch().unwrap();
            assert!(!stats.converged);
            // #2825: the epoch's frames are NO LONGER invariant to the scale it starts
            // from, and that is the corrected behaviour rather than a regression. The
            // support step admits a block only when it lowers `‖x − γ Σ P_g x‖²`, an
            // objective that carries γ; the old invariance held because the top-k gate
            // rule ignored γ entirely, which is the same blindness that let the support
            // step raise the objective the frame and γ steps lower. The scale-free
            // alternative was measured and does NOT close the `K ≫ rank` certificates
            // this rule closes, so the invariance is what gives way.
            //
            // What every seed must still do is what this test is named for: refresh the
            // scale and then DESCEND AT THE REFRESHED SCALE, for either gauge
            // orientation. Both are asserted below, per seed.
            // Recompute the tied codes at each candidate, independently of the
            // stream's moments, holding this epoch's support fixed.
            let loss = |candidate: &Array2<f32>| -> f64 {
                x.outer_iter()
                    .map(|row| {
                        let weights: Vec<f64> = candidate
                            .outer_iter()
                            .map(|direction| {
                                direction
                                    .iter()
                                    .zip(row.iter())
                                    .map(|(&d, &x)| d as f64 * x as f64)
                                    .sum()
                            })
                            .collect();
                        (0..x.ncols())
                            .map(|feature| {
                                let reconstruction: f64 = (0..2)
                                    .map(|block| {
                                        stats.gamma as f64
                                            * weights[block]
                                            * candidate[[block, feature]] as f64
                                    })
                                    .sum();
                                (row[feature] as f64 - reconstruction).powi(2)
                            })
                            .sum::<f64>()
                    })
                    .sum()
            };
            assert!(loss(&state.decoder) < loss(&decoder));
        }
    }
}

#[test]
fn tied_frame_update_descends_on_the_frozen_code_ascent_witness_2825() {
    let x = array![[0.4_f32, 4.0], [0.3, -3.0], [-0.3, 3.0]];
    let mut decoder = array![[1.0_f32, -10.0], [10.0, 3.0]];
    for mut row in decoder.outer_iter_mut() {
        let norm = row.iter().map(|&v| (v as f64).powi(2)).sum::<f64>().sqrt();
        row.mapv_inplace(|v| (v as f64 / norm) as f32);
    }
    let (_, _, config) = coupled_fixture();
    let mut state = BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
    let x64 = x.mapv(f64::from);
    let loss = |frame: &Array2<f32>, gamma: f32| {
        let d = frame.mapv(f64::from);
        let residual = &x64 - x64.dot(&d.t()).dot(&d) * gamma as f64;
        residual.iter().map(|v| v * v).sum::<f64>()
    };
    state.partial_fit(x.view()).unwrap();
    let first = state.end_epoch().unwrap();
    let proposal = state.decoder.clone();
    let before = loss(&decoder, first.gamma);
    let after = loss(&proposal, first.gamma);
    assert!(
        after < before,
        "actual tied RSS ascended: {before} -> {after}"
    );
    state.partial_fit(x.view()).unwrap();
    let second = state.end_epoch().unwrap();
    let profiled_after = loss(&proposal, second.gamma);
    assert!(profiled_after <= after + 1e-12);
    assert!(profiled_after < before);
}

#[test]
fn tied_projector_moments_match_actual_loss_directional_derivatives() {
    // Exact orthonormal axes make horizontal perturbations independent of
    // f32 frame roundoff. The blocks overlap and have rank two.
    let decoder = array![
        [1.0_f32, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0],
    ];
    let x = Array2::from_shape_fn((11, 4), |(i, j)| ((i * 7 + j * 3) as f32).sin());
    let x64 = x.mapv(f64::from);
    let mut config = BlockSparseConfig::new(2, 2);
    config.block_topk = 2;
    config.aux_k = 0;
    for gamma in [0.3, 1.4] {
        let mut state = BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
        state.partial_fit(x.view()).unwrap();
        let data_scale = 2.0 * gamma - 2.0 * gamma * gamma;
        for block in 0..2 {
            let d = decoder.mapv(f64::from);
            let frame = d.slice(ndarray::s![block * 2..(block + 1) * 2, ..]);
            let raw = Array2::from_shape_fn((2, 4), |(i, j)| ((i * 5 + j + 1) as f64).cos());
            let tangent = &raw - raw.dot(&frame.t()).dot(&frame);
            let moment =
                &state.data_cross[block] * data_scale + &state.coupling[block] * (gamma * gamma);
            let analytic = -2.0
                * moment
                    .iter()
                    .zip(tangent.t().iter())
                    .map(|(a, b)| a * b)
                    .sum::<f64>();
            // The streamed moment is the derivative of the loss AT FIXED SUPPORT, so
            // the numerical reference has to hold the same support. `x D' D` sums every
            // block for every row, which was the same thing only while the router took
            // its top `k` unconditionally; a row now declines a block that does not
            // lower its loss (#2825), and admission boundaries are kinks this central
            // difference would otherwise straddle. Take the support ONCE from the
            // unperturbed frames, then re-project it at each candidate.
            // Route at the scale the STREAM routed at, not the scale this loop prices
            // the moments with: the accumulators were built under `state.gamma`, and a
            // reference routed at a different scale would fix a different support.
            let (blocks, gates, _) = crate::sparse_dict::block_sparse_dictionary_transform(
                x.view(),
                decoder.view(),
                state.gamma,
                config.block_size,
                config.block_topk,
                config.block_tile,
            )
            .expect("route the fixed support from the unperturbed frames");
            let loss = |step: f64| {
                let mut candidate = d.clone();
                let perturbed = &frame.to_owned() + &tangent * step;
                candidate
                    .slice_mut(ndarray::s![block * 2..(block + 1) * 2, ..])
                    .assign(&perturbed);
                (0..x64.nrows())
                    .map(|row| {
                        let mut reconstruction = vec![0.0f64; x64.ncols()];
                        for slot in 0..blocks.ncols() {
                            if gates[[row, slot]] == 0.0 {
                                continue;
                            }
                            let selected = blocks[[row, slot]] as usize;
                            for axis in 0..config.block_size {
                                let direction = candidate.row(selected * config.block_size + axis);
                                let weight: f64 = direction
                                    .iter()
                                    .zip(x64.row(row).iter())
                                    .map(|(&u, &value)| u * value)
                                    .sum();
                                for (accumulated, &u) in
                                    reconstruction.iter_mut().zip(direction.iter())
                                {
                                    *accumulated += gamma * weight * u;
                                }
                            }
                        }
                        x64.row(row)
                            .iter()
                            .zip(reconstruction.iter())
                            .map(|(value, fitted)| (value - fitted).powi(2))
                            .sum::<f64>()
                    })
                    .sum::<f64>()
            };
            let step = 1e-5;
            let numerical = (loss(step) - loss(-step)) / (2.0 * step);
            assert!(
                (analytic - numerical).abs() < 1e-7 * (1.0 + analytic.abs()),
                "block {block}, gamma {gamma}: analytic {analytic}, numerical {numerical}"
            );
        }
    }
}

#[test]
fn parallel_stream_rejects_selected_duplicate_birth_using_complete_baseline() {
    let x = Array2::from_shape_fn(
        (17, 2),
        |(_, column)| if column == 0 { 1.0_f32 } else { 0.0 },
    );
    let baseline = array![[0.0_f32, 0.0], [1.0, 0.0]];
    let candidate = array![[1.0_f32, 0.0], [1.0, 0.0]];
    let mut config = BlockSparseConfig::new(2, 1);
    config.block_topk = 1;
    config.minibatch = 5;
    config.aux_k = 1;
    config.frame_ridge = 0.0;
    let mut state = BlockSparseStreamState::new_with_decoder(candidate, &config).unwrap();
    state.pending_birth = Some(super::PendingBlockBirth {
        block: 0,
        baseline_decoder: baseline.clone(),
        baseline_gamma: 1.0,
        baseline_rss: 0.0,
        baseline_rows: 0,
        baseline_usage: vec![0; 2],
        baseline_second: (0..2).map(|_| Array2::zeros((1, 1))).collect(),
    });
    state.partial_fit(x.view()).unwrap();
    assert_eq!(
        state.usage,
        vec![17, 0],
        "the duplicate must win the routing tie"
    );
    let pending = state.pending_birth.as_ref().unwrap();
    assert_eq!(pending.baseline_usage, vec![0, 17]);
    assert_eq!(pending.baseline_rss, state.rss);
    let stats = state.end_epoch().unwrap();
    assert_eq!(stats.accepted_births, 0);
    assert!(
        !stats.converged,
        "rejection requires a measured baseline pass"
    );
    assert_eq!(state.decoder, baseline);
    assert_eq!(state.last_usage, vec![0, 17]);
    assert_eq!(state.last_second[1][[0, 0]], 17.0);
}

#[test]
fn a_large_spectral_shift_cannot_certify_a_nonstationary_frame() {
    let (x, decoder, mut config) = coupled_fixture();
    // Force a proposal below f32 resolution while leaving the actual tied
    // objective unchanged. The gradient certificate must see through it.
    config.frame_ridge = 1e20;
    let mut first = BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
    first.partial_fit(x.view()).unwrap();
    let measured = first.end_epoch().unwrap();
    let mut state = BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
    state.gamma = measured.gamma;
    state.prev_ev = measured.explained_variance;
    state.epochs_run = 1;
    state.partial_fit(x.view()).unwrap();
    let stats = state.end_epoch().unwrap();
    let displacement = crate::sparse_dict::block::frame_fixed_point_residual(
        decoder.view(),
        state.decoder.view(),
        2,
        1,
    )
    .unwrap();
    assert!(displacement <= config.tolerance);
    assert!(stats.gamma_residual <= config.tolerance);
    assert!(stats.frame_residual > config.tolerance);
    assert!(!stats.converged);
    assert!(state.finalize().is_err());
}

#[test]
fn equal_ev_cannot_certify_changing_gamma_or_frames_2825() {
    let (x, decoder, config) = coupled_fixture();
    let mut first = BlockSparseStreamState::new_with_decoder(decoder.clone(), &config).unwrap();
    first.partial_fit(x.view()).unwrap();
    let measured = first.end_epoch().unwrap();
    let mut state = BlockSparseStreamState::new_with_decoder(decoder, &config).unwrap();
    state.prev_ev = measured.explained_variance;
    state.epochs_run = 1;
    state.partial_fit(x.view()).unwrap();
    let stats = state.end_epoch().unwrap();
    assert_eq!(stats.explained_variance, measured.explained_variance);
    // The contract this control exists for is the one its NAME states: an equal EV
    // cannot certify while gamma OR the frames are still moving. Requiring BOTH to be
    // open was incidental to a router that re-derived a different support every epoch;
    // with the support step descending the shared objective the scale settles first, so
    // the conjunction would pin the old churn rather than the contract. Name which
    // residual is holding the certificate open, so a future change that closes BOTH —
    // which WOULD make this fixture vacuous — fails here instead of passing silently.
    let gamma_open = stats.gamma_residual > config.tolerance;
    let frame_open = stats.frame_residual > config.tolerance;
    assert!(
        gamma_open || frame_open,
        "equal EV with both residuals closed is a converged fit, not this control: \
         gamma_residual={} frame_residual={} tol={}",
        stats.gamma_residual,
        stats.frame_residual,
        config.tolerance
    );
    assert!(!stats.converged);
    assert!(state.finalize().is_err());
}

#[test]
fn coupled_stream_certifies_the_returned_frames_gamma_and_fresh_ev_2825() {
    let (x, decoder, config) = coupled_fixture();
    let mut state = BlockSparseStreamState::new_with_decoder(decoder, &config).unwrap();
    let mut last = None;
    for _ in 0..config.max_epochs {
        state.partial_fit(x.view()).unwrap();
        let stats = state.end_epoch().unwrap();
        last = Some(stats);
        if stats.converged {
            break;
        }
    }
    let last = last.unwrap();
    assert!(last.converged, "{last:?}");
    assert!(last.gamma_residual <= config.tolerance);
    assert!(last.frame_residual <= config.tolerance);
    let artifact = state.finalize().unwrap();
    let mean = x.mean_axis(ndarray::Axis(0)).unwrap();
    let mut rss = 0.0;
    let mut tss = 0.0;
    for row in x.outer_iter() {
        let mut reconstructed = vec![0.0_f64; x.ncols()];
        for direction in artifact.decoder.outer_iter() {
            let projection: f64 = direction
                .iter()
                .zip(row.iter())
                .map(|(&d, &x)| d as f64 * x as f64)
                .sum();
            for (feature, value) in reconstructed.iter_mut().enumerate() {
                *value += artifact.gamma as f64 * projection * direction[feature] as f64;
            }
        }
        for feature in 0..x.ncols() {
            rss += (row[feature] as f64 - reconstructed[feature]).powi(2);
            tss += (row[feature] as f64 - mean[feature] as f64).powi(2);
        }
    }
    assert!((artifact.explained_variance - (1.0 - rss / tss)).abs() < 1e-6);
    assert!(artifact.explained_variance > 0.99);
    state.partial_fit(x.view()).unwrap();
    assert!(
        state.finalize().is_err(),
        "a new unclosed pass invalidates the certificate"
    );
}

#[test]
fn overcomplete_stream_accepts_one_evidence_birth_then_dead_tail_is_quiescent_2023() {
    // Rank-2 data with G=16 reproduces the K≫intrinsic-rank boundary behind
    // #2023. Block 0 starts on e0 and every other frame is dead. Exactly one e1
    // residual birth is warranted; after it commits the remaining fourteen dead
    // blocks must stay quiescent so the stream can certify instead of reseeding
    // them forever.
    let (rows, p, g, b) = (64usize, 2usize, 16usize, 1usize);
    let x = Array2::<f32>::from_shape_fn(
        (rows, p),
        |(row, column)| {
            if column == row % 2 { 1.0 } else { 0.0 }
        },
    );
    let mut decoder = Array2::<f32>::zeros((g * b, p));
    decoder[[0, 0]] = 1.0;
    let cfg = BlockSparseConfig {
        n_blocks: g,
        block_size: b,
        block_topk: 1,
        max_epochs: 8,
        minibatch: rows,
        block_tile: g,
        frame_ridge: 0.0,
        aux_k: g,
        matryoshka_prefix: false,
        tolerance: 0.0,
    };
    let mut state = BlockSparseStreamState::new_with_decoder(decoder, &cfg).expect("stream state");
    let mut accepted_total = 0usize;
    let mut saw_pending = false;
    let mut final_stats = None;
    for _ in 0..cfg.max_epochs {
        state.partial_fit(x.view()).expect("stream rank-2 corpus");
        let stats = state.end_epoch().expect("close rank-2 epoch");
        accepted_total += stats.accepted_births;
        saw_pending |= stats.birth_pending;
        final_stats = Some(stats);
        if stats.converged {
            break;
        }
    }
    let final_stats = final_stats.expect("at least one epoch");
    assert!(saw_pending, "a residual-row birth must be staged for e1");
    assert_eq!(
        accepted_total, 1,
        "only the missing rank-1 direction has positive exact evidence"
    );
    assert!(final_stats.converged, "dead tail prevented certification");
    assert!(!final_stats.birth_pending);
    assert_eq!(final_stats.dead, g - 2);

    let artifact = state.finalize().expect("quiescent overcomplete artifact");
    assert_eq!(
        artifact
            .block_utilization
            .iter()
            .filter(|&&value| value > 0.0)
            .count(),
        2,
    );
    assert!((artifact.explained_variance - 1.0).abs() <= f64::EPSILON);
}