gam-sae 0.3.153

//! #2027 — deterministic repro for the K≥2 whitened dictionary CO-COLLAPSE, and
//! the regression guard for the disjoint-subspace / ownership-anchor / reseed-
//! hysteresis fix.
//!
//! Two planted circles live in DISJOINT 2-planes of an ambient `p`-dim cloud; the
//! per-column-standardized ("whitened") target is their sum, so a faithful K=2
//! reconstruction REQUIRES both atoms to carry signal on different subspaces.
//! Before the fix the joint decoder refit at the co-collapse reseed re-spread one
//! residual direction across both atoms and the gate let them trade rows, so the
//! dictionary re-symmetrised into a single shared basin and the reconstruction EV
//! collapsed to the no-signal level. With the greedy disjoint-subspace decoder
//! refit + soft row-ownership anchor + reseed cooldown the two atoms hold distinct
//! territories and the fit recovers a materially positive EV.

use super::tests::deterministic_circle_noise;
use super::*;

/// Whitened two-circle target: circle A lives on the even ambient columns, circle
/// B on the odd ones (disjoint deterministic near-orthonormal 2-frames), driven by
/// two INCOMMENSURATE phases so the circles are not row-aligned. Each column is
/// standardized to zero mean / unit variance (the whitening proxy that puts both
/// circles on a common scale, the regime the real-data co-collapse lives in).
fn two_circle_whitened_target(n: usize, p: usize, sigma: f64) -> Array2<f64> {
    let mut fa = Array2::<f64>::zeros((2, p));
    let mut fb = Array2::<f64>::zeros((2, p));
    for j in 0..p {
        if j % 2 == 0 {
            fa[[0, j]] = deterministic_circle_noise(j, 0);
            fa[[1, j]] = deterministic_circle_noise(j, 1);
        } else {
            fb[[0, j]] = deterministic_circle_noise(j, 2);
            fb[[1, j]] = deterministic_circle_noise(j, 3);
        }
    }
    for f in [&mut fa, &mut fb] {
        for r in 0..2 {
            let nrm = (0..p).map(|j| f[[r, j]] * f[[r, j]]).sum::<f64>().sqrt();
            for j in 0..p {
                f[[r, j]] /= nrm.max(1.0e-300);
            }
        }
    }
    let mut z = Array2::<f64>::zeros((n, p));
    for row in 0..n {
        let ta = std::f64::consts::TAU * (row as f64) / (n as f64);
        let tb = std::f64::consts::TAU * (2.0 * row as f64 + 0.37) / (n as f64);
        let (ca, sa) = (ta.cos(), ta.sin());
        let (cb, sb) = (tb.cos(), tb.sin());
        for j in 0..p {
            z[[row, j]] = ca * fa[[0, j]]
                + sa * fa[[1, j]]
                + cb * fb[[0, j]]
                + sb * fb[[1, j]]
                + sigma * deterministic_circle_noise(row, j + 7);
        }
    }
    for j in 0..p {
        let mut mean = 0.0_f64;
        for row in 0..n {
            mean += z[[row, j]];
        }
        mean /= n as f64;
        let mut var = 0.0_f64;
        for row in 0..n {
            let d = z[[row, j]] - mean;
            var += d * d;
        }
        let sd = (var / n as f64).sqrt().max(1.0e-12);
        for row in 0..n {
            z[[row, j]] = (z[[row, j]] - mean) / sd;
        }
    }
    z
}

/// Build a fresh K=2 periodic term seeded (production PCA seed) from the whitened
/// two-circle target, decoders cold at zero.
fn two_circle_k2_term(n: usize, p: usize, m: usize) -> (SaeManifoldTerm, Array2<f64>) {
    let d = 1usize;
    let k = 2usize;
    let target = two_circle_whitened_target(n, p, 0.05);
    let basis_kinds = vec![SaeAtomBasisKind::Periodic; k];
    let dims = vec![d; k];
    let seed = sae_pca_seed_initial_coords(target.view(), &basis_kinds, &dims)
        .expect("the PCA seed covers every declared atom basis kind and dim");
    let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(m)
        .expect("a periodic harmonic evaluator exists for this odd basis count"));

    let mut basis_values = Array3::<f64>::zeros((k, n, m));
    let mut basis_jacobian = Array4::<f64>::zeros((k, n, m, d));
    let decoder = Array3::<f64>::zeros((k, m, p));
    let mut penalties = Array3::<f64>::zeros((k, m, m));
    let mut coords_vec: Vec<Array2<f64>> = Vec::new();
    for atom in 0..k {
        let coords = seed.slice(s![atom, .., 0..d]).to_owned();
        let (phi, jet) = evaluator.evaluate(coords.view())
            .expect("the periodic evaluator accepts the seeded coordinate block");
        basis_values.slice_mut(s![atom, .., ..]).assign(&phi);
        basis_jacobian.slice_mut(s![atom, .., .., ..]).assign(&jet);
        penalties
            .slice_mut(s![atom, .., ..])
            .assign(&Array2::<f64>::eye(m));
        coords_vec.push(coords);
    }
    let logits = Array2::<f64>::zeros((n, k));
    let mut evaluators: Vec<Option<Arc<dyn SaeBasisSecondJet>>> = Vec::new();
    for _ in 0..k {
        evaluators.push(Some(evaluator.clone()));
    }
    let term = term_from_padded_blocks_with_mode(
        n,
        p,
        &basis_kinds,
        basis_values.view(),
        basis_jacobian.view(),
        &vec![m; k],
        &dims,
        decoder.view(),
        penalties.view(),
        logits.view(),
        &coords_vec,
        AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        &evaluators,
    )
    .expect("the fixture assignment blocks match the declared mode");
    (term, target)
}

/// The K=2 whitened two-circle fit must recover a materially positive
/// reconstruction EV — the concrete bar this test asserts is `EV > 0.20`. Two disjoint
/// circles together span a rank-4 subspace of the whitened cloud, so an honest K=2
/// dictionary explains a large fraction of the variance (the torch proxy reaches
/// ≈0.47 on the sibling nursery experiment). The disjoint-subspace decoder refit,
/// row-ownership anchor, and reseed cooldown keep the two atoms on distinct
/// territories through the joint solve.
#[test]
pub(crate) fn two_circle_whitened_k2_recovers_disjoint_signal_2027() {
    let n = 96usize;
    let p = 16usize;
    let m = 5usize; // [1, sin2πt, cos2πt, sin4πt, cos4πt]
    let (mut term, target) = two_circle_k2_term(n, p, m);

    let mut rho = SaeManifoldRho::new(
        0.0,
        -6.0,
        vec![Array1::<f64>::zeros(1), Array1::<f64>::zeros(1)],
    );
    let loss = term
        .run_joint_fit_arrow_schur(target.view(), &mut rho, None, 60, 0.05, 1.0e-3, 1.0e-3)
        .expect("the joint arrow-Schur fit converges on this fixture");
    assert!(loss.total().is_finite(), "loss must stay finite");

    let ev = term
        .dictionary_reconstruction_ev(target.view(), &rho)
        .expect("reconstruction EV is defined for a fitted term");
    eprintln!(
        "[#2027 repro] K=2 whitened two-circle EV = {ev:.4}, cocollapse_reseeds = {}",
        term.dictionary_cocollapse_reseeds
    );
    assert!(
        ev > 0.20,
        "K=2 whitened two-circle dictionary co-collapsed: EV={ev:.4} (expected > 0.20; \
         two disjoint circles span a rank-4 subspace that an honest K=2 fit recovers)"
    );
}

/// The greedy disjoint-subspace decoder refit must, on a co-collapsed reseed,
/// leave BOTH atoms carrying material decoder norm — never let one atom take all
/// the residual while the other stays ≈0 (the relative-norm collapse the joint
/// refit permitted). A direct unit check of `refit_decoder_sequential_deflation`.
#[test]
pub(crate) fn sequential_deflation_gives_both_atoms_material_norm_2027() {
    let n = 96usize;
    let p = 16usize;
    let m = 5usize;
    let (mut term, target) = two_circle_k2_term(n, p, m);
    term.refit_decoder_sequential_deflation(target.view())
        .expect("sequential deflation refits the fixture decoders");
    let mut norms = [0.0_f64; 2];
    for (atom_idx, atom) in term.atoms.iter().enumerate() {
        norms[atom_idx] = atom
            .decoder_coefficients()
            .iter()
            .map(|v| v * v)
            .sum::<f64>()
            .sqrt();
    }
    let (lo, hi) = if norms[0] <= norms[1] {
        (norms[0], norms[1])
    } else {
        (norms[1], norms[0])
    };
    eprintln!("[#2027 repro] deflation decoder norms = {norms:?}");
    assert!(hi > 0.0, "at least one atom must carry decoder norm");
    assert!(
        lo > 1.0e-3 * hi,
        "both atoms must carry material decoder norm after deflation: norms={norms:?}"
    );
}

/// #2027 WIDTH-SCALING + STRUCTURE-RECOVERY guard — the discriminating test.
///
/// Two facts from the sibling nursery evidence make raw EV an INSUFFICIENT guard:
///   1. The pathology is WIDTH-dependent — the REML control converges at `p = 16`
///      but hangs / co-collapses at `p ≈ 96`. A fix must be checked at BOTH widths:
///      the narrow arm must stay healthy, the wide arm is the regime being rescued.
///   2. Its fingerprint is a co-collapse that posts a DECENT reconstruction EV while
///      recovering NEITHER planted circle — both atoms pile into one shared subspace
///      and ride one circle plus noise (torch proxy: EV 0.63, adjacency 0.43/0.25).
///      Asserting EV alone therefore passes a co-collapsed fit.
///
/// The two circles are planted on DISJOINT ambient column PARITIES (circle A on the
/// even output channels, circle B on the odd), so an honest K=2 dictionary MUST
/// separate: one atom's decoder concentrates its Frobenius energy on the even
/// channels, the other on the odd. Co-collapse piles both atoms onto the same
/// channels — detected here as both atoms landing on the SAME side of the 0.5
/// even-energy split. We require, at both widths: finite loss (no thrash), a
/// materially positive EV, and the two atoms SEPARATED onto opposite-parity
/// subspaces (the structure the disjoint-deflation + ownership-anchor fix restores).
///
/// NOTE: this exercises the INNER joint solve (`run_joint_fit_arrow_schur` at a
/// fixed ρ) — the co-collapse / structure-recovery layer the seeding + anchoring fix
/// lives in — not the outer penalized quasi-Laplace ρ-search whose non-PD-Hessian retries are the
/// separate Python-side "hang" at wide `p`.
#[test]
pub(crate) fn two_circle_separates_at_narrow_and_wide_widths_2027() {
    let m = 5usize;
    // (n, p): a NARROW arm that must stay healthy and the WIDE arm being rescued.
    for &(n, p) in &[(96usize, 16usize), (120usize, 96usize)] {
        let (mut term, target) = two_circle_k2_term(n, p, m);
        let mut rho = SaeManifoldRho::new(
            0.0,
            -6.0,
            vec![Array1::<f64>::zeros(1), Array1::<f64>::zeros(1)],
        );
        let loss = term
            .run_joint_fit_arrow_schur(target.view(), &mut rho, None, 60, 0.05, 1.0e-3, 1.0e-3)
            .expect("the joint arrow-Schur fit converges on this fixture");
        assert!(
            loss.total().is_finite(),
            "p={p}: joint fit must return a finite loss (no thrash / NaN)"
        );
        let ev = term
            .dictionary_reconstruction_ev(target.view(), &rho)
            .expect("reconstruction EV is defined for a fitted term");

        // Per-atom decoder energy split across even vs odd output channels: circle A
        // lives on even channels, circle B on odd.
        let mut even_frac = [0.0_f64; 2];
        for (atom_idx, atom) in term.atoms.iter().enumerate() {
            let b = atom.decoder_coefficients(); // (m × p)
            let mut e_even = 0.0_f64;
            let mut e_odd = 0.0_f64;
            for col in 0..b.nrows() {
                for out in 0..p {
                    let v = b[[col, out]] * b[[col, out]];
                    if out % 2 == 0 {
                        e_even += v;
                    } else {
                        e_odd += v;
                    }
                }
            }
            even_frac[atom_idx] = e_even / (e_even + e_odd).max(1.0e-300);
        }
        eprintln!(
            "[#2027 repro] p={p}: EV={ev:.4}, per-atom even-energy fraction={even_frac:?}, \
             cocollapse_reseeds={}",
            term.dictionary_cocollapse_reseeds
        );
        assert!(
            ev > 0.20,
            "p={p}: dictionary co-collapsed to the null floor (EV={ev:.4} <= 0.20)"
        );
        // STRUCTURE RECOVERY: the atoms must land on OPPOSITE planted subspaces — one
        // even-dominant, one odd-dominant. Both on the same side of 0.5 is the
        // co-collapse signature (acceptable EV, neither circle recovered).
        let (lo, hi) = if even_frac[0] <= even_frac[1] {
            (even_frac[0], even_frac[1])
        } else {
            (even_frac[1], even_frac[0])
        };
        assert!(
            lo < 0.5 && hi > 0.5,
            "p={p}: atoms did NOT separate onto the two planted circles \
             (even-energy fractions {even_frac:?} both on one side of 0.5 = co-collapse: \
             EV looks fine but neither circle is recovered)"
        );
    }
}

/// #2082/#2132/#1893 — the STRUCTURAL coherence detector fires on FUNCTIONAL
/// REDUNDANCY (two atoms that reconstruct the SAME rows — a genuine duplicate),
/// NOT on mere output-subspace sharing. Three cases pin the contract:
///
///  (1) ORTHOGONAL output subspaces → frames don't even overlap → NOT flagged.
///  (2) SAME output subspace but DIFFERENT charts (identical decoder, distinct
///      phases) → the atoms decode DIFFERENT rows, so their gated contributions
///      `Y_k = diag(a)ΦB` are NOT collinear → benign, NOT flagged. This is the
///      over-complete (`K > rank`) regime the old frame-coherence detector
///      false-positived on (the `ordered_beta_bernoulli_default_alpha` regression: healthy EV≈0.99,
///      frame coherence ≈1, contribution cosine ≈ the independence null): several
///      curved atoms MUST share the ≤`p`-dim output space while encoding distinct
///      structure.
///  (3) TRUE DUPLICATE (identical decoder AND identical chart) → `Y_0 ∝ Y_1` →
///      contribution cosine → 1 → FLAGGED.
#[test]
pub(crate) fn structural_coherence_detector_fires_on_duplicate_not_orthogonal_2082() {
    let n = 48usize;
    let p = 8usize;
    let m = 5usize;

    // (1) ORTHOGONAL output subspaces: atom 0 decodes only EVEN output channels,
    // atom 1 only ODD → orthogonal frames → not a candidate → NOT flagged.
    let (mut term, _target) = two_circle_k2_term(n, p, m);
    for atom in 0..2 {
        let mut b = Array2::<f64>::zeros((m, p));
        for col in 0..m {
            let out = (if atom == 0 { 0 } else { 1 }) + 2 * (col % (p / 2));
            if out < p {
                b[[col, out]] = 1.0;
            }
        }
        term.atoms[atom].set_decoder_coefficients(b).expect("decoder matches its atom basis");
    }
    assert!(
        term.structural_coherence_collapse_detected()
            .expect("collapse detection is defined for a fully built term")
            .is_none(),
        "orthogonal-subspace atoms must NOT be flagged as structurally collapsed"
    );

    // (2) SAME output subspace, DIFFERENT charts: identical decoder, but the two
    // atoms keep their distinct PCA-seeded phases → same output frame (coherence
    // ≈1) yet DIFFERENT per-row contributions → NOT functional redundancy → the
    // functional-redundancy detector must stay SILENT (the old frame-only detector
    // wrongly fired here).
    let mut dup = Array2::<f64>::zeros((m, p));
    dup[[1, 0]] = 1.0;
    dup[[2, 1]] = 1.0;
    term.atoms[0].set_decoder_coefficients(dup.clone()).expect("decoder matches its atom basis");
    term.atoms[1].set_decoder_coefficients(dup.clone()).expect("decoder matches its atom basis");
    let mut shifted = term.assignment.coords[0].as_matrix().to_owned();
    for t in shifted.iter_mut() {
        *t = (*t + 0.25).rem_euclid(1.0);
    }
    let shifted_flat: Array1<f64> = shifted.iter().copied().collect();
    term.assignment.coords[1].set_flat(shifted_flat.view());
    term.atoms[1].refresh_basis(shifted.view())
        .expect("the atom basis refreshes at the supplied coords");
    assert!(
        term.structural_coherence_collapse_detected()
            .expect("collapse detection is defined for a fully built term")
            .is_none(),
        "same output subspace with DIFFERENT charts is benign over-completeness and \
         must NOT be flagged (the ordered_beta_bernoulli_default_alpha false positive)"
    );

    // (3) TRUE DUPLICATE: identical decoder AND identical chart (copy atom 0's
    // coords onto atom 1) → the two atoms reconstruct the SAME rows → contribution
    // cosine → 1 → FLAGGED as the genuine high-EV co-collapse.
    let coords0 = term.assignment.coords[0].as_matrix().to_owned();
    let flat0: Array1<f64> = coords0.iter().copied().collect();
    term.assignment.coords[1].set_flat(flat0.view());
    term.atoms[1]
        .refresh_basis(term.assignment.coords[1].as_matrix().view())
        .expect("the atom basis refreshes at the supplied coords");
    let hit = term
        .structural_coherence_collapse_detected()
        .expect("collapse detection is defined for a fully built term")
        .expect("a true duplicate (identical decoder AND chart) must be flagged");
    assert_eq!((hit.0, hit.1), (0, 1), "the offending pair is (0, 1)");
    assert!(
        hit.2 > 0.9,
        "true-duplicate contribution cosine must be ~1, got {}",
        hit.2
    );
}

/// #2132 #2b — build a K=3 periodic term whose three atoms ALL decode into the
/// SAME 2-D output plane (the first-harmonic sin/cos map onto output columns 0
/// and 1 of a `p`-dim output), so the union output-frame rank `R = 2 < K = 3`
/// and the dictionary is OVERCOMPLETE. `duplicate = true` makes atoms 0 and 1 a
/// TRUE duplicate (identical decoder AND identical chart/phase) with atom 2 on a
/// distinct phase; `duplicate = false` gives all three DISTINCT phases (identical
/// decoder, different charts) — benign pigeonhole sharing.
fn overcomplete_k3_planar_term(n: usize, p: usize, m: usize, duplicate: bool) -> SaeManifoldTerm {
    let d = 1usize;
    let k = 3usize;
    let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(m)
        .expect("a periodic harmonic evaluator exists for this odd basis count"));
    let mut basis_values = Array3::<f64>::zeros((k, n, m));
    let mut basis_jacobian = Array4::<f64>::zeros((k, n, m, d));
    let mut decoder = Array3::<f64>::zeros((k, m, p));
    let mut penalties = Array3::<f64>::zeros((k, m, m));
    let mut coords_vec: Vec<Array2<f64>> = Vec::new();
    for atom in 0..k {
        // atoms 0,1 share phase 0 when duplicating (atom 2 shifted); otherwise all
        // three phases are distinct.
        let phase = if duplicate {
            if atom == 2 { 1.0 / 3.0 } else { 0.0 }
        } else {
            atom as f64 / k as f64
        };
        let mut coords = Array2::<f64>::zeros((n, d));
        for row in 0..n {
            coords[[row, 0]] = ((row as f64) / (n as f64) + phase).rem_euclid(1.0);
        }
        let (phi, jet) = evaluator.evaluate(coords.view())
            .expect("the periodic evaluator accepts the seeded coordinate block");
        basis_values.slice_mut(s![atom, .., ..]).assign(&phi);
        basis_jacobian.slice_mut(s![atom, .., .., ..]).assign(&jet);
        penalties
            .slice_mut(s![atom, .., ..])
            .assign(&Array2::<f64>::eye(m));
        // Every atom decodes the first harmonic into the (e0, e1) output plane
        // (sin → col 0, cos → col 1), so all three output frames span the same
        // 2-D subspace ⇒ union rank 2 < k = 3 (overcomplete).
        decoder[[atom, 1, 0]] = 1.0;
        decoder[[atom, 2, 1]] = 1.0;
        coords_vec.push(coords);
    }
    let logits = Array2::<f64>::zeros((n, k));
    let mut evaluators: Vec<Option<Arc<dyn SaeBasisSecondJet>>> = Vec::new();
    for _ in 0..k {
        evaluators.push(Some(evaluator.clone()));
    }
    term_from_padded_blocks_with_mode(
        n,
        p,
        &vec![SaeAtomBasisKind::Periodic; k],
        basis_values.view(),
        basis_jacobian.view(),
        &vec![m; k],
        &vec![d; k],
        decoder.view(),
        penalties.view(),
        logits.view(),
        &coords_vec,
        AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        &evaluators,
    )
    .expect("the fixture assignment blocks match the declared mode")
}

/// #2132 #2b — an OVERCOMPLETE (`K > R`) true duplicate must be detected. The old
/// union-frame-rank gate returned NO collapsed pairs whenever `K > R`, which
/// DISABLED the detector exactly in the overcomplete regime where duplicates are
/// most likely (the second-stage contribution-cosine verdict was never reached).
/// #2b keeps the detector alive overcomplete — it drops only the pigeonhole-
/// forced PASS-1 frame prune — so PASS 2 flags the true duplicate (0,1). This
/// fixture is genuinely overcomplete (three atoms in a 2-D output plane, R=2<K=3).
#[test]
pub(crate) fn overcomplete_duplicate_is_detected_past_the_k_gt_r_gate_2132() {
    let term = overcomplete_k3_planar_term(96, 4, 5, true);
    let hit = term
        .structural_coherence_collapse_detected()
        .expect("collapse detection is defined for a fully built term")
        .expect(
            "an overcomplete (K>R) true duplicate must be flagged now that the detector \
             reaches PASS 2 past the K>R gate",
        );
    assert_eq!((hit.0, hit.1), (0, 1), "the duplicate pair is (0, 1)");
    assert!(
        hit.2 > 0.9,
        "overcomplete duplicate contribution cosine must be ~1, got {}",
        hit.2
    );
}

/// #2132 #2b — the overcomplete detector must NOT false-fire on benign pigeonhole
/// sharing. Same K=3-in-a-2-plane geometry but all three atoms carry DISTINCT
/// phases (identical decoder, different charts): every output frame overlaps by
/// pigeonhole, yet the PASS-2 contribution cosine sits at the independence null,
/// so nothing is flagged (the `ordered_beta_bernoulli` false-positive the K>R gate was
/// originally meant to avoid — now avoided by the verdict stage, not by disabling
/// the detector).
#[test]
pub(crate) fn overcomplete_distinct_phases_stay_silent_2132() {
    let term = overcomplete_k3_planar_term(96, 4, 5, false);
    assert!(
        term.structural_coherence_collapse_detected()
            .expect("collapse detection is defined for a fully built term")
            .is_none(),
        "benign overcomplete pigeonhole sharing (distinct phases) must NOT be flagged"
    );
}

/// Two circles of UNEQUAL amplitude on disjoint output-channel parities: circle A
/// (unit amplitude, winding 1) on the even channels {0, 2}, circle B (`amp_b < 1`)
/// on the odd channels {1, 3}. Circle B winds THREE times per revolution of A, a
/// harmonic index BEYOND the atoms' order-2 (`m = 5`) chart span, so neither
/// circle lies in the other's harmonic reach: atom A's decoder cannot absorb B as
/// one of its own harmonics (which an absorbable 2× winding would let it do,
/// collapsing the residual to zero and making the second reseed correctly
/// terminal). A DOMINATES, so a co-collapse reseed that seeds both atoms from the
/// same residual reads circle A for BOTH (re-collision); only a sequential-
/// deflation reseed — peel A onto atom 0, then seed atom 1 from what A genuinely
/// leaves behind (circle B) — separates them onto the two disjoint circles.
fn two_amplitude_circle_target(n: usize, amp_b: f64) -> Array2<f64> {
    let p = 4usize;
    let mut z = Array2::<f64>::zeros((n, p));
    for row in 0..n {
        let ta = std::f64::consts::TAU * (row as f64) / (n as f64);
        let tb = std::f64::consts::TAU * (3.0 * row as f64 + 0.37) / (n as f64);
        z[[row, 0]] = ta.cos();
        z[[row, 2]] = ta.sin();
        z[[row, 1]] = amp_b * tb.cos();
        z[[row, 3]] = amp_b * tb.sin();
    }
    z
}

/// Build a fresh K=2 periodic term (production PCA seed, decoders cold at zero)
/// from an arbitrary target — the general form of [`two_circle_k2_term`].
fn k2_periodic_term_from_target(target: &Array2<f64>, m: usize) -> SaeManifoldTerm {
    let n = target.nrows();
    let p = target.ncols();
    let d = 1usize;
    let k = 2usize;
    let basis_kinds = vec![SaeAtomBasisKind::Periodic; k];
    let dims = vec![d; k];
    let seed = sae_pca_seed_initial_coords(target.view(), &basis_kinds, &dims)
        .expect("the PCA seed covers every declared atom basis kind and dim");
    let evaluator = Arc::new(PeriodicHarmonicEvaluator::new(m)
        .expect("a periodic harmonic evaluator exists for this odd basis count"));

    let mut basis_values = Array3::<f64>::zeros((k, n, m));
    let mut basis_jacobian = Array4::<f64>::zeros((k, n, m, d));
    let decoder = Array3::<f64>::zeros((k, m, p));
    let mut penalties = Array3::<f64>::zeros((k, m, m));
    let mut coords_vec: Vec<Array2<f64>> = Vec::new();
    for atom in 0..k {
        let coords = seed.slice(s![atom, .., 0..d]).to_owned();
        let (phi, jet) = evaluator.evaluate(coords.view())
            .expect("the periodic evaluator accepts the seeded coordinate block");
        basis_values.slice_mut(s![atom, .., ..]).assign(&phi);
        basis_jacobian.slice_mut(s![atom, .., .., ..]).assign(&jet);
        penalties
            .slice_mut(s![atom, .., ..])
            .assign(&Array2::<f64>::eye(m));
        coords_vec.push(coords);
    }
    let logits = Array2::<f64>::zeros((n, k));
    let mut evaluators: Vec<Option<Arc<dyn SaeBasisSecondJet>>> = Vec::new();
    for _ in 0..k {
        evaluators.push(Some(evaluator.clone()));
    }
    term_from_padded_blocks_with_mode(
        n,
        p,
        &basis_kinds,
        basis_values.view(),
        basis_jacobian.view(),
        &vec![m; k],
        &dims,
        decoder.view(),
        penalties.view(),
        logits.view(),
        &coords_vec,
        AssignmentMode::ordered_beta_bernoulli(1.0, 1.0, false),
        &evaluators,
    )
    .expect("the fixture assignment blocks match the declared mode")
}

/// #2132 births — the SEQUENTIAL-DEFLATION birth reseed must separate co-collapsed
/// CURVED atoms onto DISJOINT structure. With cold decoders the reconstruction
/// residual is the full two-circle target; reseeding both atoms must peel the
/// dominant circle A onto atom 0 and land atom 1 on circle B (what atom 0 left
/// behind), so their provisional decoders concentrate on OPPOSITE output-channel
/// parities. A simultaneous seed reads circle A for both atoms (both even-dominant
/// = re-collision), which this opposite-parity assertion rejects.
#[test]
pub(crate) fn birth_reseed_sequential_deflation_separates_curved_atoms_2132() {
    let n = 96usize;
    let m = 5usize;
    let target = two_amplitude_circle_target(n, 0.4);
    let p = target.ncols();
    let mut term = k2_periodic_term_from_target(&target, m);
    let rho = SaeManifoldRho::new(
        0.0,
        -6.0,
        vec![Array1::<f64>::zeros(1), Array1::<f64>::zeros(1)],
    );
    // Cold decoders ⇒ reconstruction residual == target (both circles uncovered),
    // so both atoms carry distinct signal and the sequence reseeds both.
    let reseeded = term
        .reseed_curved_atoms_sequential_deflation(&[0, 1], target.view(), &rho)
        .expect("reseeding is defined for the cold-decoder fixture");
    assert_eq!(
        reseeded,
        vec![0, 1],
        "both curved atoms sit above the residual noise floor, so both reseed"
    );

    // Per-atom provisional-decoder energy split across even vs odd output channels:
    // circle A lives on the even channels, circle B on the odd.
    let mut even_frac = [0.0_f64; 2];
    for (atom_idx, atom) in term.atoms.iter().enumerate() {
        let b = atom.decoder_coefficients(); // (m × p)
        let (mut e_even, mut e_odd) = (0.0_f64, 0.0_f64);
        for col in 0..b.nrows() {
            for out in 0..p {
                let v = b[[col, out]] * b[[col, out]];
                if out % 2 == 0 {
                    e_even += v;
                } else {
                    e_odd += v;
                }
            }
        }
        even_frac[atom_idx] = e_even / (e_even + e_odd).max(1.0e-300);
    }
    eprintln!("[#2132 birth reseed] per-atom even-energy fraction = {even_frac:?}");
    let (lo, hi) = if even_frac[0] <= even_frac[1] {
        (even_frac[0], even_frac[1])
    } else {
        (even_frac[1], even_frac[0])
    };
    assert!(
        lo < 0.5 && hi > 0.5,
        "sequential-deflation birth reseed did NOT separate the two curved atoms onto \
         disjoint structure: even-energy fractions {even_frac:?} both on one side of 0.5 \
         (a simultaneous seed re-reads the dominant circle for both atoms = re-collision)"
    );
}