gam_sae/atom_codes.rs
1//! Per-point sparse atom codes for multi-manifold reconstruction.
2//!
3//! This module owns the storage of per-observation soft assignments over a
4//! library of `K` candidate manifold-atoms (see [`crate::assignment::SaeAssignment`]
5//! for the surrounding selection/gate layer). The two key types are:
6//!
7//! * [`BitVec`] — a minimal dependency-free bitset used to record the *active
8//! support* `S_n ⊆ {0, …, K−1}` of each observation. We avoid pulling in
9//! the external `bitvec` crate to keep this module aligned with the rest of
10//! `gam`'s "no extra deps for new primitives" policy.
11//! * [`SparseAtomCode`] — the per-point pair `(active_mask, weights)` whose
12//! semantics are documented on the type. Reconstruction at point `n` is
13//!
14//! ```text
15//! Ẑ_n = Σ_{k ∈ S_n} w_{n,k} · decoder_k(t_{n,k})
16//! ```
17//!
18//! so `weights[k]` is meaningful only when `active_mask.get(k) == true`.
19//! We store `weights` densely (`Vec<f64>` of length `K`) rather than
20//! sparsely; for the typical SAE workload `K` is small (tens to low
21//! hundreds), and the dense layout lets us reuse [`ndarray`] views and
22//! simple BLAS-shaped loops downstream. The mask carries the discrete
23//! active-set information; the weights carry the soft amplitudes.
24//!
25//! ## Per-point block locality (arrow structure)
26//!
27//! Each [`SparseAtomCode`] is the per-row ext-coordinate block for observation `n`
28//! restricted to the `K` atoms. Combined with the per-atom on-manifold
29//! coordinate `t_{n,k} ∈ ℝ^{d_k}` (held in the per-atom latent-coordinate
30//! blocks of [`crate::assignment::SaeAssignment`]), the row-local
31//! ext-coordinate vector is
32//!
33//! ```text
34//! ext_n = ( a_{n,1..K} ; t_{n,1,·} ; … ; t_{n,K,·} )
35//! ```
36//!
37//! whose interaction graph with the shared decoder coefficients `B_1..B_K`
38//! is exactly the arrow / bordered-Hessian pattern from `latent_coord.md`
39//! §2.2. The Schur complement that Piece 1 uses to eliminate β before the
40//! per-row solve generalises here with one change: the row-`n` block now
41//! couples to *only the active subset* `S_n` of decoder borders, not to all
42//! K of them. That is the structural fact this module records.
43
44use ndarray::Array1;
45
46/// Minimal bit-vector. Backing storage is `Vec<u64>` words.
47///
48/// We expose only the operations the atom-selection layer needs: construction,
49/// `get`, `set`, `count_ones`, and iteration of set indices. This is
50/// deliberately tiny — adding the external `bitvec` crate would be overkill
51/// for a few hundred bits per observation.
52#[derive(Debug, Clone, PartialEq, Eq)]
53pub struct BitVec {
54 words: Vec<u64>,
55 len: usize,
56}
57
58impl BitVec {
59 /// All-zero bitset of length `len`.
60 pub fn zeros(len: usize) -> Self {
61 let words = vec![0u64; len.div_ceil(64)];
62 Self { words, len }
63 }
64
65 /// All-ones bitset of length `len`.
66 pub fn ones(len: usize) -> Self {
67 let mut bv = Self::zeros(len);
68 for i in 0..len {
69 bv.set(i, true);
70 }
71 bv
72 }
73
74 pub fn len(&self) -> usize {
75 self.len
76 }
77
78 pub fn is_empty(&self) -> bool {
79 self.len == 0
80 }
81
82 #[inline]
83 pub fn get(&self, i: usize) -> bool {
84 assert!(
85 i < self.len,
86 "BitVec::get index {i} out of bounds {}",
87 self.len
88 );
89 let (w, b) = (i / 64, i % 64);
90 (self.words[w] >> b) & 1 == 1
91 }
92
93 #[inline]
94 pub fn set(&mut self, i: usize, v: bool) {
95 assert!(
96 i < self.len,
97 "BitVec::set index {i} out of bounds {}",
98 self.len
99 );
100 let (w, b) = (i / 64, i % 64);
101 if v {
102 self.words[w] |= 1u64 << b;
103 } else {
104 self.words[w] &= !(1u64 << b);
105 }
106 }
107
108 /// Number of set bits.
109 pub fn count_ones(&self) -> usize {
110 self.words.iter().map(|w| w.count_ones() as usize).sum()
111 }
112
113 /// Iterator over set indices in ascending order.
114 pub fn iter_ones(&self) -> impl Iterator<Item = usize> + '_ {
115 (0..self.len).filter(move |&i| self.get(i))
116 }
117
118 /// Zero all bits in place.
119 pub fn clear(&mut self) {
120 for w in self.words.iter_mut() {
121 *w = 0;
122 }
123 }
124}
125
126/// Per-point sparse code over `K` candidate atoms.
127///
128/// Invariants (checked in debug builds):
129///
130/// * `active_mask.len() == weights.len() == K`.
131/// * For any `k` with `active_mask.get(k) == false`, the value `weights[k]`
132/// is a nuisance — it must not influence reconstruction. Selection
133/// strategies that lower a weight to zero (e.g. an L¹-relaxed gate after
134/// thresholding) are responsible for clearing the
135/// corresponding mask bit *and* zeroing `weights[k]`.
136///
137/// We do not require `weights[k] >= 0`; some strategies (entropic softmax,
138/// TopK projection) keep the simplex, while others (L¹-relaxed) only enforce
139/// non-negativity at the active-set step. The owning
140/// gate/selection layer ([`crate::assignment::SaeAssignment`]) documents which
141/// invariant it maintains.
142#[derive(Debug, Clone)]
143pub struct SparseAtomCode {
144 /// Length-`K` bitmask of active atoms for this point.
145 pub active_mask: BitVec,
146 /// Length-`K` dense weight vector. Only entries at active indices are
147 /// semantically meaningful.
148 pub weights: Vec<f64>,
149}
150
151impl SparseAtomCode {
152 /// Cold-start: no atoms active, all weights zero.
153 pub fn empty(k_atoms: usize) -> Self {
154 Self {
155 active_mask: BitVec::zeros(k_atoms),
156 weights: vec![0.0; k_atoms],
157 }
158 }
159
160 /// Total number of candidate atoms `K` this code is sized for.
161 pub fn k_atoms(&self) -> usize {
162 self.weights.len()
163 }
164
165 /// Cardinality of the active support `|S_n|`.
166 pub fn n_active(&self) -> usize {
167 self.active_mask.count_ones()
168 }
169
170 /// Sum of active weights. For simplex-projected codes this should be ≈ 1.
171 pub fn active_weight_sum(&self) -> f64 {
172 self.active_mask.iter_ones().map(|k| self.weights[k]).sum()
173 }
174
175 /// Set the weight for atom `k` and mark it active.
176 pub fn assign(&mut self, k: usize, w: f64) {
177 assert!(k < self.k_atoms());
178 self.active_mask.set(k, true);
179 self.weights[k] = w;
180 }
181
182 /// Deactivate atom `k` and zero its stored weight.
183 pub fn deactivate(&mut self, k: usize) {
184 assert!(k < self.k_atoms());
185 self.active_mask.set(k, false);
186 self.weights[k] = 0.0;
187 }
188
189 /// Materialize the *effective* weight vector (zeros at inactive indices)
190 /// as an owned `Array1`. Useful for matmul-shaped downstream code.
191 pub fn effective_weights(&self) -> Array1<f64> {
192 let mut out = Array1::<f64>::zeros(self.k_atoms());
193 for k in self.active_mask.iter_ones() {
194 out[k] = self.weights[k];
195 }
196 out
197 }
198}
199
200/// Storage for the per-row codes of all `N` observations.
201///
202/// Held column-of-structs rather than struct-of-columns: each row's
203/// `(active_mask, weights)` lives together because the atom-selection
204/// strategies all touch a single row at a time. Cross-row vectorization
205/// happens through ndarray views built on demand.
206#[derive(Debug, Clone)]
207pub struct SparseAtomCodes {
208 codes: Vec<SparseAtomCode>,
209 k_atoms: usize,
210}
211
212impl SparseAtomCodes {
213 /// Allocate `n_obs` empty codes, each sized for `k_atoms`.
214 pub fn empty(n_obs: usize, k_atoms: usize) -> Self {
215 let codes = (0..n_obs).map(|_| SparseAtomCode::empty(k_atoms)).collect();
216 Self { codes, k_atoms }
217 }
218
219 pub fn n_obs(&self) -> usize {
220 self.codes.len()
221 }
222
223 pub fn k_atoms(&self) -> usize {
224 self.k_atoms
225 }
226
227 pub fn row(&self, n: usize) -> &SparseAtomCode {
228 &self.codes[n]
229 }
230
231 pub fn row_mut(&mut self, n: usize) -> &mut SparseAtomCode {
232 &mut self.codes[n]
233 }
234
235 pub fn iter(&self) -> impl Iterator<Item = &SparseAtomCode> {
236 self.codes.iter()
237 }
238
239 pub fn iter_mut(&mut self) -> impl Iterator<Item = &mut SparseAtomCode> {
240 self.codes.iter_mut()
241 }
242
243 /// Flatten weights into a single `(N, K)` array, with zeros where the
244 /// mask is unset. Allocates; intended for diagnostic / post-fit use.
245 pub fn weights_matrix(&self) -> ndarray::Array2<f64> {
246 let n = self.n_obs();
247 let k = self.k_atoms();
248 let mut out = ndarray::Array2::<f64>::zeros((n, k));
249 for n_idx in 0..n {
250 let code = &self.codes[n_idx];
251 for kk in code.active_mask.iter_ones() {
252 out[[n_idx, kk]] = code.weights[kk];
253 }
254 }
255 out
256 }
257
258 /// Co-activation statistics for one atom pair `(a, b)` — the #976
259 /// code-dependence trigger. Pure popcount ratios over the active masks:
260 /// `P(a|b) = #{rows: a∧b} / #{rows: b}` and symmetrically.
261 ///
262 /// Two derived readings drive the structure search:
263 ///
264 /// * [`CoactivationStats::dependence`] (symmetric, the FUSION trigger) —
265 /// independent atoms with marginal activation rates `π_a, π_b` co-activate
266 /// at rate `π_a·π_b`, so both conditionals stay near the marginals; a
267 /// shattered curved family re-encoded as several near-duplicate atoms
268 /// pushes *both* conditionals toward 1.
269 /// * [`CoactivationStats::absorption_asymmetry`] (the ABSORPTION-audit
270 /// trigger) — an A⇒B hierarchy where sparsity folded B's content into A
271 /// shows `P(parent|child) ≈ 1` without the converse, so a large asymmetry
272 /// with one conditional near 1 flags the pair for the within-atom
273 /// substructure audit (#907 race on the atom's own code distribution).
274 ///
275 /// These are *triggers*, not decisions: they rank move proposals
276 /// deterministically; acceptance is owned by the e-process gates in
277 /// [`gam_solve::structure_search`].
278 pub fn coactivation(&self, a: usize, b: usize) -> CoactivationStats {
279 assert!(
280 a < self.k_atoms && b < self.k_atoms,
281 "SparseAtomCodes::coactivation: atoms ({a}, {b}) out of range K={}",
282 self.k_atoms
283 );
284 let n_obs = self.n_obs();
285 let mut n_a = 0usize;
286 let mut n_b = 0usize;
287 let mut n_joint = 0usize;
288 for code in &self.codes {
289 let on_a = code.active_mask.get(a);
290 let on_b = code.active_mask.get(b);
291 n_a += usize::from(on_a);
292 n_b += usize::from(on_b);
293 n_joint += usize::from(on_a && on_b);
294 }
295 let cond = |joint: usize, marg: usize| {
296 if marg == 0 {
297 0.0
298 } else {
299 joint as f64 / marg as f64
300 }
301 };
302 let lift = if n_a == 0 || n_b == 0 || n_obs == 0 {
303 0.0
304 } else {
305 (n_joint as f64 * n_obs as f64) / (n_a as f64 * n_b as f64)
306 };
307 CoactivationStats {
308 n_obs,
309 n_a,
310 n_b,
311 n_joint,
312 p_a_given_b: cond(n_joint, n_b),
313 p_b_given_a: cond(n_joint, n_a),
314 lift,
315 weight_correlation: self.weight_codependence(a, b),
316 }
317 }
318
319 /// #976 — the AMPLITUDE half of the fusion criterion: the Pearson
320 /// correlation of the two atoms' activation WEIGHTS over the rows where both
321 /// are active. Support co-activation ([`CoactivationStats::dependence`]) only
322 /// says the two atoms fire together; it cannot distinguish a single curved
323 /// family SHATTERED across two near-duplicate atoms (where moving along the
324 /// family smoothly trades amplitude between the pair, so their weights are
325 /// strongly — typically negatively — correlated on the joint support) from
326 /// two GENUINELY INDEPENDENT atoms that merely happen to co-fire on the same
327 /// input class (weights uncorrelated). The magnitude `|ρ|` of this
328 /// correlation is the interaction-evidence the issue's fusion trigger pairs
329 /// with code dependence: high support-overlap AND high `|weight_correlation|`
330 /// is the shattering signature ("dependent codes + joint interaction
331 /// evidence"), whereas high overlap with `|ρ|≈0` is two independent features
332 /// that should NOT be fused.
333 ///
334 /// Returns `0.0` when fewer than two rows are jointly active or when either
335 /// atom's weight is constant on the joint support (an undefined correlation
336 /// is, for the trigger, "no amplitude dependence detected").
337 pub fn weight_codependence(&self, a: usize, b: usize) -> f64 {
338 assert!(
339 a < self.k_atoms && b < self.k_atoms,
340 "SparseAtomCodes::weight_codependence: atoms ({a}, {b}) out of range K={}",
341 self.k_atoms
342 );
343 let mut wa = Vec::new();
344 let mut wb = Vec::new();
345 for code in &self.codes {
346 if code.active_mask.get(a) && code.active_mask.get(b) {
347 wa.push(code.weights[a]);
348 wb.push(code.weights[b]);
349 }
350 }
351 let m = wa.len();
352 if m < 2 {
353 return 0.0;
354 }
355 let inv = 1.0 / m as f64;
356 let mean_a: f64 = wa.iter().sum::<f64>() * inv;
357 let mean_b: f64 = wb.iter().sum::<f64>() * inv;
358 let mut cov = 0.0_f64;
359 let mut var_a = 0.0_f64;
360 let mut var_b = 0.0_f64;
361 for i in 0..m {
362 let da = wa[i] - mean_a;
363 let db = wb[i] - mean_b;
364 cov += da * db;
365 var_a += da * da;
366 var_b += db * db;
367 }
368 if !(var_a > 0.0 && var_b > 0.0) {
369 return 0.0;
370 }
371 let rho = cov / (var_a.sqrt() * var_b.sqrt());
372 // Numerical clamp: accumulation can nudge a perfect ±1 a hair past the
373 // bound.
374 rho.clamp(-1.0, 1.0)
375 }
376}
377
378/// Pairwise co-activation summary for two atoms (see
379/// [`SparseAtomCodes::coactivation`]). All probabilities are empirical
380/// popcount ratios over the active-support masks.
381#[derive(Clone, Copy, Debug, PartialEq)]
382pub struct CoactivationStats {
383 /// Total number of observations the codes cover.
384 pub n_obs: usize,
385 /// Rows where atom `a` is active.
386 pub n_a: usize,
387 /// Rows where atom `b` is active.
388 pub n_b: usize,
389 /// Rows where both are active.
390 pub n_joint: usize,
391 /// `P(a active | b active)`; `0` when `b` is never active.
392 pub p_a_given_b: f64,
393 /// `P(b active | a active)`; `0` when `a` is never active.
394 pub p_b_given_a: f64,
395 /// `P(a∧b) / (P(a)·P(b))`; `1` for independent atoms, `0` when either
396 /// marginal is empty.
397 pub lift: f64,
398 /// Pearson correlation of the two atoms' activation WEIGHTS over the
399 /// jointly-active rows (see [`SparseAtomCodes::weight_codependence`]) — the
400 /// amplitude/interaction half of the fusion criterion. `0` when the joint
401 /// support is too small or a weight is constant there.
402 pub weight_correlation: f64,
403}
404
405impl CoactivationStats {
406 /// Symmetric code dependence `min(P(a|b), P(b|a))` — the canonical-order
407 /// trigger for FUSION proposals (descending). Near 0 for independent or
408 /// disjoint atoms; near 1 only when the two supports essentially coincide,
409 /// which is the shattering signature.
410 pub fn dependence(&self) -> f64 {
411 self.p_a_given_b.min(self.p_b_given_a)
412 }
413
414 /// Conditional asymmetry `|P(a|b) − P(b|a)|` — large when one atom's
415 /// support nests inside the other's (the A⇒B absorption signature, where
416 /// `P(parent|child) ≈ 1` but not conversely). Flags the pair for a
417 /// targeted within-atom substructure audit; it is never itself an
418 /// acceptance criterion.
419 pub fn absorption_asymmetry(&self) -> f64 {
420 (self.p_a_given_b - self.p_b_given_a).abs()
421 }
422
423 /// #976 — the combined FUSION evidence: `dependence · |weight_correlation|`.
424 /// A fusion proposal needs BOTH halves — the atoms must co-activate (support
425 /// dependence) AND their amplitudes must be dependent on the joint support
426 /// (the interaction evidence that a single curved family was shattered). Two
427 /// independent atoms that happen to co-fire score near 0 on the second factor
428 /// and so are NOT proposed for fusion even at high support overlap; a genuine
429 /// shattered pair scores high on both. This is the scalar the canonical-order
430 /// fusion ranking ("fusions by code dependence descending") should sort on,
431 /// and the threshold the e-process acceptance gate guards.
432 pub fn fusion_evidence(&self) -> f64 {
433 self.dependence() * self.weight_correlation.abs()
434 }
435}
436
437#[cfg(test)]
438mod tests {
439 use super::*;
440
441 #[test]
442 fn bitvec_basic() {
443 let mut bv = BitVec::zeros(70);
444 assert_eq!(bv.len(), 70);
445 assert!(!bv.get(5));
446 bv.set(5, true);
447 bv.set(64, true);
448 assert!(bv.get(5));
449 assert!(bv.get(64));
450 assert_eq!(bv.count_ones(), 2);
451 let ones: Vec<usize> = bv.iter_ones().collect();
452 assert_eq!(ones, vec![5, 64]);
453 bv.set(5, false);
454 assert_eq!(bv.count_ones(), 1);
455 }
456
457 #[test]
458 fn sparse_code_assign() {
459 let mut c = SparseAtomCode::empty(8);
460 c.assign(2, 0.7);
461 c.assign(5, 0.3);
462 assert_eq!(c.n_active(), 2);
463 assert!((c.active_weight_sum() - 1.0).abs() < 1e-12);
464 c.deactivate(2);
465 assert_eq!(c.n_active(), 1);
466 assert_eq!(c.weights[2], 0.0);
467 }
468
469 #[test]
470 fn codes_matrix_roundtrip() {
471 let mut codes = SparseAtomCodes::empty(3, 4);
472 codes.row_mut(0).assign(1, 0.5);
473 codes.row_mut(2).assign(3, 0.9);
474 let m = codes.weights_matrix();
475 assert_eq!(m[[0, 1]], 0.5);
476 assert_eq!(m[[2, 3]], 0.9);
477 assert_eq!(m[[1, 0]], 0.0);
478 }
479
480 /// Co-activation triggers separate the three planted regimes: independent
481 /// atoms (low dependence), a shattered duplicate pair (dependence ≈ 1,
482 /// symmetric), and an absorption hierarchy (high asymmetry, parent
483 /// conditional ≈ 1).
484 #[test]
485 fn coactivation_separates_independent_shattered_and_absorbed() {
486 let n = 100usize;
487 let mut codes = SparseAtomCodes::empty(n, 4);
488 for row in 0..n {
489 // Atom 0: active on even rows; atom 1: active on rows ≡ 0 (mod 5)
490 // — independent-ish supports (joint = rows ≡ 0 mod 10).
491 if row % 2 == 0 {
492 codes.row_mut(row).assign(0, 1.0);
493 }
494 if row % 5 == 0 {
495 codes.row_mut(row).assign(1, 1.0);
496 }
497 // Atoms 2 and 3: a nested pair — 3 (child) active on rows ≡ 0
498 // (mod 4), 2 (parent) active whenever 3 is plus half of the rest.
499 if row % 4 == 0 || row % 2 == 1 {
500 codes.row_mut(row).assign(2, 1.0);
501 }
502 if row % 4 == 0 {
503 codes.row_mut(row).assign(3, 1.0);
504 }
505 }
506
507 // Independent pair: P(0|1) = 0.5 (even rows among multiples of 5),
508 // P(1|0) = 10/50 = 0.2 → low symmetric dependence, lift = 1.
509 let indep = codes.coactivation(0, 1);
510 assert_eq!(indep.n_joint, 10);
511 assert!((indep.p_a_given_b - 0.5).abs() < 1e-12);
512 assert!((indep.p_b_given_a - 0.2).abs() < 1e-12);
513 assert!((indep.lift - 1.0).abs() < 1e-12);
514 assert!(indep.dependence() < 0.25);
515
516 // Nested (absorption-suspect) pair: P(parent|child) = 1, converse
517 // small → near-maximal asymmetry.
518 let nested = codes.coactivation(2, 3);
519 assert!((nested.p_a_given_b - 1.0).abs() < 1e-12);
520 assert!(nested.p_b_given_a < 0.5);
521 assert!(nested.absorption_asymmetry() > 0.6);
522
523 // Shattered pair: identical supports → dependence = 1, asymmetry = 0.
524 let mut dup = SparseAtomCodes::empty(n, 2);
525 for row in (0..n).step_by(3) {
526 dup.row_mut(row).assign(0, 1.0);
527 dup.row_mut(row).assign(1, 1.0);
528 }
529 let shat = dup.coactivation(0, 1);
530 assert!((shat.dependence() - 1.0).abs() < 1e-12);
531 assert!(shat.absorption_asymmetry() < 1e-12);
532
533 // Empty marginals are total, not NaN.
534 let empty = SparseAtomCodes::empty(4, 2).coactivation(0, 1);
535 assert_eq!(empty.dependence(), 0.0);
536 assert_eq!(empty.lift, 0.0);
537 }
538
539 /// #976 — the fusion criterion's discriminating power. Support co-activation
540 /// alone cannot tell a SHATTERED curved family (one manifold smeared across
541 /// two near-duplicate atoms) from two GENUINELY INDEPENDENT atoms that fire
542 /// on the same input class: BOTH can have identical supports (dependence = 1).
543 /// The amplitude half — [`SparseAtomCodes::weight_codependence`] — separates
544 /// them: a shattered pair trades activation weight smoothly as it moves along
545 /// the family, so the pair's weights are strongly correlated on the joint
546 /// support, while independent co-active atoms have uncorrelated weights.
547 /// `fusion_evidence` (dependence · |weight_correlation|) must therefore fire
548 /// on the planted shatter and stay near zero for the independent pair.
549 #[test]
550 fn fusion_criterion_distinguishes_shattered_from_independent_coactive() {
551 let n = 120usize;
552
553 // Planted SHATTER: a 1-D family parametrised by t ∈ [0,1] re-encoded as
554 // two near-duplicate atoms whose amplitudes are a smooth partition of
555 // unity — atom 0 carries `t`, atom 1 carries `1 − t`. Moving along the
556 // family trades weight between them, so on the (identical) joint support
557 // their weights are PERFECTLY anti-correlated (|ρ| = 1).
558 let mut shattered = SparseAtomCodes::empty(n, 2);
559 for row in 0..n {
560 let t = (row as f64 + 0.5) / n as f64;
561 shattered.row_mut(row).assign(0, t);
562 shattered.row_mut(row).assign(1, 1.0 - t);
563 }
564 let shat = shattered.coactivation(0, 1);
565 assert!(
566 (shat.dependence() - 1.0).abs() < 1e-12,
567 "shattered pair shares support: dependence={}",
568 shat.dependence()
569 );
570 assert!(
571 shat.weight_correlation < -0.99,
572 "shattered family's partition-of-unity weights are anti-correlated on \
573 the joint support: weight_correlation={}",
574 shat.weight_correlation
575 );
576 assert!(
577 shat.fusion_evidence() > 0.95,
578 "fusion evidence must FIRE on a planted shatter: {}",
579 shat.fusion_evidence()
580 );
581
582 // INDEPENDENT co-active pair: same identical support (dependence = 1),
583 // but the two atoms' weights are drawn independently — a deterministic,
584 // mutually-uncorrelated pair of sequences (one from a low-frequency
585 // sinusoid, one from a coprime-frequency sinusoid, phase-offset) so the
586 // joint-support weight correlation is ≈ 0.
587 let mut independent = SparseAtomCodes::empty(n, 2);
588 for row in 0..n {
589 let x = row as f64;
590 let wa = 0.5 + 0.4 * (2.0 * std::f64::consts::PI * x / 7.0).sin();
591 let wb = 0.5 + 0.4 * (2.0 * std::f64::consts::PI * x / 11.0 + 1.3).cos();
592 independent.row_mut(row).assign(0, wa);
593 independent.row_mut(row).assign(1, wb);
594 }
595 let indep = independent.coactivation(0, 1);
596 assert!(
597 (indep.dependence() - 1.0).abs() < 1e-12,
598 "independent pair was constructed with identical support: dependence={}",
599 indep.dependence()
600 );
601 assert!(
602 indep.weight_correlation.abs() < 0.3,
603 "independent co-active atoms have ~uncorrelated weights: \
604 weight_correlation={}",
605 indep.weight_correlation
606 );
607
608 // The criterion SEPARATES them despite identical support overlap.
609 assert!(
610 shat.fusion_evidence() > 3.0 * indep.fusion_evidence().max(1e-6),
611 "fusion evidence must rank the shattered pair far above the independent \
612 pair: shattered={}, independent={}",
613 shat.fusion_evidence(),
614 indep.fusion_evidence()
615 );
616
617 // Degenerate joint supports give a defined (zero) amplitude reading.
618 let mut tiny = SparseAtomCodes::empty(3, 2);
619 tiny.row_mut(0).assign(0, 1.0);
620 tiny.row_mut(0).assign(1, 1.0);
621 assert_eq!(
622 tiny.weight_codependence(0, 1),
623 0.0,
624 "a single jointly-active row carries no amplitude correlation"
625 );
626 }
627}