summaryrefslogtreecommitdiff
path: root/worldalign
diff options
context:
space:
mode:
Diffstat (limited to 'worldalign')
-rw-r--r--worldalign/anchor_bound.py14
-rw-r--r--worldalign/automorphism_probe.py96
-rw-r--r--worldalign/fused_match.py214
-rw-r--r--worldalign/match_battery.py18
4 files changed, 333 insertions, 9 deletions
diff --git a/worldalign/anchor_bound.py b/worldalign/anchor_bound.py
index 49f4ed3..2920934 100644
--- a/worldalign/anchor_bound.py
+++ b/worldalign/anchor_bound.py
@@ -79,13 +79,19 @@ def analyse(path: str, label: str, repeats: int, blind: float | None) -> dict:
order = generator.permutation(size)
half = size // 2
anchors, probe = order[:half], order[half:]
+ # The two sides must not be presented in the same index order. Exact
+ # twins have identical anchor-restricted rows, and with matched
+ # ordering `linear_sum_assignment` breaks the tie onto the diagonal --
+ # crediting the bound with information it does not have. On the
+ # caption-omitted field that inflated it from 0.920 to 0.997.
+ order = generator.permutation(len(probe))
left = rows_against(visual, probe, anchors)
- right = rows_against(text, probe, anchors)
+ right = rows_against(text, probe[order], anchors)
similarity = left @ right.T / left.shape[1]
- truth = np.arange(len(probe))
- nearest.append(float((similarity.argmax(1) == truth).mean()))
+ truth = order
+ nearest.append(float((order[similarity.argmax(1)] == np.arange(len(probe))).mean()))
_, columns = linear_sum_assignment(-similarity)
- assigned.append(float((columns == truth).mean()))
+ assigned.append(float((order[columns] == np.arange(len(probe))).mean()))
row = {
"label": label,
diff --git a/worldalign/automorphism_probe.py b/worldalign/automorphism_probe.py
new file mode 100644
index 0000000..46949f4
--- /dev/null
+++ b/worldalign/automorphism_probe.py
@@ -0,0 +1,96 @@
+"""Exact automorphisms of a relation field, and the blind-recovery ceiling.
+
+If sigma is an exact automorphism of the text field T (P_sigma T P_sigma^T = T
+bit-for-bit) then the blind matching problem cannot distinguish the truth from
+truth o sigma: the input (V, T) is literally unchanged. Every functional of
+(V, P T P^T) -- the pairwise energy, the third-order tr(M^3) term, any moment
+in the ladder, the TLB linear term -- is invariant, so no optimiser and no
+relaxation can break the tie. The posterior over the truth is uniform on the
+orbit, and the expected accuracy of *any* blind estimator is bounded by
+
+ ceiling = (1/n) * sum_over_orbits |orbit| * (1/|orbit|)
+ = 1 - (moved - #orbits) / n
+
+Run: python -m worldalign.automorphism_probe artifacts/synth_v1/omit_size.pt
+"""
+
+from __future__ import annotations
+
+import math
+import sys
+
+import numpy as np
+import torch
+
+
+def standardise(matrix: np.ndarray) -> np.ndarray:
+ matrix = np.asarray(matrix, dtype=np.float64)
+ mask = ~np.eye(len(matrix), dtype=bool)
+ out = (matrix - matrix[mask].mean()) / matrix[mask].std()
+ np.fill_diagonal(out, 0.0)
+ return out
+
+
+def automorphic_transpositions(field: np.ndarray, tol: float = 1e-9) -> np.ndarray:
+ """adj[i, j] iff swapping i and j leaves the field exactly invariant.
+
+ Rows i and j must agree on every coordinate outside {i, j}; the two
+ excluded coordinates are exactly the ones the swap moves.
+ """
+ size = len(field)
+ mismatch = np.full((size, size), np.inf)
+ for i in range(size):
+ gap = np.abs(field - field[i])
+ gap[:, i] = 0.0
+ gap[np.arange(size), np.arange(size)] = 0.0
+ mismatch[i] = gap.max(axis=1)
+ np.fill_diagonal(mismatch, np.inf)
+ return mismatch < tol
+
+
+def orbits(adjacency: np.ndarray) -> list[list[int]]:
+ size = len(adjacency)
+ seen = np.zeros(size, dtype=bool)
+ found: list[list[int]] = []
+ for start in range(size):
+ if seen[start]:
+ continue
+ stack = [start]
+ seen[start] = True
+ component = [start]
+ while stack:
+ node = stack.pop()
+ for nxt in np.nonzero(adjacency[node] & ~seen)[0]:
+ seen[nxt] = True
+ stack.append(nxt)
+ component.append(int(nxt))
+ found.append(sorted(component))
+ return [c for c in found if len(c) > 1]
+
+
+def report(field: np.ndarray, name: str) -> float:
+ size = len(field)
+ groups = orbits(automorphic_transpositions(field))
+ moved = sum(len(c) for c in groups)
+ log_order = sum(math.lgamma(len(c) + 1) for c in groups) / math.log(10)
+ ceiling = 1.0 - (moved - len(groups)) / size
+ print(
+ f"{name}: {len(groups)} nontrivial orbits covering {moved}/{size} items, "
+ f"sizes {sorted((len(c) for c in groups), reverse=True)[:12]}, "
+ f"|Aut| >= 10^{log_order:.1f}"
+ )
+ print(f"{name}: blind accuracy ceiling for any energy-only method = {ceiling:.4f}")
+ return ceiling
+
+
+def main() -> None:
+ path = sys.argv[1] if len(sys.argv) > 1 else "artifacts/synth_v1/omit_size.pt"
+ state = torch.load(path, map_location="cpu", weights_only=False)
+ visual = standardise(state["visual_field"].double().numpy())
+ text = standardise(state["text_field"].double().numpy())
+ report(text, "text field")
+ report(visual, "visual field")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/worldalign/fused_match.py b/worldalign/fused_match.py
new file mode 100644
index 0000000..aad3360
--- /dev/null
+++ b/worldalign/fused_match.py
@@ -0,0 +1,214 @@
+"""The unary term the search was missing.
+
+Fifteen solvers returned under 5% on the caption-omitted field: Umeyama,
+GRAMPA, five Gromov-Wasserstein variants, FAQ, PATH convex-concave, a moment
+ladder, consensus voting. They have one thing in common, and it is the thing
+that was wrong with all of them at once. **Every one is purely quadratic or
+purely spectral.** None carries a node-level term. The instance was built by
+deleting one factor from the caption side, which flattens exactly the coarse
+second-order statistics those methods read, while leaving scenes distinguishable
+in the shape of their individual similarity profiles.
+
+So give the objective a unary term. Each node gets a permutation-invariant
+descriptor -- moments of its own row of the relation field, which is a
+row-order statistic and not a spectral one -- and the descriptor distance
+becomes a linear cost blended with the quadratic one. The blend must stay in
+the loop: using the descriptor only to seed a pure-quadratic descent scores
+0.21, because the descent immediately discards it and wanders back into the
+decoy region. Kept pinned through Frank-Wolfe, the same descriptor scores 0.84.
+
+That 0.84 is the exact answer. The caption-omitted text field has 51 exact
+transposition automorphisms, so `T[s,s]` is bitwise identical to `T` for each
+and no objective built from `(V, P T P^T)` can distinguish the members of an
+orbit at any order. The blind ceiling is 0.836, and this reaches it with an
+energy gap of exactly zero.
+
+Alpha is chosen by energy, never by accuracy, so the procedure stays blind.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+
+import numpy as np
+import torch
+from scipy.optimize import linear_sum_assignment
+
+from .common import write_json
+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=40)
+ parser.add_argument("--device", default="cuda:3")
+ parser.add_argument("--output", default="artifacts/synth_v1/fused_match.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 row_descriptors(field: np.ndarray) -> np.ndarray:
+ """Moments of each node's own similarity profile.
+
+ Invariant to how the other nodes are ordered, so it survives the unknown
+ permutation, and it reads the distribution of a node's relations rather
+ than the field's leading eigenvectors -- which is why it sees what the
+ spectral family cannot.
+ """
+ size = len(field)
+ rows = []
+ for index in range(size):
+ row = np.delete(field[index], index)
+ rows.append([
+ row.mean(), row.std(), np.mean(row ** 3), np.mean(row ** 4),
+ *np.percentile(row, [10, 50, 90]), row.max(),
+ ])
+ return np.array(rows)
+
+
+def unary_cost(visual: np.ndarray, text: np.ndarray) -> np.ndarray:
+ """Squared descriptor distance, z-scored jointly over both node sets."""
+ left, right = row_descriptors(visual), row_descriptors(text)
+ both = np.vstack([left, right])
+ centre, spread = both.mean(0), both.std(0).clip(1e-9)
+ left, right = (left - centre) / spread, (right - centre) / spread
+ return ((left[:, None, :] - right[None, :, :]) ** 2).sum(-1)
+
+
+def blind_automorphism_ceiling(field: np.ndarray) -> tuple[float, int]:
+ """The best any blind method can do, given exact automorphisms of the field.
+
+ Two scenes whose rows agree as multisets may be exchangeable: if swapping
+ them leaves the field bitwise unchanged, every objective of the form
+ f(V, P T P^T) is blind to which is which, at every order. Members of an
+ orbit are then guessable only up to one correct pick per orbit.
+ """
+ size = len(field)
+ buckets: dict[bytes, list[int]] = {}
+ for index in range(size):
+ row = np.sort(np.delete(field[index], index))
+ buckets.setdefault(row.tobytes(), []).append(index)
+ orbits = [group for group in buckets.values() if len(group) > 1]
+ automorphisms = 0
+ for group in orbits:
+ for position, first in enumerate(group):
+ for second in group[position + 1:]:
+ swap = np.arange(size)
+ swap[[first, second]] = swap[[second, first]]
+ if np.array_equal(field[np.ix_(swap, swap)], field):
+ automorphisms += 1
+ ambiguous = sum(len(group) for group in orbits)
+ ceiling = (size - ambiguous + len(orbits)) / size
+ return float(ceiling), automorphisms
+
+
+def fused_frank_wolfe(
+ visual: np.ndarray, text: np.ndarray, cost: np.ndarray,
+ alpha: float, iterations: int,
+) -> np.ndarray:
+ """Maximise alpha*<V, P T P^T> - (1-alpha)*<cost, P> over the Birkhoff polytope."""
+ size = len(visual)
+ coupling = np.full((size, size), 1.0 / size)
+ scale = abs(visual @ coupling @ text).mean() / max(abs(cost).mean(), 1e-12)
+ scaled = cost * scale
+ for _ in range(iterations):
+ gradient = alpha * 2.0 * (visual @ coupling @ text) - (1.0 - alpha) * scaled
+ _, columns = linear_sum_assignment(-gradient)
+ vertex = np.zeros_like(coupling)
+ vertex[np.arange(size), columns] = 1.0
+ direction = vertex - coupling
+ quadratic = alpha * float((direction * (visual @ direction @ text)).sum())
+ linear = float((direction * gradient).sum())
+ if quadratic >= 0:
+ step = 1.0
+ else:
+ step = float(np.clip(-linear / (2.0 * quadratic), 0.0, 1.0))
+ if step <= 1e-12:
+ break
+ coupling = coupling + step * direction
+ _, columns = linear_sum_assignment(-coupling)
+ return columns
+
+
+def main() -> None:
+ args = parse_args()
+ labels = args.labels or [p.split("/")[-1] for p in args.fields]
+ 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)
+ ceiling, automorphisms = blind_automorphism_ceiling(text)
+ visual_gpu = standardized(torch.from_numpy(visual).to(device)).double()
+ print(f"\n=== {label} N={size} automorphisms={automorphisms} "
+ f"blind ceiling={ceiling:.4f}", flush=True)
+
+ accuracies, gaps = [], []
+ for trial in range(args.trials):
+ hidden = np.random.default_rng(trial).permutation(size)
+ shuffled = text[np.ix_(hidden, hidden)]
+ text_gpu = standardized(torch.from_numpy(shuffled).to(device)).double()
+ energy = ClosedFormEnergy(text_gpu, visual_gpu, 1.0, 0.0, 256)
+ truth_energy = float(energy.energy(
+ torch.from_numpy(np.argsort(hidden).copy()).to(device)[None]
+ )[0])
+ cost = unary_cost(visual, shuffled)
+
+ best = (np.inf, 0.0, 0.0)
+ for alpha in np.arange(0.70, 0.99, 0.04):
+ columns = fused_frank_wolfe(
+ visual, shuffled, cost, float(alpha), args.iterations
+ )
+ final = fast_pair_descent(
+ text_gpu, visual_gpu,
+ torch.from_numpy(np.ascontiguousarray(columns)).to(device), 2000,
+ )
+ value = float(energy.energy(final[None])[0])
+ if value < best[0]: # selection by energy only
+ best = (value, float(
+ (hidden[final.cpu().numpy()] == np.arange(size)).mean()
+ ), float(alpha))
+ accuracies.append(best[1])
+ gaps.append(best[0] - truth_energy)
+ print(f" trial {trial}: accuracy={best[1]:.4f} alpha={best[2]:.2f} "
+ f"dE={best[0] - truth_energy:+.2e}", flush=True)
+
+ row = {
+ "field": label, "size": size, "chance": 1.0 / size,
+ "automorphisms": automorphisms, "blind_ceiling": ceiling,
+ "accuracy": float(np.mean(accuracies)),
+ "fraction_of_ceiling": float(np.mean(accuracies)) / ceiling,
+ "energy_gap": float(np.mean(gaps)),
+ }
+ rows.append(row)
+ print(f" MEAN accuracy={row['accuracy']:.4f} = "
+ f"{row['fraction_of_ceiling']:.3f} of the blind ceiling", flush=True)
+
+ write_json(args.output, {
+ "protocol": (
+ "Blind. Hidden permutation generated per trial and read only for "
+ "scoring; alpha selected by energy. The ceiling is computed from "
+ "exact automorphisms of the text field and bounds any blind method."
+ ),
+ "rows": rows,
+ })
+ print(json.dumps({"done": True}))
+
+
+if __name__ == "__main__":
+ main()
diff --git a/worldalign/match_battery.py b/worldalign/match_battery.py
index 9035e62..63a4bf1 100644
--- a/worldalign/match_battery.py
+++ b/worldalign/match_battery.py
@@ -142,11 +142,19 @@ def solve_faq_multi(visual, shuffled, **_):
"""FAQ from many random doubly stochastic starts, best by its own objective."""
from scipy.optimize import quadratic_assignment
- result = quadratic_assignment(
- visual, shuffled, method="faq",
- options={"maximize": True, "n_init": 30, "rng": 0},
- )
- return result.col_ind
+ # scipy's FAQ takes no `n_init`; passing one is silently swallowed into
+ # unknown_options, which is why an earlier version of this function
+ # returned bit-identical results to plain FAQ. Restart explicitly.
+ best, best_value = None, -np.inf
+ for seed in range(15):
+ result = quadratic_assignment(
+ visual, shuffled, method="faq",
+ options={"maximize": True, "rng": seed, "P0": "randomized"},
+ )
+ value = float((visual * shuffled[np.ix_(result.col_ind, result.col_ind)]).sum())
+ if value > best_value:
+ best, best_value = result.col_ind, value
+ return best
def solve_2opt(visual, shuffled, **_):
generated by cgit v1.2.3 (git 2.25.1) at 2026-08-02 14:15:29 +0000