Test the new instruments; fix an O(1/n) bias in the degree decomposition

Four tests around today's additions. Two failed on first run and both
were worth having.

The degree decomposition left an O(1/n) residual on a field that is
purely additive: excluding the diagonal makes the two-way design
unbalanced, so one pass of row and column means does not remove a pure
degree effect. Swept to convergence instead. At N=256 the correction
moves the reported variance shares by under 0.001, so the refutation of
the hubness hypothesis stands unchanged -- but the instrument that
produced it now does what it claims.

The other failure was the test's own scale: two random 16-dimensional
subspaces of R^64 overlap above 0.7 by chance, which is why the real
measurements are made at N=256 where the null sits at 1.0.

Co-Authored-By: Claude <noreply@anthropic.com>

1 files changed, 22 insertions, 7 deletions

diff --git a/worldalign/field_anatomy.py b/worldalign/field_anatomy.py
index 6156018..9e04ec6 100644
--- a/worldalign/field_anatomy.py
+++ b/worldalign/field_anatomy.py
@@ -49,23 +49,38 @@ def correlate(first: np.ndarray, second: np.ndarray) -> float:
     return float(np.corrcoef(offdiagonal(first), offdiagonal(second))[0, 1])
 
 
-def degree_part(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+def degree_part(
+    matrix: np.ndarray, sweeps: int = 200
+) -> tuple[np.ndarray, np.ndarray]:
     """Split into the additive row/column-mean model and its residual.
 
     The two-way additive fit m + r_i + c_j is what a pure busyness effect
     produces: every entry explained by how prominent each of its two scenes
     is, with nothing said about the pair.
+
+    The fit is swept to convergence rather than taken in one pass. With the
+    diagonal excluded the design is unbalanced, so a single pass of row and
+    column means leaves an O(1/n) piece of a pure degree effect behind --
+    small, but it is exactly the quantity this function exists to remove.
+    Alternating the two sweeps converges to the correct fit.
     """
     size = len(matrix)
     mask = ~np.eye(size, dtype=bool)
-    work = matrix.copy().astype(np.float64)
-    np.fill_diagonal(work, np.nan)
+    work = np.where(mask, matrix.astype(np.float64), np.nan)
     grand = np.nanmean(work)
-    rows = np.nanmean(work, axis=1, keepdims=True) - grand
-    cols = np.nanmean(work, axis=0, keepdims=True) - grand
+    rows = np.zeros((size, 1))
+    cols = np.zeros((1, size))
+    for _ in range(sweeps):
+        residual = work - (grand + rows + cols)
+        step = np.nanmean(residual, axis=1, keepdims=True)
+        rows = rows + step
+        residual = work - (grand + rows + cols)
+        step_columns = np.nanmean(residual, axis=0, keepdims=True)
+        cols = cols + step_columns
+        if max(np.abs(step).max(), np.abs(step_columns).max()) < 1e-12:
+            break
     fitted = grand + rows + cols
-    residual = np.where(mask, work - fitted, 0.0)
-    return np.where(mask, fitted, 0.0), residual
+    return np.where(mask, fitted, 0.0), np.where(mask, work - fitted, 0.0)
 
 
 def spectral_strip(matrix: np.ndarray, components: int) -> np.ndarray: