Retire the correlation gate: the shared spectrum governs recovery

A controlled truncation refutes the project's central go/no-go rule.
Projecting the recovering synthetic fields to rank r holds the field
correlation at 0.902-0.929 while recovery moves 6.2% -> 12.9% -> 95.6%
across ranks 4, 8, 16. A field past the supposed 0.9 threshold recovers
13%, so correlation neither predicts nor forbids recovery and the width
of the shared spectrum is what moves it.

The gate becomes a joint condition on correlation and shared width,
measured by principal angles against a scene-shuffled null. Neither
suffices alone: 18 shared directions at 0.508 fails, 11 at 0.902 fails.

With the old gate retired, natural data was finally searched: 0.0000
against 0.0039 chance. The old verdict was right, its reasoning was not.

Also closes route D by measurement. rho_IT ~ sqrt(4 log N / N) rises as N
falls, and at N = 16 through 96 the deepest state a strong searcher
reaches is deeper than the truth in 3/3 replicates at every size.

Free gains: eigenvalue-weighted projection over a wide basis with
128-dim text vectors takes the correlation 0.656 -> 0.716 and shared
width 10 -> 16. Hubness refuted as an inflation hypothesis.

Moving the per-image segmentation eigendecomposition onto the GPU cut
batch time from 130s to 1.9s.

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

9 files changed, 1097 insertions, 30 deletions

diff --git a/worldalign/diversity_subset.py b/worldalign/diversity_subset.py
new file mode 100644
index 0000000..fd74d0c
--- /dev/null
+++ b/worldalign/diversity_subset.py
@@ -0,0 +1,171 @@
+"""Does choosing mutually dissimilar scenes buy anything, and at what price?
+
+If the cross-modal agreement lives in a narrow subspace, then which scenes
+are matched matters as much as how many. Scenes that sit close together in
+that subspace are the decoys: nothing in either field separates them, so a
+solver has no way to tell them apart. Spreading the population out should
+help more than enlarging it.
+
+Selection has to be one-sided to stay legal -- picking diverse images and
+diverse captions independently yields two sets that are not in bijection,
+which the permutation formulation cannot express. So this runs as a
+**diagnostic**: scenes are chosen by visual dissimilarity alone and their
+partners are taken from the hidden pairing. That is a protocol violation and
+the result is an upper bound, not a method. It is run first because it is
+cheap and it prices the idea: if diversity does not help even when the
+counterpart set is handed over, the legal unbalanced version is not worth
+building.
+
+Reported against a random-subset control of the same size and the same
+pipeline.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+
+from .common import write_json
+from .spectral_match import grampa
+from .synth_fast_gate import ClosedFormEnergy, all_swaps, steepest_descent
+from .synth_triangle_gate import standardized
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", default="artifacts/vg_5k/natural_pipeline_w64.pt")
+    parser.add_argument("--size", type=int, default=128)
+    parser.add_argument("--trials", type=int, default=5)
+    parser.add_argument("--device", default="cuda:3")
+    parser.add_argument("--output", default="artifacts/vg_5k/diversity_subset.json")
+    return parser.parse_args()
+
+
+def offdiagonal(matrix: np.ndarray) -> np.ndarray:
+    return matrix[~np.eye(len(matrix), dtype=bool)]
+
+
+def normalise(matrix: np.ndarray) -> np.ndarray:
+    values = offdiagonal(matrix)
+    out = (matrix - values.mean()) / values.std()
+    np.fill_diagonal(out, 0.0)
+    return out
+
+
+def farthest_point(field: np.ndarray, size: int, start: int) -> np.ndarray:
+    """Greedy farthest-point selection under the visual field's similarity.
+
+    High field entry means similar, so dissimilarity is the negated entry and
+    each new scene is the one whose closest already-chosen neighbour is as far
+    as possible.
+    """
+    chosen = [start]
+    closest = field[start].copy()
+    for _ in range(size - 1):
+        candidate = int(np.argmin(np.where(np.isin(np.arange(len(field)), chosen),
+                                           np.inf, closest)))
+        chosen.append(candidate)
+        closest = np.maximum(closest, field[candidate])
+    return np.array(chosen)
+
+
+def evaluate(visual: np.ndarray, text: np.ndarray, device, trials: int) -> dict:
+    size = len(visual)
+    visual = normalise(visual)
+    text = normalise(text)
+    mask = ~np.eye(size, dtype=bool)
+    rho = float(np.corrcoef(visual[mask], text[mask])[0, 1])
+    swaps = all_swaps(size, device)
+    accuracies = []
+    for trial in range(trials):
+        generator = np.random.default_rng(trial)
+        hidden = generator.permutation(size)
+        shuffled = text[np.ix_(hidden, hidden)]
+        columns = grampa(visual, shuffled, 1.0)
+        energy = ClosedFormEnergy(
+            standardized(torch.from_numpy(shuffled).to(device)).float(),
+            standardized(torch.from_numpy(visual).to(device)).float(),
+            1.0, 1.0, 256,
+        )
+        final, _ = steepest_descent(
+            energy, torch.from_numpy(columns.copy()).to(device), swaps, 3000
+        )
+        accuracies.append(float((hidden[final.cpu().numpy()] == np.arange(size)).mean()))
+    # how much of the correlation survives stripping the leading directions
+    symmetric = (visual + visual.T) / 2.0
+    values, vectors = np.linalg.eigh(symmetric)
+    order = np.argsort(np.abs(values))[::-1][10:]
+    stripped_visual = (vectors[:, order] * values[order]) @ vectors[:, order].T
+    symmetric_text = (text + text.T) / 2.0
+    values_t, vectors_t = np.linalg.eigh(symmetric_text)
+    order_t = np.argsort(np.abs(values_t))[::-1][10:]
+    stripped_text = (vectors_t[:, order_t] * values_t[order_t]) @ vectors_t[:, order_t].T
+    retained = float(
+        np.corrcoef(offdiagonal(stripped_visual), offdiagonal(stripped_text))[0, 1]
+    )
+    return {
+        "correlation": rho,
+        "correlation_after_strip10": retained,
+        "recovery": float(np.mean(accuracies)),
+        "chance": 1.0 / size,
+    }
+
+
+def main() -> None:
+    args = parse_args()
+    state = torch.load(args.fields, map_location="cpu", weights_only=False)
+    visual_full = state["visual_field"].double().numpy()
+    text_full = state["text_field"].double().numpy()
+    population = len(visual_full)
+    device = torch.device(args.device)
+    size = min(args.size, population)
+
+    rows = []
+    for start in range(3):
+        index = farthest_point(visual_full, size, start * 37 % population)
+        result = evaluate(
+            visual_full[np.ix_(index, index)], text_full[np.ix_(index, index)],
+            device, args.trials,
+        )
+        result["selection"] = f"diverse (start {start})"
+        rows.append(result)
+        print(json.dumps(result), flush=True)
+
+    for trial in range(3):
+        generator = np.random.default_rng(100 + trial)
+        index = generator.choice(population, size, replace=False)
+        result = evaluate(
+            visual_full[np.ix_(index, index)], text_full[np.ix_(index, index)],
+            device, args.trials,
+        )
+        result["selection"] = f"random (seed {trial})"
+        rows.append(result)
+        print(json.dumps(result), flush=True)
+
+    diverse = [r for r in rows if r["selection"].startswith("diverse")]
+    random_rows = [r for r in rows if r["selection"].startswith("random")]
+    summary = {
+        "protocol": (
+            "DIAGNOSTIC, NOT A METHOD: scenes selected by visual dissimilarity "
+            "and their partners taken from the hidden pairing. One-sided "
+            "selection cannot be realised without the pairing, so this is an "
+            "upper bound used to price the idea before building the legal "
+            "unbalanced version."
+        ),
+        "population": population,
+        "subset_size": size,
+        "rows": rows,
+        "diverse_mean_correlation": float(np.mean([r["correlation"] for r in diverse])),
+        "random_mean_correlation": float(np.mean([r["correlation"] for r in random_rows])),
+        "diverse_mean_recovery": float(np.mean([r["recovery"] for r in diverse])),
+        "random_mean_recovery": float(np.mean([r["recovery"] for r in random_rows])),
+    }
+    print(json.dumps({k: v for k, v in summary.items() if k != "rows"}))
+    write_json(args.output, summary)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/field_anatomy.py b/worldalign/field_anatomy.py
new file mode 100644
index 0000000..6156018
--- /dev/null
+++ b/worldalign/field_anatomy.py
@@ -0,0 +1,148 @@
+"""What is the field correlation made of, and how much of it can matching use?
+
+A puzzle motivates this. At 256 scenes the information-theoretic threshold is
+0.29 and the fields correlate at 0.656, so the pairing is identifiable by a
+wide margin, yet even at 64 scenes -- where the threshold is 0.51 and the
+search space is trivially small -- recovery reaches 14%. Either the theory is
+the wrong reference or the statistic is not measuring what the theory means.
+
+The suspect is nuisance structure shared by both fields. Scene similarity has
+a component that says only "this scene is busy, so it resembles everything",
+and both modalities see it. That component inflates the correlation over all
+pairs while carrying no information about which scene is which, because it is
+the same for every candidate pairing of a scene to a partner of similar
+busyness. Correlated-matching theory assumes the correlation is in the
+exchangeable noise, so a nuisance-inflated statistic overstates what the
+theory would predict.
+
+This decomposes each field into a degree component -- the rank-one part
+predicted by row and column means -- and the residual, and reports the
+correlation of each. The residual correlation is the honest input to the
+threshold, and if it falls far below 0.656 the recovery gap is explained and
+the target moves.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+
+from .common import write_json
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", default="artifacts/vg_5k/natural_pipeline.pt")
+    parser.add_argument("--label", default="natural")
+    parser.add_argument("--output", default="artifacts/vg_5k/field_anatomy.json")
+    return parser.parse_args()
+
+
+def offdiagonal(matrix: np.ndarray) -> np.ndarray:
+    return matrix[~np.eye(len(matrix), dtype=bool)]
+
+
+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]:
+    """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.
+    """
+    size = len(matrix)
+    mask = ~np.eye(size, dtype=bool)
+    work = matrix.copy().astype(np.float64)
+    np.fill_diagonal(work, np.nan)
+    grand = np.nanmean(work)
+    rows = np.nanmean(work, axis=1, keepdims=True) - grand
+    cols = np.nanmean(work, axis=0, keepdims=True) - grand
+    fitted = grand + rows + cols
+    residual = np.where(mask, work - fitted, 0.0)
+    return np.where(mask, fitted, 0.0), residual
+
+
+def spectral_strip(matrix: np.ndarray, components: int) -> np.ndarray:
+    """Remove the leading eigen-components, a stronger nuisance model."""
+    symmetric = (matrix + matrix.T) / 2.0
+    values, vectors = np.linalg.eigh(symmetric)
+    order = np.argsort(np.abs(values))[::-1]
+    keep = order[components:]
+    return (vectors[:, keep] * values[keep]) @ vectors[:, keep].T
+
+
+def main() -> None:
+    args = parse_args()
+    state = torch.load(args.fields, map_location="cpu", weights_only=False)
+    visual = state["visual_field"].double().numpy()
+    text = state["text_field"].double().numpy()
+    size = len(visual)
+
+    total = correlate(visual, text)
+    visual_fit, visual_residual = degree_part(visual)
+    text_fit, text_residual = degree_part(text)
+
+    rows = [
+        {"component": "raw field", "correlation": total},
+        {"component": "degree model only", "correlation": correlate(visual_fit, text_fit)},
+        {
+            "component": "residual after degree",
+            "correlation": correlate(visual_residual, text_residual),
+        },
+    ]
+    for components in (1, 2, 5, 10):
+        rows.append(
+            {
+                "component": f"residual after top-{components} eigen",
+                "correlation": correlate(
+                    spectral_strip(visual, components), spectral_strip(text, components)
+                ),
+            }
+        )
+
+    # how much of each field's own variance the nuisance model absorbs
+    share = {
+        "visual_degree_variance_share": float(
+            offdiagonal(visual_fit).var() / offdiagonal(visual).var()
+        ),
+        "text_degree_variance_share": float(
+            offdiagonal(text_fit).var() / offdiagonal(text).var()
+        ),
+    }
+
+    residual_rho = rows[2]["correlation"]
+    threshold = float(np.sqrt(4 * np.log(size) / size))
+    summary = {
+        "protocol": (
+            "Both fields split into an additive row/column-mean model and its "
+            "residual, then correlated component by component. Hidden pairs "
+            "are read only to align the two fields."
+        ),
+        "label": args.label,
+        "size": size,
+        "rows": rows,
+        **share,
+        "it_threshold_at_size": threshold,
+        "residual_above_it_threshold": bool(residual_rho > threshold),
+        "reading": (
+            "If the residual correlation is far below the raw correlation, the "
+            "headline statistic is inflated by nuisance structure that carries "
+            "no matching information, and the residual is the number that "
+            "should be compared against the recovery threshold."
+        ),
+    }
+    for row in rows:
+        print(f"{row['component']:<34} rho={row['correlation']:.4f}")
+    print(json.dumps(share))
+    print(f"IT threshold at N={size}: {threshold:.4f}")
+    write_json(args.output, summary)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/natural_objects.py b/worldalign/natural_objects.py
index 65e035b..6e18259 100644
--- a/worldalign/natural_objects.py
+++ b/worldalign/natural_objects.py
@@ -62,7 +62,8 @@ def parse_args() -> argparse.Namespace:
 
 
 def spectral_segments(
-    features: torch.Tensor, grid: int, segments: int, seed: int
+    features: torch.Tensor, grid: int, segments: int, seed: int,
+    device: torch.device | str = "cpu",
 ) -> np.ndarray:
     """Segment a patch grid by clustering the normalised affinity spectrum.
 
@@ -70,22 +71,42 @@ def spectral_segments(
     stay connected, then the leading eigenvectors of the normalised
     Laplacian are clustered. This is the standard unsupervised recipe over
     self-supervised features.
+
+    The eigendecomposition dominates the run: one dense symmetric problem of
+    side grid^2 per image, which on a contended CPU costs more than the
+    forward pass that produced the features. Running it on the accelerator
+    that is already holding the model turns a multi-hour extraction into a
+    short one, and the spatial prior is built once and cached rather than
+    rebuilt per image.
     """
     from sklearn.cluster import KMeans
 
-    normalised = F.normalize(features.double(), dim=-1)
-    affinity = (normalised @ normalised.T).clamp_min(0).numpy()
-    coordinates = np.stack(
-        np.meshgrid(np.arange(grid), np.arange(grid), indexing="ij"), -1
-    ).reshape(-1, 2).astype(np.float64)
-    distance = ((coordinates[:, None, :] - coordinates[None, :, :]) ** 2).sum(-1)
-    affinity = affinity * np.exp(-distance / (2 * (grid / 4.0) ** 2))
+    normalised = F.normalize(features.to(device).double(), dim=-1)
+    affinity = (normalised @ normalised.T).clamp_min(0)
+    affinity = affinity * _spatial_prior(grid, device)
     degree = affinity.sum(1)
-    laplacian = affinity / np.sqrt(np.outer(degree, degree) + 1e-9)
-    values, vectors = np.linalg.eigh(laplacian)
+    laplacian = affinity / torch.sqrt(torch.outer(degree, degree) + 1e-9)
+    values, vectors = torch.linalg.eigh(laplacian)
     embedding = vectors[:, -segments:]
-    embedding /= np.linalg.norm(embedding, axis=1, keepdims=True).clip(1e-9)
-    return KMeans(segments, n_init=10, random_state=seed).fit_predict(embedding)
+    embedding = embedding / embedding.norm(dim=1, keepdim=True).clamp_min(1e-9)
+    return KMeans(segments, n_init=10, random_state=seed).fit_predict(
+        embedding.cpu().numpy()
+    )
+
+
+_PRIOR_CACHE: dict[tuple, torch.Tensor] = {}
+
+
+def _spatial_prior(grid: int, device: torch.device | str) -> torch.Tensor:
+    """Gaussian locality weight on the patch lattice, built once per grid."""
+    key = (grid, str(device))
+    if key not in _PRIOR_CACHE:
+        axis = torch.arange(grid, dtype=torch.float64, device=device)
+        rows, cols = torch.meshgrid(axis, axis, indexing="ij")
+        coordinates = torch.stack([rows.flatten(), cols.flatten()], dim=-1)
+        distance = torch.cdist(coordinates, coordinates) ** 2
+        _PRIOR_CACHE[key] = torch.exp(-distance / (2 * (grid / 4.0) ** 2))
+    return _PRIOR_CACHE[key]
 
 
 def describe_segments(
@@ -220,7 +241,7 @@ def main() -> None:
                 )
             else:
                 labels = spectral_segments(
-                    patches[position], grid, args.segments, args.seed
+                    patches[position], grid, args.segments, args.seed, args.device
                 )
                 described = describe_segments(
                     labels, patches[position], raw[position], grid, args.min_patches
diff --git a/worldalign/natural_pipeline.py b/worldalign/natural_pipeline.py
index c30af4d..1a4f63e 100644
--- a/worldalign/natural_pipeline.py
+++ b/worldalign/natural_pipeline.py
@@ -41,9 +41,18 @@ def parse_args() -> argparse.Namespace:
     parser.add_argument("--min-count", type=int, default=60)
     parser.add_argument("--max-vocabulary", type=int, default=600)
     parser.add_argument("--context-window", type=int, default=2)
-    parser.add_argument("--keep", type=int, default=8)
+    parser.add_argument("--keep", type=int, default=64)
     parser.add_argument("--shrinkage", type=float, default=0.05)
+    parser.add_argument(
+        "--weight-power",
+        type=float,
+        default=0.5,
+        help="Scale each kept direction by its discriminant eigenvalue to "
+        "this power. 0 reproduces the unweighted projection.",
+    )
     parser.add_argument("--views", type=int, default=1)
+    parser.add_argument("--moment-degree", type=int, default=2, choices=[2, 3])
+    parser.add_argument("--third-width", type=int, default=12)
     parser.add_argument("--seed", type=int, default=0)
     parser.add_argument("--output", default="artifacts/vg_5k/natural_pipeline.pt")
     return parser.parse_args()
@@ -77,10 +86,43 @@ def word_vectors(
     return {word: embedding[row] for row, word in enumerate(vocabulary)}
 
 
+def moment_field_degree(
+    sets: list[torch.Tensor], degree: int, third_width: int
+) -> torch.Tensor:
+    """Set kernel carrying moments up to `degree`.
+
+    A kernel carrying third moments spans more directions than one carrying
+    two, so raising the degree raises the field's rank by construction --
+    which is the quantity the rank experiment showed to govern recovery. The
+    third moment is taken over a narrower basis because its size is cubic;
+    the first two are kept at full width.
+    """
+    if degree < 3:
+        return moment_field(sets)
+    phis = []
+    for members in sets:
+        first = members.mean(0)
+        second = (members[:, :, None] * members[:, None, :]).mean(0).flatten()
+        narrow = members[:, :third_width]
+        outer = narrow[:, :, None] * narrow[:, None, :]
+        third = (outer[:, :, :, None] * narrow[:, None, None, :]).mean(0).flatten()
+        phi = torch.cat([first, second, third])
+        phis.append(phi / phi.norm().clamp_min(1e-9))
+    stacked = torch.stack(phis)
+    return stacked @ stacked.T
+
+
 def content_directions(
     sets: list[np.ndarray], shrinkage: float
-) -> tuple[np.ndarray, np.ndarray]:
-    """Directions separating scenes more than they separate parts."""
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Directions separating scenes more than they separate parts.
+
+    Returns the discriminant eigenvalues alongside the basis. Scaling each
+    direction by its own eigenvalue lets a wide basis be kept without the
+    weakly separating directions drowning the strong ones, which is what
+    makes a wide keep beat a narrow one -- and a wide keep is what raises
+    the rank of the resulting field.
+    """
     means = np.stack([item.mean(0) for item in sets])
     centre = means.mean(0)
     within = np.concatenate([item - item.mean(0, keepdims=True) for item in sets])
@@ -91,7 +133,8 @@ def content_directions(
     values, vectors = eigh(
         scatter_between, scatter_within + shrinkage * trace * np.eye(len(scatter_within))
     )
-    return centre, vectors[:, np.argsort(values)[::-1]]
+    order = np.argsort(values)[::-1]
+    return centre, vectors[:, order], np.maximum(values[order], 0.0)
 
 
 def main() -> None:
@@ -144,8 +187,8 @@ def main() -> None:
     language_fit = [item for item in language_fit if item is not None and len(item) > 1]
     vision_fit = [vision_state(p["vision_node_id"], 0) for p in fit]
     vision_fit = [item for item in vision_fit if item is not None and len(item) > 1]
-    language_centre, language_basis = content_directions(language_fit, args.shrinkage)
-    vision_centre, vision_basis = content_directions(vision_fit, args.shrinkage)
+    language_fitted = content_directions(language_fit, args.shrinkage)
+    vision_fitted = content_directions(vision_fit, args.shrinkage)
 
     evaluated = [
         pair for pair in pairs[: args.samples]
@@ -153,25 +196,32 @@ def main() -> None:
     ]
     keep = args.keep
 
-    def project(raw: np.ndarray, centre: np.ndarray, basis: np.ndarray) -> torch.Tensor:
+    def project(raw: np.ndarray, fitted: tuple) -> torch.Tensor:
+        centre, basis, values = fitted
         width = min(keep, basis.shape[1])
+        scale = values[:width] ** args.weight_power if args.weight_power else 1.0
         return F.normalize(
-            torch.tensor((raw - centre) @ basis[:, :width], dtype=torch.float32), dim=-1
+            torch.tensor(
+                (raw - centre) @ basis[:, :width] * scale, dtype=torch.float32
+            ), dim=-1
         )
 
     text_sets = [
-        project(language_state(p["text_node_id"]), language_centre, language_basis)
-        for p in evaluated
+        project(language_state(p["text_node_id"]), language_fitted) for p in evaluated
     ]
     view_fields = []
     for view in range(views):
         sets = [
-            project(vision_state(p["vision_node_id"], view), vision_centre, vision_basis)
+            project(vision_state(p["vision_node_id"], view), vision_fitted)
             for p in evaluated
         ]
-        view_fields.append(moment_field(sets))
+        view_fields.append(
+            moment_field_degree(sets, args.moment_degree, args.third_width)
+        )
     visual_field = torch.stack(view_fields).mean(0)
-    text_field = moment_field(text_sets)
+    text_field = moment_field_degree(
+        text_sets, args.moment_degree, args.third_width
+    )
 
     mask = ~np.eye(len(evaluated), dtype=bool)
     correlation = float(
@@ -193,6 +243,8 @@ def main() -> None:
         "samples": len(evaluated),
         "views": views,
         "kept_directions": keep,
+        "weight_power": args.weight_power,
+        "moment_degree": args.moment_degree,
         "field_correlation_at_truth": correlation,
         "recovery_threshold": 0.9,
     }
diff --git a/worldalign/natural_recovery.py b/worldalign/natural_recovery.py
new file mode 100644
index 0000000..58de5e5
--- /dev/null
+++ b/worldalign/natural_recovery.py
@@ -0,0 +1,130 @@
+"""Blind recovery on natural data, run because the old gate was wrong.
+
+The standing gate said the correlation had to reach 0.9 and natural data sat
+at 0.656, so no recovery run was warranted and none was made. The rank
+experiment retires that rule: a synthetic field truncated to rank eight has a
+correlation of 0.906 and recovers 13%, while the same field at rank sixteen
+recovers 96%. Correlation alone neither predicts nor forbids recovery, and the
+width of the shared spectrum is the variable that moves it.
+
+By the width measure the best natural configuration now sits at sixteen shared
+directions, inside the interval where the synthetic ladder crosses from 13% to
+96%. Its correlation is far lower, so the joint condition is probably not met
+-- but the honest way to find out is to run the search rather than to infer
+the answer from a statistic that has just been shown not to govern it.
+
+Hidden pairs are generated per trial and read only for scoring.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+
+from .common import write_json
+from .spectral_match import grampa, umeyama
+from .synth_fast_gate import ClosedFormEnergy, all_swaps, steepest_descent
+from .synth_triangle_gate import standardized
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", default="artifacts/vg_5k/np_n256.pt")
+    parser.add_argument("--label", default="natural-best")
+    parser.add_argument("--trials", type=int, default=5)
+    parser.add_argument("--etas", type=float, nargs="+", default=[0.2, 0.5, 1.0, 2.0])
+    parser.add_argument("--device", default="cuda:3")
+    parser.add_argument("--output", default="artifacts/vg_5k/natural_recovery.json")
+    return parser.parse_args()
+
+
+def offdiagonal(matrix: np.ndarray) -> np.ndarray:
+    return matrix[~np.eye(len(matrix), dtype=bool)]
+
+
+def normalise(matrix: np.ndarray) -> np.ndarray:
+    values = offdiagonal(matrix)
+    out = (matrix - values.mean()) / values.std()
+    np.fill_diagonal(out, 0.0)
+    return out
+
+
+def main() -> None:
+    args = parse_args()
+    state = torch.load(args.fields, map_location="cpu", weights_only=False)
+    visual = normalise(state["visual_field"].double().numpy())
+    text = normalise(state["text_field"].double().numpy())
+    size = len(visual)
+    device = torch.device(args.device)
+    swaps = all_swaps(size, device)
+    mask = ~np.eye(size, dtype=bool)
+    rho = float(np.corrcoef(visual[mask], text[mask])[0, 1])
+
+    rows = []
+    for trial in range(args.trials):
+        generator = np.random.default_rng(trial)
+        hidden = generator.permutation(size)
+        shuffled = text[np.ix_(hidden, hidden)]
+        energy = ClosedFormEnergy(
+            standardized(torch.from_numpy(shuffled).to(device)).float(),
+            standardized(torch.from_numpy(visual).to(device)).float(),
+            1.0, 1.0, 256,
+        )
+
+        def score(columns: np.ndarray) -> float:
+            return float((hidden[columns] == np.arange(size)).mean())
+
+        # several spectral starts, keep the one the energy prefers -- model
+        # selection by energy only, never by accuracy
+        best_energy, best_columns, best_tag = np.inf, None, ""
+        candidates = [("umeyama", umeyama(visual, shuffled))]
+        candidates += [(f"grampa eta={eta}", grampa(visual, shuffled, eta))
+                       for eta in args.etas]
+        for tag, columns in candidates:
+            start = torch.from_numpy(columns.copy()).to(device)
+            value = float(energy.energy(start[None])[0])
+            if value < best_energy:
+                best_energy, best_columns, best_tag = value, columns, tag
+
+        initial = score(best_columns)
+        final, _ = steepest_descent(
+            energy, torch.from_numpy(best_columns.copy()).to(device), swaps, 5000
+        )
+        refined = score(final.cpu().numpy())
+        rows.append(
+            {
+                "trial": trial,
+                "chosen_start": best_tag,
+                "spectral_accuracy": initial,
+                "refined_accuracy": refined,
+            }
+        )
+        print(
+            f"trial {trial}: start={best_tag:<14} spectral={initial:.4f} "
+            f"refined={refined:.4f}  (chance {1/size:.4f})",
+            flush=True,
+        )
+
+    summary = {
+        "protocol": (
+            "Blind: the hidden permutation is generated per trial and used "
+            "only for scoring. The spectral start is chosen by energy, never "
+            "by accuracy."
+        ),
+        "label": args.label,
+        "size": size,
+        "field_correlation_at_truth": rho,
+        "chance": 1.0 / size,
+        "rows": rows,
+        "mean_spectral": float(np.mean([r["spectral_accuracy"] for r in rows])),
+        "mean_refined": float(np.mean([r["refined_accuracy"] for r in rows])),
+    }
+    print(json.dumps({k: v for k, v in summary.items() if k != "rows"}))
+    write_json(args.output, summary)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/projection_sweep.py b/worldalign/projection_sweep.py
new file mode 100644
index 0000000..2326461
--- /dev/null
+++ b/worldalign/projection_sweep.py
@@ -0,0 +1,240 @@
+"""Sweep the content projection, the largest measured lever on natural data.
+
+Content projection took the field correlation from 0.489 to 0.656, further
+than any other single choice, and it was adopted at one setting without a
+sweep. This prices its free parameters -- how many directions to keep, how
+hard to shrink the within-scene scatter -- and two variants of the projection
+itself.
+
+The whitened variant rescales each kept direction by its own between-scene
+spread, so a direction that separates scenes weakly is not drowned by one
+that separates them strongly. The kernel variant fits the same discriminant
+in a random Fourier feature space, testing whether the directions that
+separate scenes are linear in the encoder's coordinates at all.
+
+Word vectors and object states are computed once and reused, so the whole
+sweep costs one pipeline run.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import re
+from pathlib import Path
+
+import numpy as np
+import torch
+import torch.nn.functional as F
+from scipy.linalg import eigh
+
+from .common import seed_everything, write_json
+from .natural_families import load_phrases
+from .natural_pipeline import word_vectors
+from .synth_cc_battery import moment_field
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--vg-dir", default="artifacts/vg_5k")
+    parser.add_argument("--objects", default="artifacts/vg_5k/natural_objects.pt")
+    parser.add_argument("--samples", type=int, default=256)
+    parser.add_argument("--fit-scenes", type=int, default=1200)
+    parser.add_argument("--word-vectors", type=int, default=48)
+    parser.add_argument("--min-count", type=int, default=60)
+    parser.add_argument("--max-vocabulary", type=int, default=600)
+    parser.add_argument("--context-window", type=int, default=2)
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--rff", type=int, default=256, help="Random Fourier features for the kernel variant.")
+    parser.add_argument("--output", default="artifacts/vg_5k/projection_sweep.json")
+    return parser.parse_args()
+
+
+def discriminant(
+    sets: list[np.ndarray], shrinkage: float
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Directions separating scenes more than parts, with their eigenvalues."""
+    means = np.stack([item.mean(0) for item in sets])
+    centre = means.mean(0)
+    within = np.concatenate([item - item.mean(0, keepdims=True) for item in sets])
+    scatter_within = within.T @ within / max(len(within) - 1, 1)
+    centred = means - centre
+    scatter_between = centred.T @ centred / max(len(centred) - 1, 1)
+    trace = np.trace(scatter_within) / len(scatter_within)
+    values, vectors = eigh(
+        scatter_between, scatter_within + shrinkage * trace * np.eye(len(scatter_within))
+    )
+    order = np.argsort(values)[::-1]
+    return centre, vectors[:, order], np.maximum(values[order], 0.0)
+
+
+def lift(raw: np.ndarray, weights: np.ndarray, offsets: np.ndarray) -> np.ndarray:
+    """Random Fourier features: an explicit map into an RBF feature space."""
+    return np.cos(raw @ weights + offsets) * np.sqrt(2.0 / weights.shape[1])
+
+
+def correlation(vision_sets, text_sets) -> float:
+    visual_field = moment_field(vision_sets)
+    text_field = moment_field(text_sets)
+    mask = ~np.eye(len(vision_sets), dtype=bool)
+    return float(
+        np.corrcoef(
+            visual_field.double().numpy()[mask], text_field.double().numpy()[mask]
+        )[0, 1]
+    )
+
+
+def main() -> None:
+    args = parse_args()
+    seed_everything(args.seed)
+    vg_dir = Path(args.vg_dir)
+    truth = [
+        json.loads(line)
+        for line in (vg_dir / "ground_truth.private.jsonl")
+        .read_text(encoding="utf-8")
+        .splitlines()
+        if line.strip()
+    ]
+    records = {
+        json.loads(line)["node_id"]: json.loads(line)
+        for line in (vg_dir / "text_nodes.jsonl").read_text(encoding="utf-8").splitlines()
+        if line.strip()
+    }
+    state = torch.load(args.objects, map_location="cpu", weights_only=False)
+    index = {node: position for position, node in enumerate(state["node_ids"])}
+    vectors = word_vectors(load_phrases(vg_dir / "text_nodes.jsonl", "region_closed"), args)
+
+    def language_state(node: str) -> np.ndarray | None:
+        rows = []
+        for phrase in records[node]["region_closed"]:
+            hits = [vectors[token] for token in re.findall(r"[a-z]+", phrase.lower())
+                    if token in vectors]
+            if hits:
+                rows.append(np.mean(hits, axis=0))
+        return np.stack(rows) if rows else None
+
+    def vision_state(node: str) -> np.ndarray | None:
+        segments = state["segments"][index[node]]
+        return (
+            np.stack([item["feature"] for item in segments]).astype(np.float64)
+            if segments
+            else None
+        )
+
+    pairs = [
+        pair for pair in truth
+        if pair["vision_node_id"] in index and pair["text_node_id"] in records
+    ]
+    fit = pairs[args.samples : args.samples + args.fit_scenes]
+    language_fit = [language_state(p["text_node_id"]) for p in fit]
+    language_fit = [i for i in language_fit if i is not None and len(i) > 1]
+    vision_fit = [vision_state(p["vision_node_id"]) for p in fit]
+    vision_fit = [i for i in vision_fit if i is not None and len(i) > 1]
+    evaluated = [
+        pair for pair in pairs[: args.samples]
+        if language_state(pair["text_node_id"]) is not None
+        and vision_state(pair["vision_node_id"]) is not None
+    ]
+    language_eval = [language_state(p["text_node_id"]) for p in evaluated]
+    vision_eval = [vision_state(p["vision_node_id"]) for p in evaluated]
+    print(f"fit {len(vision_fit)} vision / {len(language_fit)} text, "
+          f"evaluated {len(evaluated)}", flush=True)
+
+    rows = []
+
+    def record(variant: str, keep, shrinkage, value: float) -> None:
+        rows.append({"variant": variant, "keep": keep, "shrinkage": shrinkage,
+                     "correlation": value})
+        print(f"{variant:<10} keep={str(keep):<5} shrink={shrinkage:<6} rho={value:.4f}",
+              flush=True)
+
+    # --- linear discriminant: sweep kept width and shrinkage ---
+    for shrinkage in (0.01, 0.05, 0.2, 1.0):
+        lang = discriminant(language_fit, shrinkage)
+        vis = discriminant(vision_fit, shrinkage)
+        for keep in (4, 8, 16, 32, 64):
+            def project(raw, fitted):
+                centre, basis, _ = fitted
+                width = min(keep, basis.shape[1])
+                return F.normalize(
+                    torch.tensor((raw - centre) @ basis[:, :width], dtype=torch.float32), dim=-1
+                )
+            value = correlation(
+                [project(v, vis) for v in vision_eval],
+                [project(t, lang) for t in language_eval],
+            )
+            record("linear", keep, shrinkage, value)
+
+    # --- eigenvalue-weighted: scale each direction by its own separation ---
+    for shrinkage in (0.05, 0.2):
+        lang = discriminant(language_fit, shrinkage)
+        vis = discriminant(vision_fit, shrinkage)
+        for keep in (8, 16, 32, 64):
+            for power in (0.5, 1.0):
+                def project(raw, fitted):
+                    centre, basis, values = fitted
+                    width = min(keep, basis.shape[1])
+                    scale = values[:width] ** power
+                    return F.normalize(
+                        torch.tensor(
+                            (raw - centre) @ basis[:, :width] * scale, dtype=torch.float32
+                        ), dim=-1
+                    )
+                value = correlation(
+                    [project(v, vis) for v in vision_eval],
+                    [project(t, lang) for t in language_eval],
+                )
+                record(f"weighted^{power}", keep, shrinkage, value)
+
+    # --- kernel discriminant in a random Fourier feature space ---
+    generator = np.random.default_rng(args.seed)
+
+    def kernel_space(fit_sets, eval_sets, gamma_scale):
+        stacked = np.concatenate(fit_sets)
+        spread = np.median(np.linalg.norm(stacked - stacked.mean(0), axis=1))
+        gamma = gamma_scale / max(spread ** 2, 1e-9)
+        weights = generator.normal(
+            0.0, np.sqrt(2 * gamma), size=(stacked.shape[1], args.rff)
+        )
+        offsets = generator.uniform(0, 2 * np.pi, size=args.rff)
+        return (
+            [lift(item, weights, offsets) for item in fit_sets],
+            [lift(item, weights, offsets) for item in eval_sets],
+        )
+
+    for gamma_scale in (0.25, 1.0):
+        lang_fit_k, lang_eval_k = kernel_space(language_fit, language_eval, gamma_scale)
+        vis_fit_k, vis_eval_k = kernel_space(vision_fit, vision_eval, gamma_scale)
+        lang = discriminant(lang_fit_k, 0.05)
+        vis = discriminant(vis_fit_k, 0.05)
+        for keep in (8, 16, 32):
+            def project(raw, fitted):
+                centre, basis, _ = fitted
+                width = min(keep, basis.shape[1])
+                return F.normalize(
+                    torch.tensor((raw - centre) @ basis[:, :width], dtype=torch.float32), dim=-1
+                )
+            value = correlation(
+                [project(v, vis) for v in vis_eval_k],
+                [project(t, lang) for t in lang_eval_k],
+            )
+            record(f"kernel g={gamma_scale}", keep, 0.05, value)
+
+    best = max(rows, key=lambda row: row["correlation"])
+    summary = {
+        "protocol": (
+            "Projection fitted per modality on scenes outside the evaluated "
+            "set; hidden pairs read only for the correlation. Baseline is "
+            "linear, keep=8, shrinkage=0.05."
+        ),
+        "evaluated": len(evaluated),
+        "rows": rows,
+        "best": best,
+        "recovery_threshold": 0.9,
+    }
+    print(json.dumps(best))
+    write_json(args.output, summary)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/rank_causal.py b/worldalign/rank_causal.py
new file mode 100644
index 0000000..09ae55e
--- /dev/null
+++ b/worldalign/rank_causal.py
@@ -0,0 +1,141 @@
+"""Is the effective rank of the shared field the binding variable, not its magnitude?
+
+Stripping ten eigen-directions costs the synthetic fields that recover almost
+nothing -- 0.929 to 0.870 -- and costs the natural fields almost everything --
+0.656 to 0.080. The natural shared signal is therefore concentrated in about
+ten directions while the synthetic one is spread across the spectrum. That is
+a correlation between rank and recovery, and this makes it causal by
+intervention: take the synthetic fields that do recover, project them to
+successively lower rank, and watch recovery as a function of rank while the
+correlation between the two truncated fields stays high.
+
+If recovery dies at low rank while the correlation is still far above the
+threshold, then rank is the binding variable and the headline correlation is
+the wrong target. The theory agrees in advance: correlated-matching thresholds
+assume the agreement lives in exchangeable full-rank noise, so N^2/2 entries
+each carry an independent constraint. A rank-r shared component supplies about
+rN, and at 256 scenes with r near ten that is an eighth of what the threshold
+calculation assumes.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+
+from .common import write_json
+from .spectral_match import grampa
+from .synth_fast_gate import ClosedFormEnergy, all_swaps, steepest_descent
+from .synth_triangle_gate import standardized
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", default="artifacts/synth_v1/fields_tier0_ws_256.pt")
+    parser.add_argument("--label", default="synthetic-watershed")
+    parser.add_argument("--trials", type=int, default=3)
+    parser.add_argument("--device", default="cuda:3")
+    parser.add_argument("--output", default="artifacts/synth_v1/rank_causal.json")
+    return parser.parse_args()
+
+
+def offdiagonal(matrix: np.ndarray) -> np.ndarray:
+    return matrix[~np.eye(len(matrix), dtype=bool)]
+
+
+def normalise(matrix: np.ndarray) -> np.ndarray:
+    values = offdiagonal(matrix)
+    out = (matrix - values.mean()) / values.std()
+    np.fill_diagonal(out, 0.0)
+    return out
+
+
+def truncate(matrix: np.ndarray, rank: int) -> np.ndarray:
+    """Keep the leading `rank` eigen-components of the symmetrised field."""
+    symmetric = (matrix + matrix.T) / 2.0
+    values, vectors = np.linalg.eigh(symmetric)
+    order = np.argsort(np.abs(values))[::-1][:rank]
+    return (vectors[:, order] * values[order]) @ vectors[:, order].T
+
+
+def effective_rank(matrix: np.ndarray) -> float:
+    """Participation ratio of the eigen-spectrum: how many directions carry it."""
+    values = np.abs(np.linalg.eigvalsh((matrix + matrix.T) / 2.0))
+    weights = values / values.sum()
+    weights = weights[weights > 0]
+    return float(np.exp(-(weights * np.log(weights)).sum()))
+
+
+def main() -> None:
+    args = parse_args()
+    state = torch.load(args.fields, map_location="cpu", weights_only=False)
+    visual_full = state["visual_field"].double().numpy()
+    text_full = state["text_field"].double().numpy()
+    size = len(visual_full)
+    device = torch.device(args.device)
+    swaps = all_swaps(size, device)
+
+    rows = []
+    ranks = [4, 8, 16, 32, 64, 128, size]
+    for rank in ranks:
+        visual = normalise(truncate(visual_full, rank) if rank < size else visual_full)
+        text = normalise(truncate(text_full, rank) if rank < size else text_full)
+        mask = ~np.eye(size, dtype=bool)
+        rho = float(np.corrcoef(visual[mask], text[mask])[0, 1])
+
+        accuracies = []
+        for trial in range(args.trials):
+            generator = np.random.default_rng(trial)
+            hidden = generator.permutation(size)
+            shuffled = text[np.ix_(hidden, hidden)]
+            columns = grampa(visual, shuffled, 1.0)
+            energy = ClosedFormEnergy(
+                standardized(torch.from_numpy(shuffled).to(device)).float(),
+                standardized(torch.from_numpy(visual).to(device)).float(),
+                1.0, 1.0, 256,
+            )
+            final, _ = steepest_descent(
+                energy, torch.from_numpy(columns.copy()).to(device), swaps, 3000
+            )
+            accuracies.append(
+                float((hidden[final.cpu().numpy()] == np.arange(size)).mean())
+            )
+        row = {
+            "rank": rank,
+            "correlation": rho,
+            "recovery": float(np.mean(accuracies)),
+            "visual_effective_rank": effective_rank(visual),
+        }
+        rows.append(row)
+        print(
+            f"rank={rank:<5} rho={rho:.4f}  recovery={row['recovery']:.3f}  "
+            f"eff_rank={row['visual_effective_rank']:.1f}",
+            flush=True,
+        )
+
+    summary = {
+        "protocol": (
+            "Both fields truncated to the same rank, then matched blind from a "
+            "spectral start with exact refinement. Hidden pairing generated "
+            "per trial and read only for scoring."
+        ),
+        "label": args.label,
+        "size": size,
+        "chance": 1.0 / size,
+        "rows": rows,
+        "reading": (
+            "Recovery collapsing at low rank while the correlation stays above "
+            "the recovery threshold shows the headline correlation is not "
+            "sufficient and the shared spectrum's width is the binding "
+            "constraint."
+        ),
+    }
+    write_json(args.output, summary)
+    print(json.dumps({"chance": 1.0 / size}))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/shared_rank.py b/worldalign/shared_rank.py
new file mode 100644
index 0000000..629ae0b
--- /dev/null
+++ b/worldalign/shared_rank.py
@@ -0,0 +1,135 @@
+"""Which side is narrow? Per-modality rank against the shared rank.
+
+Truncating a recovering synthetic field to rank four leaves its correlation at
+0.90 and its recovery at 6%, so the width of the shared spectrum is a separate
+and binding constraint that the correlation does not express. The engineering
+question follows immediately: which modality is narrow. Upgrading the vision
+encoder is worth its compute only if vision is what limits the shared width,
+and if the text side is the narrow one no vision backbone will help.
+
+Three numbers per field pair. Each field's own effective rank, from the
+participation ratio of its spectrum, says how many directions the modality
+uses at all. The principal angles between the two leading eigenspaces say how
+many of those directions the two modalities share. And the per-direction
+correlation profile says where along the spectrum agreement dies.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+
+from .common import write_json
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--fields", nargs="+", required=True)
+    parser.add_argument("--labels", nargs="+", default=None)
+    parser.add_argument("--width", type=int, default=32)
+    parser.add_argument("--output", default="artifacts/vg_5k/shared_rank.json")
+    return parser.parse_args()
+
+
+def spectrum(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    symmetric = (matrix + matrix.T) / 2.0
+    values, vectors = np.linalg.eigh(symmetric)
+    order = np.argsort(np.abs(values))[::-1]
+    return values[order], vectors[:, order]
+
+
+def effective_rank(values: np.ndarray) -> float:
+    weights = np.abs(values) / np.abs(values).sum()
+    weights = weights[weights > 1e-15]
+    return float(np.exp(-(weights * np.log(weights)).sum()))
+
+
+def analyse(path: str, label: str, width: int) -> dict:
+    state = torch.load(path, map_location="cpu", weights_only=False)
+    visual = state["visual_field"].double().numpy()
+    text = state["text_field"].double().numpy()
+    size = len(visual)
+    mask = ~np.eye(size, dtype=bool)
+
+    visual_values, visual_vectors = spectrum(visual)
+    text_values, text_vectors = spectrum(text)
+    width = min(width, size)
+
+    # principal angles between the two leading eigenspaces: cosines near one
+    # count as directions the two modalities agree on. The raw count inflates
+    # with width -- two random w-dimensional subspaces of R^N already overlap
+    # at cosines near sqrt(w/N) -- so a scene-shuffled null is subtracted and
+    # every comparison is made at matched width.
+    overlap = visual_vectors[:, :width].T @ text_vectors[:, :width]
+    cosines = np.clip(np.linalg.svd(overlap, compute_uv=False), 0.0, 1.0)
+    shared = int((cosines > 0.7).sum())
+
+    null_counts = []
+    for trial in range(5):
+        generator = np.random.default_rng(trial)
+        order = generator.permutation(size)
+        shuffled = text[np.ix_(order, order)]
+        _, shuffled_vectors = spectrum(shuffled)
+        null_overlap = visual_vectors[:, :width].T @ shuffled_vectors[:, :width]
+        null_cosines = np.clip(np.linalg.svd(null_overlap, compute_uv=False), 0.0, 1.0)
+        null_counts.append(int((null_cosines > 0.7).sum()))
+    null = float(np.mean(null_counts))
+
+    # where along the spectrum does agreement die
+    profile = []
+    for index in range(min(width, 16)):
+        visual_direction = visual_vectors[:, index]
+        best = float(np.abs(text_vectors[:, :width].T @ visual_direction).max())
+        profile.append({"index": index, "best_text_overlap": best})
+
+    return {
+        "label": label,
+        "size": size,
+        "correlation": float(np.corrcoef(visual[mask], text[mask])[0, 1]),
+        "visual_effective_rank": effective_rank(visual_values),
+        "text_effective_rank": effective_rank(text_values),
+        "shared_directions_cos_gt_0.7": shared,
+        "shuffled_null": null,
+        "shared_above_null": shared - null,
+        "principal_cosines_top8": [round(float(c), 3) for c in cosines[:8]],
+        "spectrum_profile": profile,
+    }
+
+
+def main() -> None:
+    args = parse_args()
+    labels = args.labels or [path.split("/")[-1] for path in args.fields]
+    rows = [analyse(path, label, args.width)
+            for path, label in zip(args.fields, labels)]
+    for row in rows:
+        print(
+            f"{row['label']:<26} rho={row['correlation']:.3f}  "
+            f"eff_rank vision={row['visual_effective_rank']:6.1f} "
+            f"text={row['text_effective_rank']:6.1f}  "
+            f"shared={row['shared_directions_cos_gt_0.7']:3d} "
+            f"(null {row['shuffled_null']:4.1f}, net {row['shared_above_null']:5.1f})",
+            flush=True,
+        )
+    summary = {
+        "protocol": (
+            "Effective rank from the participation ratio of each field's "
+            "spectrum; shared directions from principal angles between the "
+            "two leading eigenspaces. No pairing enters beyond the alignment "
+            "the fields already carry."
+        ),
+        "width": args.width,
+        "rows": rows,
+        "reading": (
+            "The narrow modality is the one to upgrade. A shared count far "
+            "below both effective ranks means each side is rich but they are "
+            "rich about different things."
+        ),
+    }
+    write_json(args.output, summary)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/worldalign/tier0_pipeline.py b/worldalign/tier0_pipeline.py
index 778e20f..01c78d2 100644
--- a/worldalign/tier0_pipeline.py
+++ b/worldalign/tier0_pipeline.py
@@ -56,6 +56,12 @@ def parse_args() -> argparse.Namespace:
     parser.add_argument("--seed", type=int, default=0)
     parser.add_argument("--output", required=True)
     parser.add_argument("--states-output", default="")
+    parser.add_argument(
+        "--text-omit",
+        default="",
+        help="Comma-separated factors the captions never state "
+        "(colour, count, size). Vision still sees them.",
+    )
     return parser.parse_args()
 
 
@@ -188,7 +194,17 @@ def factor_vector(colour: int, count: int, size: int, classes: int) -> torch.Ten
     return vector
 
 
-def encode_caption(caption: str, colour_rank: dict[str, int], classes: int) -> torch.Tensor:
+def encode_caption(
+    caption: str, colour_rank: dict[str, int], classes: int, omit: frozenset = frozenset()
+) -> torch.Tensor:
+    """Encode a caption, optionally suppressing factors the text never states.
+
+    Suppressing a factor makes the caption describe strictly less of the scene
+    while vision continues to see all of it, which is the controlled form of
+    the situation on photographs: a region description and a patch descriptor
+    overlap on some aspects and not others. It is the intervention that tests
+    whether corpus overlap is what sets the width of the shared spectrum.
+    """
     vectors = []
     for phrase in parse_group_phrases(caption):
         tokens = phrase.split()
@@ -201,10 +217,20 @@ def encode_caption(caption: str, colour_rank: dict[str, int], classes: int) -> t
             )
         )
         size = 0 if "small" in tokens else (2 if "large" in tokens else 1)
-        vectors.append(factor_vector(colour, count, size, classes))
+        vector = factor_vector(colour, count, size, classes)
+        if "colour" in omit:
+            vector[:classes] = 0.0
+        if "count" in omit:
+            vector[classes : classes + 4] = 0.0
+        if "size" in omit:
+            vector[classes + 4 :] = 0.0
+        vectors.append(vector)
     if not vectors:
         vectors = [factor_vector(0, 1, 1, classes)]
-    return F.normalize(torch.stack(vectors), dim=-1)
+    stacked = torch.stack(vectors)
+    if stacked.abs().sum() == 0:
+        stacked = stacked + 1.0
+    return F.normalize(stacked, dim=-1)
 
 
 def moment_state(states: torch.Tensor) -> torch.Tensor:
@@ -264,7 +290,10 @@ def main() -> None:
             view_states = [moment_state(item) for item in sets]
     visual_field = torch.stack(per_view_fields).mean(0)
 
-    text_sets = [encode_caption(captions[row][0], colour_rank, classes) for row in rows]
+    omit = frozenset(f for f in args.text_omit.split(",") if f)
+    text_sets = [
+        encode_caption(captions[row][0], colour_rank, classes, omit) for row in rows
+    ]
     text_field = moment_field(text_sets)
 
     mask = ~np.eye(len(rows), dtype=bool)