diff options
Diffstat (limited to 'worldalign')
| -rw-r--r-- | worldalign/degeneracy_gate.py | 130 | ||||
| -rw-r--r-- | worldalign/field_geometry.py | 222 | ||||
| -rw-r--r-- | worldalign/match_battery.py | 205 | ||||
| -rw-r--r-- | worldalign/split_half.py | 190 | ||||
| -rw-r--r-- | worldalign/synth_fast_gate.py | 51 |
5 files changed, 798 insertions, 0 deletions
diff --git a/worldalign/degeneracy_gate.py b/worldalign/degeneracy_gate.py new file mode 100644 index 0000000..4ba3dc1 --- /dev/null +++ b/worldalign/degeneracy_gate.py @@ -0,0 +1,130 @@ +"""Is the truth the optimum, tied with it, or beaten? Three descents, not seven hundred. + +The question that decides the next move is narrow, and the first attempt at it +was not. Suppressing a caption factor leaves the information intact -- the +anchor bound is 0.997 -- while blind matching returns 5.6%, and there are three +possible reasons with three different responses: + + E(truth) < E(best found) the searcher simply missed it -> better optimiser + E(truth) == E(best found) the truth is one of many minima -> quotient the symmetry + E(truth) > E(best found) the truth is not the optimum -> information is gone + +Answering it does not need a strong searcher. It needs the energy of the truth +and the energy of whatever the standard pipeline actually converges to, which +is one descent per trial rather than the sixty-one the first version ran -- at +N=256 a single descent over all 32,640 transpositions costs minutes, so the +difference is hours against seconds. +""" + +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, fast_pair_descent +from .synth_triangle_gate import standardized + + +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("--trials", type=int, default=3) + parser.add_argument("--iterations", type=int, default=600) + parser.add_argument("--device", default="cuda:3") + parser.add_argument("--tolerance", type=float, default=1e-3) + parser.add_argument("--output", default="artifacts/synth_v1/degeneracy_gate.json") + return parser.parse_args() + + +def standardise(matrix: np.ndarray) -> np.ndarray: + mask = ~np.eye(len(matrix), dtype=bool) + values = matrix[mask] + out = (matrix - values.mean()) / values.std() + np.fill_diagonal(out, 0.0) + return out + + +def analyse(path: str, label: str, args: argparse.Namespace) -> dict: + state = torch.load(path, map_location="cpu", weights_only=False) + visual = standardise(state["visual_field"].double().numpy()) + text = standardise(state["text_field"].double().numpy()) + size = len(visual) + device = torch.device(args.device) + + gaps, accuracies, truths = [], [], [] + for trial in range(args.trials): + generator = np.random.default_rng(trial) + hidden = generator.permutation(size) + shuffled = text[np.ix_(hidden, hidden)] + shuffled_gpu = standardized(torch.from_numpy(shuffled).to(device)).double() + visual_gpu = standardized(torch.from_numpy(visual).to(device)).double() + energy = ClosedFormEnergy(shuffled_gpu, visual_gpu, 1.0, 0.0, 256) + truth = torch.from_numpy(np.argsort(hidden).copy()).to(device) + truth_energy = float(energy.energy(truth[None])[0]) + + start = torch.from_numpy(grampa(visual, shuffled, 1.0).copy()).to(device) + final = fast_pair_descent(shuffled_gpu, visual_gpu, start, args.iterations) + found_energy = float(energy.energy(final[None])[0]) + # descent from the truth itself: if it moves, the truth is not a local + # minimum either, which separates "not optimal" from "merely not found" + from_truth = fast_pair_descent( + shuffled_gpu, visual_gpu, truth.clone(), args.iterations + ) + truth_local = float(energy.energy(from_truth[None])[0]) + + gaps.append(found_energy - truth_energy) + truths.append(truth_energy - truth_local) + accuracies.append(float((hidden[final.cpu().numpy()] == np.arange(size)).mean())) + + gap = float(np.mean(gaps)) + scale = abs(float(np.mean([g for g in gaps]))) + 1.0 + if gap < -args.tolerance: + verdict = "truth is NOT the optimum -- information deficit" + elif abs(gap) <= args.tolerance: + verdict = "truth is TIED with what is found -- degenerate, quotient the symmetry" + else: + verdict = "truth is DEEPER than what is found -- optimiser failure" + return { + "label": label, + "size": size, + "energy_gap_found_minus_truth": gap, + "truth_descends_further_by": float(np.mean(truths)), + "blind_accuracy": float(np.mean(accuracies)), + "verdict": verdict, + } + + +def main() -> None: + args = parse_args() + labels = args.labels or [p.split("/")[-1] for p in args.fields] + rows = [] + for path, label in zip(args.fields, labels): + row = analyse(path, label, args) + rows.append(row) + print( + f"{row['label']:<20} gap(found-truth)={row['energy_gap_found_minus_truth']:+.4f} " + f"truth-descends={row['truth_descends_further_by']:+.4f} " + f"acc={row['blind_accuracy']:.3f}\n" + f"{'':<20} -> {row['verdict']}", + flush=True, + ) + write_json(args.output, { + "protocol": ( + "One spectral-start descent per trial, plus one descent started at " + "the truth. A positive gap means the searcher stopped above the " + "truth; zero means it found something the energy cannot " + "distinguish from the truth; negative means the truth is beaten." + ), + "rows": rows, + }) + print(json.dumps({"done": True})) + + +if __name__ == "__main__": + main() diff --git a/worldalign/field_geometry.py b/worldalign/field_geometry.py new file mode 100644 index 0000000..9360427 --- /dev/null +++ b/worldalign/field_geometry.py @@ -0,0 +1,222 @@ +"""The space and the operators, which we never varied. + +Everything measured so far fixed one geometry without saying so. Scene states +are compared by Euclidean inner product, the relation field holds those raw +values, and two fields are compared by Pearson correlation -- three +commitments to flatness made by default rather than by argument. The +measured deficit is that the two modalities disagree about which scenes are +similar, and "similar" is exactly what those three choices define. + +The mismatch that costs the most is likely the cheapest to remove. Two +independently trained encoders have no reason to put their similarities on +the same scale: if vision similarity is any monotone but nonlinear function +of text similarity, the two fields describe identical relational structure +and Pearson correlation still reports disagreement. A rank transform is +invariant to every monotone reparametrisation and costs one sort. + +Density is the second. A scene sitting in a crowded region of one modality +and a sparse region of the other has systematically shifted similarities to +everything, which local scaling and the diffusion kernel both remove -- the +latter by replacing similarity with how probability spreads, a quantity +defined by the graph rather than by whatever units the encoder happened to +emit. + +Each transform is applied to each modality independently, so nothing here +uses the pairing or crosses the gauge. +""" + +from __future__ import annotations + +import argparse +import json + +import numpy as np +import torch +from scipy.optimize import linear_sum_assignment +from scipy.stats import rankdata + +from .common import write_json + + +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("--repeats", type=int, default=5) + parser.add_argument("--output", default="artifacts/vg_5k/field_geometry.json") + return parser.parse_args() + + +def offdiagonal_mask(size: int) -> np.ndarray: + return ~np.eye(size, dtype=bool) + + +def standardise(matrix: np.ndarray) -> np.ndarray: + mask = offdiagonal_mask(len(matrix)) + values = matrix[mask] + out = (matrix - values.mean()) / values.std() + np.fill_diagonal(out, 0.0) + return out + + +# --- geometries, each applied within one modality ------------------------- + +def identity(field: np.ndarray) -> np.ndarray: + return standardise(field) + + +def global_rank(field: np.ndarray) -> np.ndarray: + """Copula transform: keep the ordering of similarities, discard the scale.""" + size = len(field) + mask = offdiagonal_mask(size) + out = np.zeros_like(field) + out[mask] = rankdata(field[mask]) / mask.sum() + return standardise(out) + + +def row_rank(field: np.ndarray) -> np.ndarray: + """Rank within each scene's own neighbourhood, removing per-scene scale.""" + size = len(field) + out = np.zeros_like(field) + for row in range(size): + others = np.delete(np.arange(size), row) + out[row, others] = rankdata(field[row, others]) / len(others) + out = (out + out.T) / 2.0 + return standardise(out) + + +def local_scaling(field: np.ndarray, neighbours: int = 10) -> np.ndarray: + """Self-tuning: divide by each scene's own neighbourhood scale (CSLS-like).""" + size = len(field) + work = field.copy() + np.fill_diagonal(work, -np.inf) + order = np.sort(work, axis=1)[:, ::-1] + scale = order[:, :neighbours].mean(1) + out = field - 0.5 * (scale[:, None] + scale[None, :]) + np.fill_diagonal(out, 0.0) + return standardise(out) + + +def diffusion(field: np.ndarray, steps: int = 3) -> np.ndarray: + """Replace similarity by how probability spreads on the similarity graph.""" + size = len(field) + affinity = field - field.min() + np.fill_diagonal(affinity, 0.0) + degree = affinity.sum(1).clip(1e-9) + walk = affinity / degree[:, None] + powered = np.linalg.matrix_power(walk, steps) + out = powered @ powered.T + np.fill_diagonal(out, 0.0) + return standardise(out) + + +def heat_kernel(field: np.ndarray, time: float = 1.0) -> np.ndarray: + """Heat flow on the normalised Laplacian: a metric the graph defines itself.""" + size = len(field) + affinity = field - field.min() + np.fill_diagonal(affinity, 0.0) + degree = affinity.sum(1).clip(1e-9) + normalised = affinity / np.sqrt(np.outer(degree, degree)) + values, vectors = np.linalg.eigh((normalised + normalised.T) / 2.0) + out = (vectors * np.exp(time * (values - values.max()))) @ vectors.T + np.fill_diagonal(out, 0.0) + return standardise(out) + + +def hyperbolic(field: np.ndarray, dimension: int = 16) -> np.ndarray: + """Embed each field in hyperbolic space and re-derive distances there. + + If scene similarity is hierarchical -- broad categories containing finer + ones -- a flat inner product cannot hold it without distortion, while + negative curvature can. The field is embedded by eigenmap and read back + as Lorentzian distance. + """ + size = len(field) + values, vectors = np.linalg.eigh((field + field.T) / 2.0) + order = np.argsort(values)[::-1][:dimension] + coordinates = vectors[:, order] * np.sqrt(np.abs(values[order])) + scale = np.abs(coordinates).max().clip(1e-9) + spatial = coordinates / scale * 0.9 + time = np.sqrt(1.0 + (spatial ** 2).sum(1)) + product = np.outer(time, time) - spatial @ spatial.T + out = -np.arccosh(np.clip(product, 1.0, None)) + np.fill_diagonal(out, 0.0) + return standardise(out) + + +GEOMETRIES = { + "euclidean (current)": identity, + "global rank (copula)": global_rank, + "row rank": row_rank, + "local scaling": local_scaling, + "diffusion t=3": diffusion, + "heat kernel": heat_kernel, + "hyperbolic d=16": hyperbolic, +} + + +# --- comparison operators ------------------------------------------------- + +def anchor_bound(first: np.ndarray, second: np.ndarray, repeats: int, + spearman: bool = False) -> float: + size = len(first) + scores = [] + for repeat in range(repeats): + generator = np.random.default_rng(repeat) + order = generator.permutation(size) + half = size // 2 + anchors, probe = order[:half], order[half:] + + def profile(field): + block = field[np.ix_(probe, anchors)] + if spearman: + block = np.apply_along_axis(rankdata, 1, block) + block = block - block.mean(1, keepdims=True) + return block / block.std(1, keepdims=True).clip(1e-9) + + similarity = profile(first) @ profile(second).T / half + _, columns = linear_sum_assignment(-similarity) + scores.append(float((columns == np.arange(len(probe))).mean())) + return float(np.mean(scores)) + + +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() + mask = offdiagonal_mask(len(visual)) + + rows = [] + print("%-24s %8s %10s %12s" % ("geometry", "rho", "bound", "bound+rank-op")) + for name, transform in GEOMETRIES.items(): + try: + first, second = transform(visual), transform(text) + except Exception as error: # a geometry that cannot be formed is a result + print("%-24s failed: %s" % (name, error)) + continue + rho = float(np.corrcoef(first[mask], second[mask])[0, 1]) + plain = anchor_bound(first, second, args.repeats) + ranked = anchor_bound(first, second, args.repeats, spearman=True) + rows.append({"geometry": name, "correlation": rho, + "anchor_bound": plain, "anchor_bound_rank_operator": ranked}) + print("%-24s %8.3f %10.4f %12.4f" % (name, rho, plain, ranked)) + + best = max(rows, key=lambda row: max(row["anchor_bound"], + row["anchor_bound_rank_operator"])) + write_json(args.output, { + "protocol": ( + "Each geometry is applied to each modality independently, so no " + "transform crosses the gauge or uses the pairing. The bound is the " + "anchor upper bound; the last column re-computes it with a rank " + "comparison operator instead of a linear one." + ), + "label": args.label, + "rows": rows, + "best": best, + }) + print(json.dumps(best)) + + +if __name__ == "__main__": + main() diff --git a/worldalign/match_battery.py b/worldalign/match_battery.py new file mode 100644 index 0000000..45fa26e --- /dev/null +++ b/worldalign/match_battery.py @@ -0,0 +1,205 @@ +"""A matching benchmark, and the solvers we never tried on it. + +Everything the project has used to search is one narrow family: spectral +initialisation on the pairwise field, exact steepest descent over +transpositions, tempering, and Sinkhorn. All of them act on the permutation +group with a pairwise objective in flat space. The measured cost of that +narrowness is large and specific -- on a field whose anchor bound is 0.997, +meaning the information is intact almost perfectly, blind matching returns +5.6%. + +That makes a clean benchmark, because unlike natural data the answer is known +to be reachable. Three synthetic fields span the range: one that already +recovers, one whose captions omit a factor, one truncated to rank eight. +Natural data is included as the case where no method should be expected to +work, so a method that scores there is reporting a bug rather than a result. + +Gromov-Wasserstein is the canonical formulation for exactly our problem -- +align two metric spaces with no shared embedding, by matching their internal +distance structures -- and it has never been run on these fields. Its entropic +form with an annealed regulariser is the landscape-deformation method that the +rank ladder argued for and that we never reached. Each solver is scored blind +and then re-scored after exact refinement, since composition is what made the +spectral solver work. +""" + +from __future__ import annotations + +import argparse +import json +import time + +import numpy as np +import torch + +from .common import write_json +from .spectral_match import grampa, umeyama +from .synth_fast_gate import ClosedFormEnergy, fast_pair_descent +from .synth_triangle_gate import standardized + + +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("--trials", type=int, default=3) + parser.add_argument("--iterations", type=int, default=600) + parser.add_argument("--device", default="cuda:3") + parser.add_argument("--methods", nargs="+", default=None) + parser.add_argument("--output", default="artifacts/synth_v1/match_battery.json") + return parser.parse_args() + + +def standardise(matrix: np.ndarray) -> np.ndarray: + mask = ~np.eye(len(matrix), dtype=bool) + values = matrix[mask] + out = (matrix - values.mean()) / values.std() + np.fill_diagonal(out, 0.0) + return out + + +def coupling_to_permutation(coupling: np.ndarray) -> np.ndarray: + from scipy.optimize import linear_sum_assignment + + _, columns = linear_sum_assignment(-coupling) + return columns + + +# --- solvers: each returns column assignment for rows of `visual` --------- + +def solve_grampa(visual, shuffled, **_): + return grampa(visual, shuffled, 1.0) + + +def solve_umeyama(visual, shuffled, **_): + return umeyama(visual, shuffled) + + +def solve_gw(visual, shuffled, **_): + import ot + + size = len(visual) + weights = ot.unif(size) + coupling = ot.gromov.gromov_wasserstein( + visual, shuffled, weights, weights, "square_loss", + max_iter=200, tol_rel=1e-9, + ) + return coupling_to_permutation(coupling) + + +def solve_gw_annealed(visual, shuffled, **_): + """Entropic GW with the regulariser annealed from smooth to sharp.""" + import ot + + size = len(visual) + weights = ot.unif(size) + coupling = None + for epsilon in (5e-2, 2e-2, 1e-2, 5e-3, 2e-3): + try: + coupling = ot.gromov.entropic_gromov_wasserstein( + visual, shuffled, weights, weights, "square_loss", + epsilon=epsilon, G0=coupling, max_iter=200, tol=1e-9, + solver="PGD", + ) + except Exception: + break + if not np.isfinite(coupling).all(): + break + if coupling is None or not np.isfinite(coupling).all(): + raise RuntimeError("entropic GW diverged at every regulariser") + return coupling_to_permutation(coupling) + + +def solve_bapg(visual, shuffled, **_): + """Bregman alternating projected gradient: a newer GW solver.""" + import ot + + size = len(visual) + weights = ot.unif(size) + coupling = ot.gromov.BAPG_gromov_wasserstein( + visual, shuffled, weights, weights, "square_loss", + epsilon=1e-1, max_iter=1000, tol=1e-9, + ) + return coupling_to_permutation(coupling) + + +SOLVERS = { + "grampa (current)": solve_grampa, + "umeyama": solve_umeyama, + "gromov-wasserstein": solve_gw, + "entropic GW annealed": solve_gw_annealed, + "BAPG GW": solve_bapg, +} + + +def main() -> None: + args = parse_args() + labels = args.labels or [p.split("/")[-1] for p in args.fields] + chosen = args.methods or list(SOLVERS) + device = torch.device(args.device) + rows = [] + + for path, label in zip(args.fields, labels): + state = torch.load(path, map_location="cpu", weights_only=False) + visual = standardise(state["visual_field"].double().numpy()) + text = standardise(state["text_field"].double().numpy()) + size = len(visual) + text_gpu = standardized(torch.from_numpy(text).to(device)).double() + visual_gpu = standardized(torch.from_numpy(visual).to(device)).double() + print(f"\n=== {label} (N={size}, chance={1/size:.4f})", flush=True) + + for name in chosen: + solver = SOLVERS[name] + raw, refined, seconds, failures = [], [], [], [] + for trial in range(args.trials): + generator = np.random.default_rng(trial) + hidden = generator.permutation(size) + shuffled = text[np.ix_(hidden, hidden)] + start = time.time() + try: + columns = solver(visual=visual, shuffled=shuffled) + except Exception as error: + failures.append(str(error)[:80]) + continue + seconds.append(time.time() - start) + raw.append(float((hidden[columns] == np.arange(size)).mean())) + + shuffled_gpu = standardized( + torch.from_numpy(shuffled).to(device) + ).double() + final = fast_pair_descent( + shuffled_gpu, visual_gpu, + torch.from_numpy(np.ascontiguousarray(columns)).to(device), + args.iterations, + ) + refined.append( + float((hidden[final.cpu().numpy()] == np.arange(size)).mean()) + ) + row = { + "field": label, "method": name, + "raw": float(np.mean(raw)) if raw else None, + "refined": float(np.mean(refined)) if refined else None, + "seconds": float(np.mean(seconds)) if seconds else None, + "failures": failures, + } + rows.append(row) + if raw: + print(" %-22s raw=%.4f refined=%.4f (%.1fs)" + % (name, row["raw"], row["refined"], row["seconds"]), flush=True) + else: + print(" %-22s FAILED: %s" % (name, failures[0] if failures else "?"), + flush=True) + + write_json(args.output, { + "protocol": ( + "Blind: the hidden permutation is generated per trial and read only " + "for scoring. Each solver is scored as returned and again after " + "exact steepest-descent refinement." + ), + "rows": rows, + }) + print(json.dumps({"done": True})) + + +if __name__ == "__main__": + main() diff --git a/worldalign/split_half.py b/worldalign/split_half.py new file mode 100644 index 0000000..ee3f273 --- /dev/null +++ b/worldalign/split_half.py @@ -0,0 +1,190 @@ +"""The ceiling nobody measured: can a modality identify a scene against itself? + +Every deficit so far has been read as a cross-modal one, and that assumed each +side describes its scenes well enough to identify them. Nothing tested the +assumption. Split each scene's parts into two disjoint halves -- odd and even +regions, odd and even segments -- build a relation field from each half, and +match the halves against each other with the anchor bound. Both halves come +from the same modality and the same scene, so any failure is that modality +failing to pin down its own scenes. + +This upper-bounds anything cross-modal. If one half of an image's description +cannot pick that image out from among the others given the other half, no +vision encoder can be blamed for failing to, and no better cross-modal recipe +exists to find. It costs one pass over data already on disk, and it should have +been the first measurement on natural data rather than one of the last. +""" + +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.optimize import linear_sum_assignment + +from .common import seed_everything, write_json +from .natural_families import load_phrases +from .natural_pipeline import content_directions, 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/nobj_base_gpu.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=128) + 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=64) + parser.add_argument("--shrinkage", type=float, default=0.05) + parser.add_argument("--weight-power", type=float, default=0.5) + parser.add_argument("--seed", type=int, default=0) + parser.add_argument("--output", default="artifacts/vg_5k/split_half.json") + return parser.parse_args() + + +def standardise(matrix: np.ndarray) -> np.ndarray: + values = matrix[~np.eye(len(matrix), dtype=bool)] + out = (matrix - values.mean()) / values.std() + np.fill_diagonal(out, 0.0) + return out + + +def anchor_bound(first: np.ndarray, second: np.ndarray, repeats: int = 5) -> float: + size = len(first) + scores = [] + for repeat in range(repeats): + generator = np.random.default_rng(repeat) + order = generator.permutation(size) + half = size // 2 + anchors, probe = order[:half], order[half:] + + def profile(field): + block = field[np.ix_(probe, anchors)] + block = block - block.mean(1, keepdims=True) + return block / block.std(1, keepdims=True).clip(1e-9) + + similarity = profile(first) @ profile(second).T / half + _, columns = linear_sum_assignment(-similarity) + scores.append(float((columns == np.arange(len(probe))).mean())) + return float(np.mean(scores)) + + +def build_field(sets, fitted, keep, power) -> torch.Tensor: + centre, basis, values = fitted + width = min(keep, basis.shape[1]) + scale = values[:width] ** power if power else 1.0 + projected = [ + F.normalize( + torch.tensor((item - centre) @ basis[:, :width] * scale, dtype=torch.float32), + dim=-1, + ) + for item in sets + ] + return moment_field(projected) + + +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 text_parts(node): + rows = [] + for phrase in records[node]["region_closed"]: + hits = [vectors[t] for t in re.findall(r"[a-z]+", phrase.lower()) if t in vectors] + if hits: + rows.append(np.mean(hits, axis=0)) + return rows + + def vision_parts(node): + return [s["feature"].astype(np.float64) for s in state["segments"][index[node]]] + + pairs = [p for p in truth if p["vision_node_id"] in index and p["text_node_id"] in records] + + def collect(getter, key, rows, minimum): + out = [] + for pair in rows: + parts = getter(pair[key]) + out.append(np.stack(parts) if len(parts) >= minimum else None) + return out + + results = {} + for name, getter, key in ( + ("text", text_parts, "text_node_id"), + ("vision", vision_parts, "vision_node_id"), + ): + fit_rows = pairs[args.samples : args.samples + args.fit_scenes] + fit = [x for x in collect(getter, key, fit_rows, 2) if x is not None] + fitted = content_directions(fit, args.shrinkage) + + evaluated, first_half, second_half = [], [], [] + for pair in pairs[: args.samples]: + parts = getter(pair[key]) + if len(parts) < 4: + continue + stacked = np.stack(parts) + evaluated.append(pair) + first_half.append(stacked[0::2]) + second_half.append(stacked[1::2]) + + field_a = standardise( + build_field(first_half, fitted, args.keep, args.weight_power).double().numpy() + ) + field_b = standardise( + build_field(second_half, fitted, args.keep, args.weight_power).double().numpy() + ) + mask = ~np.eye(len(evaluated), dtype=bool) + results[name] = { + "scenes": len(evaluated), + "mean_parts": float(np.mean([len(getter(p[key])) for p in evaluated])), + "split_half_correlation": float( + np.corrcoef(field_a[mask], field_b[mask])[0, 1] + ), + "split_half_anchor_bound": anchor_bound(field_a, field_b), + } + print( + f"{name:<7} scenes={results[name]['scenes']:4d} " + f"parts={results[name]['mean_parts']:5.1f} " + f"rho={results[name]['split_half_correlation']:.3f} " + f"self-identification bound={results[name]['split_half_anchor_bound']:.3f}", + flush=True, + ) + + summary = { + "protocol": ( + "Each scene's parts split into disjoint halves within one modality; " + "a relation field is built from each half and the two are matched by " + "the anchor bound. Both halves describe the same scenes, so this " + "upper-bounds any cross-modal result." + ), + "results": results, + "cross_modal_bound_for_reference": 0.291, + } + print(json.dumps(summary["results"])) + write_json(args.output, summary) + + +if __name__ == "__main__": + main() diff --git a/worldalign/synth_fast_gate.py b/worldalign/synth_fast_gate.py index 54b6891..01fb18f 100644 --- a/worldalign/synth_fast_gate.py +++ b/worldalign/synth_fast_gate.py @@ -343,3 +343,54 @@ def main() -> None: if __name__ == "__main__": main() + + +def all_pair_swap_deltas( + permuted_text: torch.Tensor, visual: torch.Tensor +) -> torch.Tensor: + """Energy change of every transposition at once, for the pairwise term. + + The brute-force descent forms all N(N-1)/2 candidate permutations and + scores each through the full energy, which at 256 scenes means 32,640 + matrices per step and, with the triangle term active, a batched cube + trace over every one of them -- 203 seconds for a single descent. The + pairwise term does not need any of that. Writing S for the alignment sum + and C for `permuted_text @ visual`, swapping positions p and q changes S + by 2(C_pq + C_qp - C_pp - C_qq + 2 A_pq B_pq), so one matrix product + yields the entire table. + + Both fields must be symmetric with zero diagonal, which standardisation + already guarantees. Returns a dense [N, N] table of energy deltas; the + sign convention matches `ClosedFormEnergy.energy`, so negative is an + improvement. + """ + size = permuted_text.shape[-1] + product = permuted_text @ visual + diagonal = torch.diagonal(product) + gain = 2.0 * ( + product + product.T + - diagonal[:, None] - diagonal[None, :] + + 2.0 * permuted_text * visual + ) + return -2.0 * gain / (size * (size - 1)) + + +def fast_pair_descent( + text: torch.Tensor, visual: torch.Tensor, start: torch.Tensor, max_steps: int +) -> torch.Tensor: + """Steepest descent on the pairwise energy using the closed-form table.""" + size = len(visual) + current = start.clone() + triangle = torch.triu( + torch.ones(size, size, dtype=torch.bool, device=visual.device), diagonal=1 + ) + for _ in range(max_steps): + permuted = text[current[:, None], current[None, :]] + deltas = all_pair_swap_deltas(permuted, visual) + deltas = torch.where(triangle, deltas, torch.inf) + best = int(deltas.argmin()) + if float(deltas.view(-1)[best]) >= -1e-12: + break + p, q = best // size, best % size + current[[p, q]] = current[[q, p]] + return current |