gam-sae 0.3.127

Sparse-autoencoder latent-manifold terms for the gam penalized-likelihood engine
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
//! Per-point sparse atom codes for multi-manifold reconstruction.
//!
//! This module owns the storage of per-observation soft assignments over a
//! library of `K` candidate manifold-atoms (see [`crate::atom_selection`]
//! for the surrounding selection layer). The two key types are:
//!
//! * [`BitVec`] — a minimal dependency-free bitset used to record the *active
//!   support* `S_n ⊆ {0, …, K−1}` of each observation. We avoid pulling in
//!   the external `bitvec` crate to keep this module aligned with the rest of
//!   `gam`'s "no extra deps for new primitives" policy.
//! * [`SparseAtomCode`] — the per-point pair `(active_mask, weights)` whose
//!   semantics are documented on the type. Reconstruction at point `n` is
//!
//!   ```text
//!   Ẑ_n  =  Σ_{k ∈ S_n}  w_{n,k}  ·  decoder_k(t_{n,k})
//!   ```
//!
//!   so `weights[k]` is meaningful only when `active_mask.get(k) == true`.
//!   We store `weights` densely (`Vec<f64>` of length `K`) rather than
//!   sparsely; for the typical SAE workload `K` is small (tens to low
//!   hundreds), and the dense layout lets us reuse [`ndarray`] views and
//!   simple BLAS-shaped loops downstream. The mask carries the discrete
//!   active-set information; the weights carry the soft amplitudes.
//!
//! ## Per-point block locality (arrow structure)
//!
//! Each [`SparseAtomCode`] is the per-row ext-coordinate block for observation `n`
//! restricted to the `K` atoms. Combined with the per-atom on-manifold
//! coordinate `t_{n,k} ∈ ℝ^{d_k}` (held in
//! [`crate::atom_selection::AtomLibrary`]'s per-atom
//! `LatentCoordValues`), the row-local ext-coordinate vector is
//!
//! ```text
//!   ext_n  =  ( a_{n,1..K}  ;  t_{n,1,·}  ;  …  ;  t_{n,K,·} )
//! ```
//!
//! whose interaction graph with the shared decoder coefficients `B_1..B_K`
//! is exactly the arrow / bordered-Hessian pattern from `latent_coord.md`
//! §2.2. The Schur complement that Piece 1 uses to eliminate β before the
//! per-row solve generalises here with one change: the row-`n` block now
//! couples to *only the active subset* `S_n` of decoder borders, not to all
//! K of them. That is the structural fact this module records.

use ndarray::Array1;

/// Minimal bit-vector. Backing storage is `Vec<u64>` words.
///
/// We expose only the operations the atom-selection layer needs: construction,
/// `get`, `set`, `count_ones`, and iteration of set indices. This is
/// deliberately tiny — adding the external `bitvec` crate would be overkill
/// for a few hundred bits per observation.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct BitVec {
    words: Vec<u64>,
    len: usize,
}

impl BitVec {
    /// All-zero bitset of length `len`.
    pub fn zeros(len: usize) -> Self {
        let words = vec![0u64; len.div_ceil(64)];
        Self { words, len }
    }

    /// All-ones bitset of length `len`.
    pub fn ones(len: usize) -> Self {
        let mut bv = Self::zeros(len);
        for i in 0..len {
            bv.set(i, true);
        }
        bv
    }

    pub fn len(&self) -> usize {
        self.len
    }

    pub fn is_empty(&self) -> bool {
        self.len == 0
    }

    #[inline]
    pub fn get(&self, i: usize) -> bool {
        assert!(
            i < self.len,
            "BitVec::get index {i} out of bounds {}",
            self.len
        );
        let (w, b) = (i / 64, i % 64);
        (self.words[w] >> b) & 1 == 1
    }

    #[inline]
    pub fn set(&mut self, i: usize, v: bool) {
        assert!(
            i < self.len,
            "BitVec::set index {i} out of bounds {}",
            self.len
        );
        let (w, b) = (i / 64, i % 64);
        if v {
            self.words[w] |= 1u64 << b;
        } else {
            self.words[w] &= !(1u64 << b);
        }
    }

    /// Number of set bits.
    pub fn count_ones(&self) -> usize {
        self.words.iter().map(|w| w.count_ones() as usize).sum()
    }

    /// Iterator over set indices in ascending order.
    pub fn iter_ones(&self) -> impl Iterator<Item = usize> + '_ {
        (0..self.len).filter(move |&i| self.get(i))
    }

    /// Zero all bits in place.
    pub fn clear(&mut self) {
        for w in self.words.iter_mut() {
            *w = 0;
        }
    }
}

/// Per-point sparse code over `K` candidate atoms.
///
/// Invariants (checked in debug builds):
///
/// * `active_mask.len() == weights.len() == K`.
/// * For any `k` with `active_mask.get(k) == false`, the value `weights[k]`
///   is a nuisance — it must not influence reconstruction. Selection
///   strategies that lower a weight to zero (e.g. [`crate::atom_selection::AtomSelectionStrategy`]'s
///   `L1Relaxed` after thresholding) are responsible for clearing the
///   corresponding mask bit *and* zeroing `weights[k]`.
///
/// We do not require `weights[k] >= 0`; some strategies (entropic softmax,
/// TopK projection) keep the simplex, while others (L¹-relaxed) only enforce
/// non-negativity at the active-set step. The owning
/// [`crate::atom_selection::AtomSelectionStrategy`] documents which
/// invariant it maintains.
#[derive(Debug, Clone)]
pub struct SparseAtomCode {
    /// Length-`K` bitmask of active atoms for this point.
    pub active_mask: BitVec,
    /// Length-`K` dense weight vector. Only entries at active indices are
    /// semantically meaningful.
    pub weights: Vec<f64>,
}

impl SparseAtomCode {
    /// Cold-start: no atoms active, all weights zero.
    pub fn empty(k_atoms: usize) -> Self {
        Self {
            active_mask: BitVec::zeros(k_atoms),
            weights: vec![0.0; k_atoms],
        }
    }

    /// Total number of candidate atoms `K` this code is sized for.
    pub fn k_atoms(&self) -> usize {
        self.weights.len()
    }

    /// Cardinality of the active support `|S_n|`.
    pub fn n_active(&self) -> usize {
        self.active_mask.count_ones()
    }

    /// Sum of active weights. For simplex-projected codes this should be ≈ 1.
    pub fn active_weight_sum(&self) -> f64 {
        self.active_mask.iter_ones().map(|k| self.weights[k]).sum()
    }

    /// Set the weight for atom `k` and mark it active.
    pub fn assign(&mut self, k: usize, w: f64) {
        assert!(k < self.k_atoms());
        self.active_mask.set(k, true);
        self.weights[k] = w;
    }

    /// Deactivate atom `k` and zero its stored weight.
    pub fn deactivate(&mut self, k: usize) {
        assert!(k < self.k_atoms());
        self.active_mask.set(k, false);
        self.weights[k] = 0.0;
    }

    /// Materialize the *effective* weight vector (zeros at inactive indices)
    /// as an owned `Array1`. Useful for matmul-shaped downstream code.
    pub fn effective_weights(&self) -> Array1<f64> {
        let mut out = Array1::<f64>::zeros(self.k_atoms());
        for k in self.active_mask.iter_ones() {
            out[k] = self.weights[k];
        }
        out
    }
}

/// Storage for the per-row codes of all `N` observations.
///
/// Held column-of-structs rather than struct-of-columns: each row's
/// `(active_mask, weights)` lives together because the atom-selection
/// strategies all touch a single row at a time. Cross-row vectorization
/// happens through ndarray views built on demand.
#[derive(Debug, Clone)]
pub struct SparseAtomCodes {
    codes: Vec<SparseAtomCode>,
    k_atoms: usize,
}

impl SparseAtomCodes {
    /// Allocate `n_obs` empty codes, each sized for `k_atoms`.
    pub fn empty(n_obs: usize, k_atoms: usize) -> Self {
        let codes = (0..n_obs).map(|_| SparseAtomCode::empty(k_atoms)).collect();
        Self { codes, k_atoms }
    }

    pub fn n_obs(&self) -> usize {
        self.codes.len()
    }

    pub fn k_atoms(&self) -> usize {
        self.k_atoms
    }

    pub fn row(&self, n: usize) -> &SparseAtomCode {
        &self.codes[n]
    }

    pub fn row_mut(&mut self, n: usize) -> &mut SparseAtomCode {
        &mut self.codes[n]
    }

    pub fn iter(&self) -> impl Iterator<Item = &SparseAtomCode> {
        self.codes.iter()
    }

    pub fn iter_mut(&mut self) -> impl Iterator<Item = &mut SparseAtomCode> {
        self.codes.iter_mut()
    }

    /// Flatten weights into a single `(N, K)` array, with zeros where the
    /// mask is unset. Allocates; intended for diagnostic / post-fit use.
    pub fn weights_matrix(&self) -> ndarray::Array2<f64> {
        let n = self.n_obs();
        let k = self.k_atoms();
        let mut out = ndarray::Array2::<f64>::zeros((n, k));
        for n_idx in 0..n {
            let code = &self.codes[n_idx];
            for kk in code.active_mask.iter_ones() {
                out[[n_idx, kk]] = code.weights[kk];
            }
        }
        out
    }

    /// Co-activation statistics for one atom pair `(a, b)` — the #976
    /// code-dependence trigger. Pure popcount ratios over the active masks:
    /// `P(a|b) = #{rows: a∧b} / #{rows: b}` and symmetrically.
    ///
    /// Two derived readings drive the structure search:
    ///
    /// * [`CoactivationStats::dependence`] (symmetric, the FUSION trigger) —
    ///   independent atoms with marginal activation rates `π_a, π_b` co-activate
    ///   at rate `π_a·π_b`, so both conditionals stay near the marginals; a
    ///   shattered curved family re-encoded as several near-duplicate atoms
    ///   pushes *both* conditionals toward 1.
    /// * [`CoactivationStats::absorption_asymmetry`] (the ABSORPTION-audit
    ///   trigger) — an A⇒B hierarchy where sparsity folded B's content into A
    ///   shows `P(parent|child) ≈ 1` without the converse, so a large asymmetry
    ///   with one conditional near 1 flags the pair for the within-atom
    ///   substructure audit (#907 race on the atom's own code distribution).
    ///
    /// These are *triggers*, not decisions: they rank move proposals
    /// deterministically; acceptance is owned by the e-process gates in
    /// [`gam_solve::structure_search`].
    pub fn coactivation(&self, a: usize, b: usize) -> CoactivationStats {
        assert!(
            a < self.k_atoms && b < self.k_atoms,
            "SparseAtomCodes::coactivation: atoms ({a}, {b}) out of range K={}",
            self.k_atoms
        );
        let n_obs = self.n_obs();
        let mut n_a = 0usize;
        let mut n_b = 0usize;
        let mut n_joint = 0usize;
        for code in &self.codes {
            let on_a = code.active_mask.get(a);
            let on_b = code.active_mask.get(b);
            n_a += usize::from(on_a);
            n_b += usize::from(on_b);
            n_joint += usize::from(on_a && on_b);
        }
        let cond = |joint: usize, marg: usize| {
            if marg == 0 {
                0.0
            } else {
                joint as f64 / marg as f64
            }
        };
        let lift = if n_a == 0 || n_b == 0 || n_obs == 0 {
            0.0
        } else {
            (n_joint as f64 * n_obs as f64) / (n_a as f64 * n_b as f64)
        };
        CoactivationStats {
            n_obs,
            n_a,
            n_b,
            n_joint,
            p_a_given_b: cond(n_joint, n_b),
            p_b_given_a: cond(n_joint, n_a),
            lift,
            weight_correlation: self.weight_codependence(a, b),
        }
    }

    /// #976 — the AMPLITUDE half of the fusion criterion: the Pearson
    /// correlation of the two atoms' activation WEIGHTS over the rows where both
    /// are active. Support co-activation ([`CoactivationStats::dependence`]) only
    /// says the two atoms fire together; it cannot distinguish a single curved
    /// family SHATTERED across two near-duplicate atoms (where moving along the
    /// family smoothly trades amplitude between the pair, so their weights are
    /// strongly — typically negatively — correlated on the joint support) from
    /// two GENUINELY INDEPENDENT atoms that merely happen to co-fire on the same
    /// input class (weights uncorrelated). The magnitude `|ρ|` of this
    /// correlation is the interaction-evidence the issue's fusion trigger pairs
    /// with code dependence: high support-overlap AND high `|weight_correlation|`
    /// is the shattering signature ("dependent codes + joint interaction
    /// evidence"), whereas high overlap with `|ρ|≈0` is two independent features
    /// that should NOT be fused.
    ///
    /// Returns `0.0` when fewer than two rows are jointly active or when either
    /// atom's weight is constant on the joint support (an undefined correlation
    /// is, for the trigger, "no amplitude dependence detected").
    pub fn weight_codependence(&self, a: usize, b: usize) -> f64 {
        assert!(
            a < self.k_atoms && b < self.k_atoms,
            "SparseAtomCodes::weight_codependence: atoms ({a}, {b}) out of range K={}",
            self.k_atoms
        );
        let mut wa = Vec::new();
        let mut wb = Vec::new();
        for code in &self.codes {
            if code.active_mask.get(a) && code.active_mask.get(b) {
                wa.push(code.weights[a]);
                wb.push(code.weights[b]);
            }
        }
        let m = wa.len();
        if m < 2 {
            return 0.0;
        }
        let inv = 1.0 / m as f64;
        let mean_a: f64 = wa.iter().sum::<f64>() * inv;
        let mean_b: f64 = wb.iter().sum::<f64>() * inv;
        let mut cov = 0.0_f64;
        let mut var_a = 0.0_f64;
        let mut var_b = 0.0_f64;
        for i in 0..m {
            let da = wa[i] - mean_a;
            let db = wb[i] - mean_b;
            cov += da * db;
            var_a += da * da;
            var_b += db * db;
        }
        if !(var_a > 0.0 && var_b > 0.0) {
            return 0.0;
        }
        let rho = cov / (var_a.sqrt() * var_b.sqrt());
        // Numerical clamp: accumulation can nudge a perfect ±1 a hair past the
        // bound.
        rho.clamp(-1.0, 1.0)
    }
}

/// Pairwise co-activation summary for two atoms (see
/// [`SparseAtomCodes::coactivation`]). All probabilities are empirical
/// popcount ratios over the active-support masks.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CoactivationStats {
    /// Total number of observations the codes cover.
    pub n_obs: usize,
    /// Rows where atom `a` is active.
    pub n_a: usize,
    /// Rows where atom `b` is active.
    pub n_b: usize,
    /// Rows where both are active.
    pub n_joint: usize,
    /// `P(a active | b active)`; `0` when `b` is never active.
    pub p_a_given_b: f64,
    /// `P(b active | a active)`; `0` when `a` is never active.
    pub p_b_given_a: f64,
    /// `P(a∧b) / (P(a)·P(b))`; `1` for independent atoms, `0` when either
    /// marginal is empty.
    pub lift: f64,
    /// Pearson correlation of the two atoms' activation WEIGHTS over the
    /// jointly-active rows (see [`SparseAtomCodes::weight_codependence`]) — the
    /// amplitude/interaction half of the fusion criterion. `0` when the joint
    /// support is too small or a weight is constant there.
    pub weight_correlation: f64,
}

impl CoactivationStats {
    /// Symmetric code dependence `min(P(a|b), P(b|a))` — the canonical-order
    /// trigger for FUSION proposals (descending). Near 0 for independent or
    /// disjoint atoms; near 1 only when the two supports essentially coincide,
    /// which is the shattering signature.
    pub fn dependence(&self) -> f64 {
        self.p_a_given_b.min(self.p_b_given_a)
    }

    /// Conditional asymmetry `|P(a|b) − P(b|a)|` — large when one atom's
    /// support nests inside the other's (the A⇒B absorption signature, where
    /// `P(parent|child) ≈ 1` but not conversely). Flags the pair for a
    /// targeted within-atom substructure audit; it is never itself an
    /// acceptance criterion.
    pub fn absorption_asymmetry(&self) -> f64 {
        (self.p_a_given_b - self.p_b_given_a).abs()
    }

    /// #976 — the combined FUSION evidence: `dependence · |weight_correlation|`.
    /// A fusion proposal needs BOTH halves — the atoms must co-activate (support
    /// dependence) AND their amplitudes must be dependent on the joint support
    /// (the interaction evidence that a single curved family was shattered). Two
    /// independent atoms that happen to co-fire score near 0 on the second factor
    /// and so are NOT proposed for fusion even at high support overlap; a genuine
    /// shattered pair scores high on both. This is the scalar the canonical-order
    /// fusion ranking ("fusions by code dependence descending") should sort on,
    /// and the threshold the e-process acceptance gate guards.
    pub fn fusion_evidence(&self) -> f64 {
        self.dependence() * self.weight_correlation.abs()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn bitvec_basic() {
        let mut bv = BitVec::zeros(70);
        assert_eq!(bv.len(), 70);
        assert!(!bv.get(5));
        bv.set(5, true);
        bv.set(64, true);
        assert!(bv.get(5));
        assert!(bv.get(64));
        assert_eq!(bv.count_ones(), 2);
        let ones: Vec<usize> = bv.iter_ones().collect();
        assert_eq!(ones, vec![5, 64]);
        bv.set(5, false);
        assert_eq!(bv.count_ones(), 1);
    }

    #[test]
    fn sparse_code_assign() {
        let mut c = SparseAtomCode::empty(8);
        c.assign(2, 0.7);
        c.assign(5, 0.3);
        assert_eq!(c.n_active(), 2);
        assert!((c.active_weight_sum() - 1.0).abs() < 1e-12);
        c.deactivate(2);
        assert_eq!(c.n_active(), 1);
        assert_eq!(c.weights[2], 0.0);
    }

    #[test]
    fn codes_matrix_roundtrip() {
        let mut codes = SparseAtomCodes::empty(3, 4);
        codes.row_mut(0).assign(1, 0.5);
        codes.row_mut(2).assign(3, 0.9);
        let m = codes.weights_matrix();
        assert_eq!(m[[0, 1]], 0.5);
        assert_eq!(m[[2, 3]], 0.9);
        assert_eq!(m[[1, 0]], 0.0);
    }

    /// Co-activation triggers separate the three planted regimes: independent
    /// atoms (low dependence), a shattered duplicate pair (dependence ≈ 1,
    /// symmetric), and an absorption hierarchy (high asymmetry, parent
    /// conditional ≈ 1).
    #[test]
    fn coactivation_separates_independent_shattered_and_absorbed() {
        let n = 100usize;
        let mut codes = SparseAtomCodes::empty(n, 4);
        for row in 0..n {
            // Atom 0: active on even rows; atom 1: active on rows ≡ 0 (mod 5)
            // — independent-ish supports (joint = rows ≡ 0 mod 10).
            if row % 2 == 0 {
                codes.row_mut(row).assign(0, 1.0);
            }
            if row % 5 == 0 {
                codes.row_mut(row).assign(1, 1.0);
            }
            // Atoms 2 and 3: a nested pair — 3 (child) active on rows ≡ 0
            // (mod 4), 2 (parent) active whenever 3 is plus half of the rest.
            if row % 4 == 0 || row % 2 == 1 {
                codes.row_mut(row).assign(2, 1.0);
            }
            if row % 4 == 0 {
                codes.row_mut(row).assign(3, 1.0);
            }
        }

        // Independent pair: P(0|1) = 0.5 (even rows among multiples of 5),
        // P(1|0) = 10/50 = 0.2 → low symmetric dependence, lift = 1.
        let indep = codes.coactivation(0, 1);
        assert_eq!(indep.n_joint, 10);
        assert!((indep.p_a_given_b - 0.5).abs() < 1e-12);
        assert!((indep.p_b_given_a - 0.2).abs() < 1e-12);
        assert!((indep.lift - 1.0).abs() < 1e-12);
        assert!(indep.dependence() < 0.25);

        // Nested (absorption-suspect) pair: P(parent|child) = 1, converse
        // small → near-maximal asymmetry.
        let nested = codes.coactivation(2, 3);
        assert!((nested.p_a_given_b - 1.0).abs() < 1e-12);
        assert!(nested.p_b_given_a < 0.5);
        assert!(nested.absorption_asymmetry() > 0.6);

        // Shattered pair: identical supports → dependence = 1, asymmetry = 0.
        let mut dup = SparseAtomCodes::empty(n, 2);
        for row in (0..n).step_by(3) {
            dup.row_mut(row).assign(0, 1.0);
            dup.row_mut(row).assign(1, 1.0);
        }
        let shat = dup.coactivation(0, 1);
        assert!((shat.dependence() - 1.0).abs() < 1e-12);
        assert!(shat.absorption_asymmetry() < 1e-12);

        // Empty marginals are total, not NaN.
        let empty = SparseAtomCodes::empty(4, 2).coactivation(0, 1);
        assert_eq!(empty.dependence(), 0.0);
        assert_eq!(empty.lift, 0.0);
    }

    /// #976 — the fusion criterion's discriminating power. Support co-activation
    /// alone cannot tell a SHATTERED curved family (one manifold smeared across
    /// two near-duplicate atoms) from two GENUINELY INDEPENDENT atoms that fire
    /// on the same input class: BOTH can have identical supports (dependence = 1).
    /// The amplitude half — [`SparseAtomCodes::weight_codependence`] — separates
    /// them: a shattered pair trades activation weight smoothly as it moves along
    /// the family, so the pair's weights are strongly correlated on the joint
    /// support, while independent co-active atoms have uncorrelated weights.
    /// `fusion_evidence` (dependence · |weight_correlation|) must therefore fire
    /// on the planted shatter and stay near zero for the independent pair.
    #[test]
    fn fusion_criterion_distinguishes_shattered_from_independent_coactive() {
        let n = 120usize;

        // Planted SHATTER: a 1-D family parametrised by t ∈ [0,1] re-encoded as
        // two near-duplicate atoms whose amplitudes are a smooth partition of
        // unity — atom 0 carries `t`, atom 1 carries `1 − t`. Moving along the
        // family trades weight between them, so on the (identical) joint support
        // their weights are PERFECTLY anti-correlated (|ρ| = 1).
        let mut shattered = SparseAtomCodes::empty(n, 2);
        for row in 0..n {
            let t = (row as f64 + 0.5) / n as f64;
            shattered.row_mut(row).assign(0, t);
            shattered.row_mut(row).assign(1, 1.0 - t);
        }
        let shat = shattered.coactivation(0, 1);
        assert!(
            (shat.dependence() - 1.0).abs() < 1e-12,
            "shattered pair shares support: dependence={}",
            shat.dependence()
        );
        assert!(
            shat.weight_correlation < -0.99,
            "shattered family's partition-of-unity weights are anti-correlated on \
             the joint support: weight_correlation={}",
            shat.weight_correlation
        );
        assert!(
            shat.fusion_evidence() > 0.95,
            "fusion evidence must FIRE on a planted shatter: {}",
            shat.fusion_evidence()
        );

        // INDEPENDENT co-active pair: same identical support (dependence = 1),
        // but the two atoms' weights are drawn independently — a deterministic,
        // mutually-uncorrelated pair of sequences (one from a low-frequency
        // sinusoid, one from a coprime-frequency sinusoid, phase-offset) so the
        // joint-support weight correlation is ≈ 0.
        let mut independent = SparseAtomCodes::empty(n, 2);
        for row in 0..n {
            let x = row as f64;
            let wa = 0.5 + 0.4 * (2.0 * std::f64::consts::PI * x / 7.0).sin();
            let wb = 0.5 + 0.4 * (2.0 * std::f64::consts::PI * x / 11.0 + 1.3).cos();
            independent.row_mut(row).assign(0, wa);
            independent.row_mut(row).assign(1, wb);
        }
        let indep = independent.coactivation(0, 1);
        assert!(
            (indep.dependence() - 1.0).abs() < 1e-12,
            "independent pair was constructed with identical support: dependence={}",
            indep.dependence()
        );
        assert!(
            indep.weight_correlation.abs() < 0.3,
            "independent co-active atoms have ~uncorrelated weights: \
             weight_correlation={}",
            indep.weight_correlation
        );

        // The criterion SEPARATES them despite identical support overlap.
        assert!(
            shat.fusion_evidence() > 3.0 * indep.fusion_evidence().max(1e-6),
            "fusion evidence must rank the shattered pair far above the independent \
             pair: shattered={}, independent={}",
            shat.fusion_evidence(),
            indep.fusion_evidence()
        );

        // Degenerate joint supports give a defined (zero) amplitude reading.
        let mut tiny = SparseAtomCodes::empty(3, 2);
        tiny.row_mut(0).assign(0, 1.0);
        tiny.row_mut(0).assign(1, 1.0);
        assert_eq!(
            tiny.weight_codependence(0, 1),
            0.0,
            "a single jointly-active row carries no amplitude correlation"
        );
    }
}