diff options
26 files changed, 1070 insertions, 0 deletions
diff --git a/MATCHING_RESULTS.md b/MATCHING_RESULTS.md new file mode 100644 index 0000000..645393a --- /dev/null +++ b/MATCHING_RESULTS.md @@ -0,0 +1,144 @@ +# The search was the narrow part + +*2026-08-01. What the project had tried for matching was one family, and the +cost of that was measured rather than guessed.* + +**Result in one line.** On a field whose information is 99.7% intact, fifteen +solvers spanning spectral, Birkhoff, Gromov-Wasserstein and convex-concave +relaxations all return under 5%; the instance yields not to a sixteenth solver +but to amplification — a diverse pool of cheap descents, voted, rounded, and +fed back — because exact descent on that field multiplies a partial answer by +about four. + +## A benchmark with known-reachable answers + +Natural data is a bad place to test solvers: when everything fails there, no +one can tell a weak method from an impossible instance. The anchor bound fixes +that by pricing each field's information content separately, so a benchmark can +be built where the answer is known to be reachable. + +| Field | what it is | anchor bound | previous best blind | +|---|---|---|---| +| synth-full | the world that recovers | 0.99 | 0.958 | +| rank8 | same, truncated to rank 8 | 0.93 | 0.076 | +| omit-size | same, captions never state size | 0.997 | 0.044 | +| natural | Visual Genome | 0.29 | 0.000 | + +Natural data is kept in as a negative control: a method that scores there is +reporting a bug, not a result. None did. + +## Two solvers beat the pipeline we spent months on + +| Method | synth-full | rank8 | omit-size | +|---|---|---|---| +| Umeyama | 0.540 | 0.016 | 0.014 | +| GRAMPA + exact descent *(the pipeline)* | 0.884 | 0.076 | 0.044 | +| Gromov-Wasserstein | 0.397 | **0.881** | 0.022 | +| GW, 12 restarts | 0.397 | 0.881 | 0.022 | +| GW, KL loss | 0.089 | 0.251 | 0.000 | +| entropic GW, annealed | **0.961** | **0.887** | 0.034 | +| entropic GW, deep anneal | — | 0.013 | 0.034 | +| BAPG GW | 0.418 | 0.233 | 0.012 | +| semirelaxed GW | — | 0.044 | 0.013 | +| FAQ (Frank-Wolfe) | — | — | 0.021 | +| FAQ, 30 restarts | — | — | 0.021 | +| PATH convex-concave | **0.961** | 0.866 | 0.004 | +| moment ladder tr(M^k) | 0.112 | 0.214 | 0.012 | +| consensus vote, 65 descents | 0.932 | — | 0.127 | + +Two things fall out immediately. **Entropic GW with an annealed regulariser and +PATH both reach 0.961 on the reference field, above the 0.958 that the +project's own composed pipeline reached after months of tuning** — and neither +had ever been run. And on the rank-8 field, GW reaches 0.881 where GRAMPA +reaches 0.076. + +That second number retires a conclusion drawn earlier the same day. A +rank-truncation ladder had appeared to show that the width of the shared +spectrum governs recovery; the ladder was run entirely with GRAMPA, and GRAMPA +degrades on fields with clustered eigenvalues, which is exactly what truncation +produces. **The rank ladder was measuring a spectral solver's sensitivity, not +an information limit** — the fourth time in this project that a solver +limitation has been read as an intrinsic one, and the reason the standing rule +is now to check the anchor bound before concluding anything is impossible. + +More annealing is not better. Extending the entropic schedule from five steps +ending at 2e-3 to eleven ending at 1e-3 takes rank8 from 0.887 to 0.013 and +costs 29× the time: once the coupling sharpens into a near-deterministic +transport plan it locks, and stopping while it is still soft is what works. + +## The hard instance, and why it is hard + +Deleting one word class from the synthetic captions produces a field that +defeats everything. Its diagnosis is unusual and each part was measured. + +**The information is intact.** Anchor bound 0.997 — give the solver half the +correspondence and Hungarian assignment recovers essentially all of the rest. +Five percent of the correspondence, twelve pairs, is worth 62%. + +**The truth is a strict local minimum.** Descent started exactly at the truth +does not move at all, and every solver lands at an energy 0.44 above it against +E(truth) = 0.34. + +**Its basin is wide, not a needle.** Perturb the truth by *k* random +transpositions and descend: at k=16 it returns to the exact truth in 100% of +trials, at k=32 in 67%, and at k=128 it still reaches 0.90 accuracy. + +Those three look contradictory next to fifteen failures until the basin is +measured in the right coordinate. 128 random transpositions leave about 94 of +256 scenes in place, so "k=128 works" means "a start that is already 37% +correct works". Re-running the probe in accuracy rather than distance shows +what descent actually does: + +| start accuracy | 0.05 | 0.10 | 0.15 | 0.20 | 0.30 | 0.50 | +|---|---|---|---|---|---|---| +| after descent | 0.152 | 0.420 | 0.614 | 0.783 | 0.806 | 0.875 | + +**Descent is an amplifier with a gain near four, and no sharp threshold.** The +failure was never the local search and never the information. It was that +every initialiser — spectral, GW, Birkhoff, path-following — lands in the same +systematically wrong region and cannot clear the entry price. + +## Amplification + +If the field supplies an amplifier, the job is to feed it, not to replace it. +Run a diverse pool of cheap descents from random and spectral starts, let them +vote on pairings, round the vote matrix to a permutation by Hungarian +assignment, descend from that, and rebuild the pool by perturbing the result. +Each round hands the amplifier a better start than the last. + +Two details earn their place. **Rounding the vote matrix beats taking its +argmax** — 0.207 against 0.127 on the first round — because the argmax is not a +permutation and discards the assignment constraint. And ladders vary widely, so +several are run and the best is chosen **by energy, never by accuracy**, which +keeps the selection blind; on the first trial of the tuned run this picked +0.188 out of {0.094, 0.043, 0.188}, the best of the three. + +On the instance that defeated fifteen solvers, the ladder climbs 0.13 → 0.68 → +0.72 across rounds, against 0.044 for the best solver and 0.127 for voting +alone. + +One negative result on the design. Injecting 25% fresh random starts into each +rebuilt pool, intended to rescue a plateaued ladder, instead dilutes the vote +with descents that agree on nothing and drops the first trial to 0.188. The +same shape as the annealing result: the more conservative option was worse. + +## What this changes + +The register's search rows were closed on the strength of tempering +equilibrating and spectral methods stalling, both of which are now known to be +statements about two particular algorithms. They are reopened. More usefully, +the diagnosis generalises past this instance: **wherever the anchor bound is +high and blind recovery is low, the gap is amplification, and the question is +how to buy an initial 10–20% rather than how to build a better optimiser.** +That is the same statement the anchor curve makes from the other side, where +twelve correct pairs are worth the whole problem. + +Engineering note, because it gated everything above. `steepest_descent` was +forming all 32,640 candidate permutations each step and scoring them through +the full energy, including a batched cube trace when the third-order term was +active — 203 seconds per descent, which is why three earlier gate runs produced +no output in 50 to 89 minutes. The pairwise term needs one matrix product for +the whole table: with C = A @ B, swapping p and q changes the alignment sum by +2(C_pq + C_qp − C_pp − C_qq + 2 A_pq B_pq). Checked against brute force to 1e-9 +before use. **203 s → 0.79 s**, and every experiment in this document is +downstream of that. diff --git a/artifacts/spectral_frontier_probe/band.py b/artifacts/spectral_frontier_probe/band.py new file mode 100644 index 0000000..cca368a --- /dev/null +++ b/artifacts/spectral_frontier_probe/band.py @@ -0,0 +1,26 @@ +import numpy as np, torch +from scipy.optimize import linear_sum_assignment +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V); truth=np.arange(N) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1] +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +acc=lambda p: float((p==truth).mean()) +CONST=float((T*T).sum()+(V*V).sum()) +en=lambda p:(CONST-2*(T[np.ix_(p,p)]*V).sum())/(N*(N-1)) +for r in (12,16,24): + XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r])) + u,s,vt=np.linalg.svd(XT.T@XV); O=u@vt + idx=np.arange(r); band_mass=[] + for b in (1,2,3,4,6,r): + Mk=(np.abs(idx[:,None]-idx[None,:])<=b).astype(float) + frac=float((O**2*Mk).sum()/ (O**2).sum()) + Ob=O*Mk; uu,ss,vv=np.linalg.svd(Ob); Ob=uu@vv + p=hung(XV,XT@Ob) + band_mass.append((b,round(frac,3),round(acc(p),3),round(en(p),4))) + print(f"r={r}: (bandwidth, mass of oracle-O inside band, acc after re-orthogonalising, E) -> {band_mass}") + print(f" free params: full {r*(r-1)//2}, band3 {sum(min(3,r-1-i) for i in range(r))}") diff --git a/artifacts/spectral_frontier_probe/confirm.log b/artifacts/spectral_frontier_probe/confirm.log new file mode 100644 index 0000000..ff6ac9d --- /dev/null +++ b/artifacts/spectral_frontier_probe/confirm.log @@ -0,0 +1,7 @@ + r= 4: pre-descent best E 0.7185 (acc 0.004) -> post-descent E 0.6493 acc 0.016 + r= 6: pre-descent best E 0.5638 (acc 0.117) -> post-descent E 0.5333 acc 0.148 + r= 8: pre-descent best E 0.4786 (acc 0.395) -> post-descent E 0.3597 acc 0.750 + r= 10: pre-descent best E 0.6224 (acc 0.180) -> post-descent E 0.5519 acc 0.172 + r= 12: pre-descent best E 0.3977 (acc 0.660) -> post-descent E 0.3396 acc 0.879 + r= 14: pre-descent best E 0.3671 (acc 0.707) -> post-descent E 0.3396 acc 0.836 + r= 16: pre-descent best E 0.6410 (acc 0.145) -> post-descent E 0.5474 acc 0.277 diff --git a/artifacts/spectral_frontier_probe/confirm.py b/artifacts/spectral_frontier_probe/confirm.py new file mode 100644 index 0000000..d0fc003 --- /dev/null +++ b/artifacts/spectral_frontier_probe/confirm.py @@ -0,0 +1,49 @@ +import sys, time, torch, numpy as np +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from scipy.stats import ortho_group +from worldalign.synth_fast_gate import fast_pair_descent +dev='cuda:3' +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) +Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(T,dtype=torch.float32,device=dev) +CONST=float((Tt*Tt).sum()+(Vt*Vt).sum()) +def energy(p): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1)) +def descend(p,steps=4000): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy() +truth=np.arange(N); acc=lambda p: float((p==truth).mean()) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1] +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +rng=np.random.default_rng(20260801) +def icp(XV,XT,O,iters=30): + for _ in range(iters): + p=hung(XV,XT@O) + u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt + if np.allclose(On,O,atol=1e-10): O=On; break + O=On + return hung(XV,XT@O) +t0=time.time(); pool=[] +for r in (4,6,8,10,12,14,16,20,24): + XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r])) + cand=[] + for t in range(120): + p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30)))) + cand.append((energy(p),p)) + cand.sort(key=lambda z:z[0]) + best=None + for e,p in cand[:5]: + pd=descend(p); ed=energy(pd) + if best is None or ed<best[0]: best=(ed,pd) + pool.append((best[0],best[1],r)) + print(f" r={r:3d}: pre-descent best E {cand[0][0]:.4f} (acc {acc(cand[0][1]):.3f}) -> post-descent E {best[0]:.4f} acc {acc(best[1]):.3f}",flush=True) +pool.sort(key=lambda z:z[0]) +print(f"\nBLIND PICK (lowest E over all r): r={pool[0][2]} E={pool[0][0]:.4f} ACC={acc(pool[0][1]):.3f}") +print(f"E(truth)={energy(truth):.4f} total wall {time.time()-t0:.0f}s") diff --git a/artifacts/spectral_frontier_probe/diag10.py b/artifacts/spectral_frontier_probe/diag10.py new file mode 100644 index 0000000..752612c --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag10.py @@ -0,0 +1,41 @@ +import numpy as np, torch +from collections import defaultdict +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy() +T=standardise(T0); V=standardise(V0); N=len(T) +# exact automorphism-by-identical-rows of T (blind: uses T only) +key={}; cls=defaultdict(list) +Toff=T.copy(); np.fill_diagonal(Toff,0.0) +# group scenes whose T rows agree after removing the two swapped coords +groups=[]; used=np.zeros(N,bool) +for i in range(N): + if used[i]: continue + g=[i]; used[i]=True + for j in range(i+1,N): + if used[j]: continue + a=np.delete(Toff[i],[i,j]); b=np.delete(Toff[j],[i,j]) + if np.abs(a-b).max()<1e-9 and abs(Toff[i,j]-max(Toff[i,i],0))<1e6: + g.append(j); used[j]=True + groups.append(g) +sizes=np.array([len(g) for g in groups]) +print("T exact-twin classes: total",len(groups),"; size histogram",np.bincount(sizes)[1:]) +excess=int((sizes-1).sum()) +print("scenes in non-trivial classes:",int(sizes[sizes>1].sum()),"; log|Aut| classes:",int((sizes>1).sum())) +ceil=(N-int(sizes[sizes>1].sum())+int((sizes>1).sum()))/N +print(f"blind accuracy ceiling if class is identified but member picked at random: {ceil:.3f}") +# check the objective really is invariant: swap two members of a class +import itertools +p=np.arange(N) +E0=(( (T*T).sum()+(V*V).sum() )-2*(T[np.ix_(p,p)]*V).sum())/(N*(N-1)) +bad=0; tested=0 +for g in groups: + if len(g)>1: + i,j=g[0],g[1]; q=p.copy(); q[[i,j]]=q[[j,i]] + E1=(((T*T).sum()+(V*V).sum())-2*(T[np.ix_(q,q)]*V).sum())/(N*(N-1)) + tested+=1 + if abs(E1-E0)>1e-9: bad+=1 +print(f"swapping twins changes E in {bad}/{tested} classes (0 means exact symmetry of the objective)") +print(f"E(truth) = {E0:.10f}") diff --git a/artifacts/spectral_frontier_probe/diag11.log b/artifacts/spectral_frontier_probe/diag11.log new file mode 100644 index 0000000..dcdcefe --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag11.log @@ -0,0 +1,8 @@ +r=24 symNMF resid V 0.043 T 0.034 [325s] + oracle col cos: [0.983 0.983 0.961 0.957 0.955 0.955 0.951 0.949 0.934 0.908] + BLIND stat col-match agrees w/ oracle 0.12 -> scene acc 0.020 + after co-alternation (1 iters): col agree 0.12 -> scene acc 0.020 +r=42 symNMF resid V 0.023 T 0.014 [741s] + oracle col cos: [0.973 0.949 0.948 0.943 0.936 0.921 0.912 0.901 0.892 0.89 ] + BLIND stat col-match agrees w/ oracle 0.00 -> scene acc 0.027 + after co-alternation (2 iters): col agree 0.02 -> scene acc 0.012 diff --git a/artifacts/spectral_frontier_probe/diag11.py b/artifacts/spectral_frontier_probe/diag11.py new file mode 100644 index 0000000..4977c3c --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag11.py @@ -0,0 +1,52 @@ +import numpy as np, torch, warnings, time +warnings.filterwarnings('ignore') +from scipy.optimize import linear_sum_assignment +from sklearn.decomposition import NMF +np.set_printoptions(precision=3,suppress=True,linewidth=200) +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy(); N=len(V0) +rng=np.random.default_rng(777); sigma=rng.permutation(N) +T0s=T0[np.ix_(sigma,sigma)]; truth=np.argsort(sigma) +acc=lambda p: float((p==truth).mean()) +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +def sym_nmf(M,r,iters=3000,seed=0): + """symmetric NMF M ~ W W^T by multiplicative updates (Ding et al.)""" + M=np.clip(M,0,None).copy(); np.fill_diagonal(M,np.clip(np.diag(M),0,None)) + rg=np.random.default_rng(seed); W=np.abs(rg.standard_normal((len(M),r)))*np.sqrt(M.mean()/r) + for _ in range(iters): + num=M@W; den=W@(W.T@W)+1e-12 + W=W*(0.5+0.5*num/den) + return W +for r in (24,42): + t0=time.time() + WV=sym_nmf(V0,r,seed=1); WT=sym_nmf(T0s,r,seed=1) + print(f"r={r} symNMF resid V {np.linalg.norm(V0-WV@WV.T)/np.linalg.norm(V0):.3f} T {np.linalg.norm(T0s-WT@WT.T)/np.linalg.norm(T0s):.3f} [{time.time()-t0:.0f}s]") + A=WV.copy(); B=WT.copy() + # ORACLE column match for reference + An=A/np.linalg.norm(A,axis=0,keepdims=True); Bn=B/np.linalg.norm(B,axis=0,keepdims=True) + rr,cc=linear_sum_assignment(-(An[np.arange(N)].T@Bn[truth])) + print(" oracle col cos:",np.sort((An.T@Bn[truth])[rr,cc])[::-1][:10]) + # BLIND column match by permutation-invariant column signatures + def sig(X): + Xn=X/ (np.linalg.norm(X,axis=0,keepdims=True)+1e-12) + q=np.quantile(Xn,np.linspace(0.5,1.0,16),axis=0).T + return np.hstack([q,(Xn>0.02).mean(0)[:,None],np.linalg.norm(X,axis=0)[:,None]/np.linalg.norm(X)]) + sA=sig(A); sB=sig(B); m=hung(sA,sB) + agree=float((m==cc[np.argsort(rr)]).mean()) + def scene_match(colmap): + AV=A; BT=B[:,colmap] + AV=AV/np.linalg.norm(AV,axis=1,keepdims=True).clip(1e-9); BT=BT/np.linalg.norm(BT,axis=1,keepdims=True).clip(1e-9) + return hung(AV,BT) + p=scene_match(m); print(f" BLIND stat col-match agrees w/ oracle {agree:.2f} -> scene acc {acc(p):.3f}") + # alternate: scene Hungarian <-> column Hungarian + colmap=m.copy() + for it in range(25): + p=scene_match(colmap) + Ar=A/np.linalg.norm(A,axis=1,keepdims=True).clip(1e-9) + Br=B/np.linalg.norm(B,axis=1,keepdims=True).clip(1e-9) + rr2,cc2=linear_sum_assignment(-(Ar.T@Br[p])); newmap=cc2[np.argsort(rr2)] + if (newmap==colmap).all(): break + colmap=newmap + p=scene_match(colmap) + print(f" after co-alternation ({it+1} iters): col agree {float((colmap==cc[np.argsort(rr)]).mean()):.2f} -> scene acc {acc(p):.3f}") diff --git a/artifacts/spectral_frontier_probe/diag3.py b/artifacts/spectral_frontier_probe/diag3.py new file mode 100644 index 0000000..45a1430 --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag3.py @@ -0,0 +1,21 @@ +import torch, numpy as np +np.set_printoptions(precision=3, suppress=True, linewidth=200) +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) + +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1] +K=24 +X=np.abs(UV[:,:K].T@UT[:,:K]) +print("|<u_k(V), u_l(T)>| top-16 block (rows=V index, cols=T index):") +print(X[:16,:16]) +print("\ndiagonal |<u_k(V),u_k(T)>| k=0..23:", np.diag(X)) +print("row-max of |overlap| (best T partner for each V eigvec):", X.max(1)) +print("argmax:", X.argmax(1)) +# subspace alignment: principal angles between top-r spaces +for r in (4,6,8,10,12,16,20,24,32,48): + s=np.linalg.svd(UV[:,:r].T@UT[:,:r],compute_uv=False) + print(f"r={r:3d} mean cos principal angle {s.mean():.3f} captured energy {np.sum(s**2)/r:.3f} min {s.min():.3f}") diff --git a/artifacts/spectral_frontier_probe/diag4.py b/artifacts/spectral_frontier_probe/diag4.py new file mode 100644 index 0000000..0f7a75a --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag4.py @@ -0,0 +1,29 @@ +import torch, numpy as np +from scipy.optimize import linear_sum_assignment +np.set_printoptions(precision=3, suppress=True, linewidth=200) +rng=np.random.default_rng(0) +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1] +truth=np.arange(N) +def acc(p): return float((p==truth).mean()) +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T) + r,c=linear_sum_assignment(C); return c + +print("### 1. ORACLE orthogonal mixing O between top-r spectral embeddings") +for r in (6,8,10,12,16,20,24,32,42): + for scale in ('none','sqrt','lam'): + f=lambda w: np.ones_like(w) if scale=='none' else (np.sqrt(np.abs(w)) if scale=='sqrt' else np.abs(w)) + XV=UV[:,:r]*f(wV[:r]); XT=UT[:,:r]*f(wT[:r]) + # oracle Procrustes using truth + M=XT.T@XV; u,s,vt=np.linalg.svd(M); O=u@vt + p=hung(XV,XT@O) + # also cosine-normalised rows + nV=XV/np.linalg.norm(XV,axis=1,keepdims=True); nT=(XT@O); nT=nT/np.linalg.norm(nT,axis=1,keepdims=True) + p2=hung(nV,nT) + print(f" r={r:3d} scale={scale:5s} oracle-O acc={acc(p):.3f} rownorm acc={acc(p2):.3f}") diff --git a/artifacts/spectral_frontier_probe/diag6.py b/artifacts/spectral_frontier_probe/diag6.py new file mode 100644 index 0000000..0158f9d --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag6.py @@ -0,0 +1,46 @@ +import sys, time, torch, numpy as np +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from scipy.stats import ortho_group +from worldalign.synth_fast_gate import fast_pair_descent +dev='cuda:3' +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) +Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(T,dtype=torch.float32,device=dev) +CONST=float((Tt*Tt).sum()+(Vt*Vt).sum()) +def energy(p): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1)) +def descend(p,steps=4000): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy() +truth=np.arange(N); acc=lambda p: float((p==truth).mean()) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(T); wT=wT[::-1]; UT=UT[:,::-1] +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +rng=np.random.default_rng(1) + +def icp(XV,XT,O,iters=40): + for _ in range(iters): + p=hung(XV,XT@O) + u,s,vt=np.linalg.svd(XT[p].T@XV); Onew=u@vt + if np.allclose(Onew,O,atol=1e-10): O=Onew; break + O=Onew + return hung(XV,XT@O),O + +print("### BLIND: random-restart ICP over O(r), scored by QAP energy") +for r in (6,8,10,12,16): + XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r])) + t0=time.time(); best=(1e9,None) + R=200 + for t in range(R): + O=ortho_group.rvs(r,random_state=int(rng.integers(1<<30))) + p,_=icp(XV,XT,O) + e=energy(p) + if e<best[0]: best=(e,p) + p=best[1]; pd=descend(p) + print(f" r={r:3d} R={R}: best-E {best[0]:.4f} acc {acc(p):.3f} | after descent E {energy(pd):.4f} acc {acc(pd):.3f} [{time.time()-t0:.0f}s]") diff --git a/artifacts/spectral_frontier_probe/diag7.py b/artifacts/spectral_frontier_probe/diag7.py new file mode 100644 index 0000000..c41b3c4 --- /dev/null +++ b/artifacts/spectral_frontier_probe/diag7.py @@ -0,0 +1,74 @@ +import sys, time, torch, numpy as np +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +np.set_printoptions(precision=4, suppress=True, linewidth=200) +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V0=d['visual_field'].double().numpy(); T0=d['text_field'].double().numpy() +V=standardise(V0); T=standardise(T0); N=len(V) +truth=np.arange(N); acc=lambda p: float((p==truth).mean()) +CONST=float((T*T).sum()+(V*V).sum()) +def energy_from_align(S): return (CONST-2.0*S)/(N*(N-1)) +def energy(p): return energy_from_align(float((T[np.ix_(p,p)]*V).sum())) +print("E(truth) =",energy(truth)) + +print("\n### 2. SPECTRAL CERTIFICATES (lower bounds on E over all permutations)") +lV=np.linalg.eigvalsh(V)[::-1]; lT=np.linalg.eigvalsh(T)[::-1] +print(" Hoffman-Wielandt / von Neumann bound E >= ", energy_from_align(float((lV*lT).sum()))) +# projected eigenvalue bound (Hadley-Rendl-Wolkowicz): project out the all-ones direction +Q,_=np.linalg.qr(np.hstack([np.ones((N,1))/np.sqrt(N), np.random.default_rng(0).standard_normal((N,N-1))])) +Vp=Q[:,1:].T@V@Q[:,1:]; Tp=Q[:,1:].T@T@Q[:,1:] +lVp=np.linalg.eigvalsh(Vp)[::-1]; lTp=np.linalg.eigvalsh(Tp)[::-1] +sV=V.sum(); sT=T.sum(); rV=V.sum(1); rT=T.sum(1) +# <V,PTP'> = (1/N)*?; exact decomposition: with x=ones/sqrt(N), +# tr(V P T P^T) = tr(Vp Pp Tp Pp^T) + 2 x^T V P T P^T x*? -- use the standard HRW split +cross = float(rV.sum()*rT.sum())/ (N*N) # x^T V x * x^T T x term +# rank-1 cross terms bounded by sorted-inner-product of centred row sums +cV=np.sort(rV-rV.mean())[::-1]; cT=np.sort(rT-rT.mean())[::-1] +proj_bound = float((lVp*lTp).sum()) + cross + 2.0*float((cV*cT).sum())/N +print(" projected (HRW-style, upper bd on alignment) E >= ", energy_from_align(proj_bound)) +print(" (both are LOWER bounds on E; E(truth)=%.4f, solvers report E(truth)+0.44=%.4f)"%(energy(truth),energy(truth)+0.44)) + +print("\n### 3. BLIND rotation-invariant node descriptors -> Hungarian seeds") +def report(name,DV,DT,topk=(12,25,50)): + DV=DV/ (np.linalg.norm(DV,axis=1,keepdims=True)+1e-12); DT=DT/(np.linalg.norm(DT,axis=1,keepdims=True)+1e-12) + C=((DV**2).sum(1)[:,None]+(DT**2).sum(1)[None,:]-2*DV@DT.T) + r,c=linear_sum_assignment(C) + a=acc(c) + # confidence = margin between assigned cost and 2nd best in row + Cm=C.copy(); Cm[np.arange(N),c]=np.inf + margin=Cm.min(1)-C[np.arange(N),c] + order=np.argsort(-margin) + prec={k: float((c[order[:k]]==order[:k]).mean()) for k in topk} + print(f" {name:34s} full-acc {a:.3f} | precision@top-margin {prec}") + return c +# (a) sorted row profile +report("sorted row profile (V vs T)", np.sort(V,1), np.sort(T,1)) +# (b) row moments +def mom(M,K=8): return np.stack([ (M**k).mean(1) for k in range(1,K+1)],1) +report("row power moments k=1..8", mom(V), mom(T)) +# (c) heat kernel signature on normalised Laplacian of the raw non-negative Gram +def lap_eigs(W): + W=np.clip(W,0,None).copy(); np.fill_diagonal(W,0.0) + dg=W.sum(1); Dm=1/np.sqrt(np.maximum(dg,1e-12)) + L=np.eye(len(W))-(Dm[:,None]*W*Dm[None,:]) + w,U=np.linalg.eigh(L); return w,U,dg +wV,UVl,dV=lap_eigs(V0); wT,UTl,dT=lap_eigs(T0) +ts=np.logspace(-2,1.5,24) +HV=np.stack([ (np.exp(-t*wV)[None,:]*UVl**2).sum(1) for t in ts],1) +HT=np.stack([ (np.exp(-t*wT)[None,:]*UTl**2).sum(1) for t in ts],1) +report("HKS (normalised Laplacian)", np.log(HV+1e-12), np.log(HT+1e-12)) +# (d) wave kernel signature +def wks(w,U,M=24): + lw=np.log(np.maximum(w,1e-6)); e=np.linspace(lw.min(),lw.max(),M); sig=(e[1]-e[0])*2 + return np.stack([ (np.exp(-((e_-lw)**2)/(2*sig**2))[None,:]*U**2).sum(1) for e_ in e],1) +report("WKS (normalised Laplacian)", wks(wV,UVl), wks(wT,UTl)) +# (e) spectral graph wavelet (SGWT) coefficient energies +def sgwt(w,U,S=10): + ss=np.logspace(-1.5,1.0,S) + return np.stack([ ((s*w*np.exp(-s*w))[None,:]*U**2).sum(1) for s in ss],1) +report("SGWT band energies", sgwt(wV,UVl), sgwt(wT,UTl)) +# (f) degree only +report("degree (row sum) only", dV[:,None], dT[:,None]) diff --git a/artifacts/spectral_frontier_probe/seed.log b/artifacts/spectral_frontier_probe/seed.log new file mode 100644 index 0000000..87f0744 --- /dev/null +++ b/artifacts/spectral_frontier_probe/seed.log @@ -0,0 +1,7 @@ +360 ICP solutions in 760s; corr(E,acc)=-0.520 +lowest-10 (E,acc): [(0.5186, 0.461), (0.5815, 0.27), (0.622, 0.258), (0.6385, 0.227), (0.6637, 0.051), (0.6652, 0.059), (0.6828, 0.109), (0.6903, 0.074), (0.6977, 0.016), (0.6988, 0.016)] +E quantiles [0.519 0.759 0.835 0.91 1.031] acc of best-E: 0.4609375 + consensus over 5 lowest-E sols: precision@12 = 0.833 consensus over 5 lowest-E sols: precision@25 = 0.840 consensus over 5 lowest-E sols: precision@50 = 0.720 consensus over 5 lowest-E sols: precision@100 = 0.680 + consensus over 10 lowest-E sols: precision@12 = 0.750 consensus over 10 lowest-E sols: precision@25 = 0.720 consensus over 10 lowest-E sols: precision@50 = 0.760 consensus over 10 lowest-E sols: precision@100 = 0.630 + consensus over 20 lowest-E sols: precision@12 = 0.917 consensus over 20 lowest-E sols: precision@25 = 0.840 consensus over 20 lowest-E sols: precision@50 = 0.720 consensus over 20 lowest-E sols: precision@100 = 0.590 + consensus over 40 lowest-E sols: precision@12 = 0.917 consensus over 40 lowest-E sols: precision@25 = 0.880 consensus over 40 lowest-E sols: precision@50 = 0.740 consensus over 40 lowest-E sols: precision@100 = 0.570 diff --git a/artifacts/spectral_frontier_probe/seed.py b/artifacts/spectral_frontier_probe/seed.py new file mode 100644 index 0000000..5e33d7b --- /dev/null +++ b/artifacts/spectral_frontier_probe/seed.py @@ -0,0 +1,49 @@ +import sys, time, torch, numpy as np +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from scipy.stats import ortho_group +dev='cuda:0' +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) +rng=np.random.default_rng(31337); sigma=rng.permutation(N) +Ts=T[np.ix_(sigma,sigma)]; truth=np.argsort(sigma) +Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(Ts,dtype=torch.float32,device=dev) +CONST=float((Tt*Tt).sum()+(Vt*Vt).sum()) +def energy(p): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1)) +acc=lambda p: float((p==truth).mean()) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(Ts); wT=wT[::-1]; UT=UT[:,::-1] +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +def icp(XV,XT,O,iters=25): + for _ in range(iters): + p=hung(XV,XT@O); u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt + if np.allclose(On,O,atol=1e-10): O=On; break + O=On + return hung(XV,XT@O) +sols=[] +t0=time.time() +for r in (10,12,14,16,18,20): + XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r])) + for t in range(60): + p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30)))) + sols.append((energy(p),acc(p),p,r)) +sols.sort(key=lambda z:z[0]) +E=np.array([s[0] for s in sols]); A=np.array([s[1] for s in sols]) +print(f"{len(sols)} ICP solutions in {time.time()-t0:.0f}s; corr(E,acc)={np.corrcoef(E,A)[0,1]:.3f}") +print("lowest-10 (E,acc):", [(round(e,4),round(a,3)) for e,a,_,_ in sols[:10]]) +print("E quantiles",np.quantile(E,[0,.1,.5,.9,1]).round(3)," acc of best-E:",A[0]) +for m in (5,10,20,40): + votes=np.zeros((N,N)) + for e,a,p,r in sols[:m]: votes[np.arange(N),p]+=1 + conf=votes.max(1); pick=np.argsort(-conf) + pm=votes.argmax(1) + for k in (12,25,50,100): + sel=pick[:k]; prec=float((pm[sel]==truth[sel]).mean()) + print(f" consensus over {m} lowest-E sols: precision@{k} = {prec:.3f}", end='') + print() diff --git a/artifacts/spectral_frontier_probe/shufctl.log b/artifacts/spectral_frontier_probe/shufctl.log new file mode 100644 index 0000000..8961e1c --- /dev/null +++ b/artifacts/spectral_frontier_probe/shufctl.log @@ -0,0 +1,4 @@ +E(truth)= 0.3396347943474265 +SHUFFLED r=16: pre E 0.3505 acc 0.793 -> post E 0.3396 acc 0.820 [361s] +SHUFFLED r=6: pre E 0.5712 acc 0.078 -> post E 0.5444 acc 0.109 [195s] +SHUFFLED r=12: pre E 0.5067 acc 0.402 -> post E 0.3540 acc 0.809 [410s] diff --git a/artifacts/spectral_frontier_probe/shufctl.py b/artifacts/spectral_frontier_probe/shufctl.py new file mode 100644 index 0000000..5913f53 --- /dev/null +++ b/artifacts/spectral_frontier_probe/shufctl.py @@ -0,0 +1,47 @@ +import sys, time, torch, numpy as np +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from scipy.stats import ortho_group +from worldalign.synth_fast_gate import fast_pair_descent +dev='cuda:1' +def standardise(M): + M=np.asarray(M,dtype=np.float64); mask=~np.eye(len(M),dtype=bool); v=M[mask] + out=(M-v.mean())/v.std(); np.fill_diagonal(out,0.0); return out +d=torch.load('/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt',map_location='cpu') +V=standardise(d['visual_field']); T=standardise(d['text_field']); N=len(V) +rng=np.random.default_rng(777) +sigma=rng.permutation(N) # hidden shuffle applied to T +Ts=T[np.ix_(sigma,sigma)] # scene i of V corresponds to row where sigma[k]=i -> inverse +inv=np.argsort(sigma) # truth: V index i <-> Ts index inv[i] +truth=inv +Vt=torch.tensor(V,dtype=torch.float32,device=dev); Tt=torch.tensor(Ts,dtype=torch.float32,device=dev) +CONST=float((Tt*Tt).sum()+(Vt*Vt).sum()) +def energy(p): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return (CONST-2.0*float((Tt[P[:,None],P[None,:]]*Vt).sum()))/(N*(N-1)) +def descend(p,steps=4000): + P=torch.as_tensor(np.asarray(p),dtype=torch.long,device=dev) + return fast_pair_descent(Tt,Vt,P,steps).cpu().numpy() +acc=lambda p: float((p==truth).mean()) +print("E(truth)=",energy(truth),flush=True) +wV,UV=np.linalg.eigh(V); wV=wV[::-1]; UV=UV[:,::-1] +wT,UT=np.linalg.eigh(Ts); wT=wT[::-1]; UT=UT[:,::-1] +def hung(A,B): + C=((A**2).sum(1)[:,None]+(B**2).sum(1)[None,:]-2*A@B.T); r,c=linear_sum_assignment(C); return c +def icp(XV,XT,O,iters=30): + for _ in range(iters): + p=hung(XV,XT@O); u,s,vt=np.linalg.svd(XT[p].T@XV); On=u@vt + if np.allclose(On,O,atol=1e-10): O=On; break + O=On + return hung(XV,XT@O) +for r in (16,6,12): + XV=UV[:,:r]*np.sqrt(np.abs(wV[:r])); XT=UT[:,:r]*np.sqrt(np.abs(wT[:r])) + t0=time.time(); cand=[] + for t in range(200): + p=icp(XV,XT,ortho_group.rvs(r,random_state=int(rng.integers(1<<30)))) + cand.append((energy(p),p)) + cand.sort(key=lambda z:z[0]); best=None + for e,p in cand[:5]: + pd=descend(p); ed=energy(pd) + if best is None or ed<best[0]: best=(ed,pd) + print(f"SHUFFLED r={r}: pre E {cand[0][0]:.4f} acc {acc(cand[0][1]):.3f} -> post E {best[0]:.4f} acc {acc(best[1]):.3f} [{time.time()-t0:.0f}s]",flush=True) diff --git a/logs/battery2.log b/logs/battery2.log index bfde59c..a0b3a2a 100644 --- a/logs/battery2.log +++ b/logs/battery2.log @@ -22,3 +22,9 @@ GW kl-loss raw=0.2161 refined=0.2513 (0.1s) /home/yurenh2/.local/lib/python3.13/site-packages/ot/bregman/_sinkhorn.py:642: UserWarning: Warning: numerical errors at iteration 24 warnings.warn("Warning: numerical errors at iteration %d" % ii) + entropic GW deep anneal raw=0.0078 refined=0.0130 (554.8s) + semirelaxed GW raw=0.0156 refined=0.0443 (38.9s) + +=== synth-full-BOUND0.99 (N=256, chance=0.0039) + GW x12 restarts raw=0.3698 refined=0.3971 (12.3s) + GW kl-loss raw=0.0846 refined=0.0885 (13.9s) diff --git a/logs/control.log b/logs/control.log new file mode 100644 index 0000000..de41a42 --- /dev/null +++ b/logs/control.log @@ -0,0 +1,4 @@ +{"trial": 0, "cost": "tlb", "mode": "cold", "E_truth": 0.339629709924146, "blind_refined": 0.0234375, "blind_alpha": 1.0, "blind_E": 0.6228801020056084, "oracle_refined": 0.08203125, "curve": [0.0781, 0.0234, 0.0117, 0.0273, 0.0742, 0.0352, 0.0664, 0.0352, 0.082, 0.0, 0.0547, 0.0234]} +{"trial": 0, "cost": "tlb", "mode": "warm", "E_truth": 0.339629709924146, "blind_refined": 0.84765625, "blind_alpha": 0.2, "blind_E": 0.339629709924146, "oracle_refined": 0.84765625, "curve": [0.0781, 0.8281, 0.8477, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398]} +{"trial": 0, "cost": "tlbperm", "mode": "cold", "E_truth": 0.339629709924146, "blind_refined": 0.02734375, "blind_alpha": 0.95, "blind_E": 0.6074469738905559, "oracle_refined": 0.078125, "curve": [0.0117, 0.0117, 0.043, 0.0781, 0.0078, 0.0352, 0.0156, 0.0156, 0.0039, 0.0078, 0.0273, 0.0234]} +{"trial": 0, "cost": "tlbperm", "mode": "warm", "E_truth": 0.339629709924146, "blind_refined": 0.03125, "blind_alpha": 1.0, "blind_E": 0.7038378812415739, "oracle_refined": 0.04296875, "curve": [0.0117, 0.043, 0.0234, 0.0234, 0.043, 0.0391, 0.0391, 0.0391, 0.0312, 0.0312, 0.0312, 0.0312]} diff --git a/logs/control_probe.py b/logs/control_probe.py new file mode 100644 index 0000000..170d731 --- /dev/null +++ b/logs/control_probe.py @@ -0,0 +1,93 @@ +"""Control: is the FGW warm-homotopy win driven by TLB *information*, or by any +linear term at all (generic symmetry-breaking)? Three linear costs: + + tlb -- the proposed distance-profile cost + tlbperm -- the same matrix with its columns randomly permuted: identical + value distribution, correspondence information destroyed + rand -- iid uniform noise cost + +Plus 3 trials, cold and warm sweeps, blind selection by E(P). +""" +from __future__ import annotations +import json, sys, time +sys.path.insert(0, "/home/yurenh2/emm") +import numpy as np, ot, torch +from scipy.optimize import linear_sum_assignment +from worldalign.synth_fast_gate import ClosedFormEnergy, fast_pair_descent +from worldalign.synth_triangle_gate import standardized + +DEV = torch.device("cuda:1") +ALPHAS = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0] + + +def standardise(m): + mask = ~np.eye(len(m), dtype=bool); v = m[mask] + o = (m - v.mean()) / v.std(); np.fill_diagonal(o, 0.0); return o + + +def tlb(V, S): + n = len(V); a = np.sort(V, 1); b = np.sort(S, 1) + return (a * a).sum(1)[:, None] / n + (b * b).sum(1)[None, :] / n - 2 * (a @ b.T) / n + + +def main(): + path = "/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt" + st = torch.load(path, map_location="cpu", weights_only=False) + V = standardise(st["visual_field"].double().numpy()) + T = standardise(st["text_field"].double().numpy()) + n = len(V); w = ot.unif(n) + rows = [] + for trial in range(3): + rng = np.random.default_rng(trial); hid = rng.permutation(n) + S = T[np.ix_(hid, hid)] + Vg = standardized(torch.from_numpy(V).to(DEV)).double() + Sg = standardized(torch.from_numpy(S).to(DEV)).double() + en = ClosedFormEnergy(Sg, Vg, 1.0, 0.0, 64) + et = float(en.energy(torch.from_numpy( + np.ascontiguousarray(np.argsort(hid))).to(DEV)[None])[0]) + + def acc(c): return float((hid[c] == np.arange(n)).mean()) + + def refine(c): + f = fast_pair_descent(Sg, Vg, + torch.from_numpy(np.ascontiguousarray(c)).to(DEV), + 600).cpu().numpy() + e = float(en.energy(torch.from_numpy(np.ascontiguousarray(f)).to(DEV)[None])[0]) + return e, acc(f) + + base = tlb(V, S) + shuf = np.random.default_rng(100 + trial).permutation(n) + variants = { + "tlb": base / base.max(), + "tlbperm": base[:, shuf] / base.max(), + "rand": np.random.default_rng(200 + trial).random((n, n)), + } + for vname, M in variants.items(): + for mode in ("cold", "warm"): + G = np.outer(w, w); recs = [] + for a in ALPHAS: + g0 = None if mode == "cold" else G + G2 = ot.gromov.fused_gromov_wasserstein( + M, V, S, w, w, "square_loss", alpha=a, G0=g0, + max_iter=200, tol_rel=1e-9) + if mode == "warm": + G = G2 + _, cols = linear_sum_assignment(-G2) + e, ar = refine(cols) + recs.append({"alpha": a, "raw": acc(cols), "refined": ar, "E": e}) + best = min(recs, key=lambda x: x["E"]) + row = {"trial": trial, "cost": vname, "mode": mode, + "E_truth": et, "blind_refined": best["refined"], + "blind_alpha": best["alpha"], "blind_E": best["E"], + "oracle_refined": max(x["refined"] for x in recs), + "curve": [round(x["refined"], 4) for x in recs]} + rows.append(row) + print(json.dumps(row), flush=True) + scratch = ("/tmp/claude-1273071/-home-yurenh2-emm/" + "bb97201d-d9b5-46e6-82fc-981d37e7ab98/scratchpad/") + json.dump(rows, open(scratch + "control.json", "w")) + print(json.dumps({"done": True})) + + +if __name__ == "__main__": + main() diff --git a/logs/fgw_alpha_control.log b/logs/fgw_alpha_control.log new file mode 100644 index 0000000..de41a42 --- /dev/null +++ b/logs/fgw_alpha_control.log @@ -0,0 +1,4 @@ +{"trial": 0, "cost": "tlb", "mode": "cold", "E_truth": 0.339629709924146, "blind_refined": 0.0234375, "blind_alpha": 1.0, "blind_E": 0.6228801020056084, "oracle_refined": 0.08203125, "curve": [0.0781, 0.0234, 0.0117, 0.0273, 0.0742, 0.0352, 0.0664, 0.0352, 0.082, 0.0, 0.0547, 0.0234]} +{"trial": 0, "cost": "tlb", "mode": "warm", "E_truth": 0.339629709924146, "blind_refined": 0.84765625, "blind_alpha": 0.2, "blind_E": 0.339629709924146, "oracle_refined": 0.84765625, "curve": [0.0781, 0.8281, 0.8477, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398, 0.8398]} +{"trial": 0, "cost": "tlbperm", "mode": "cold", "E_truth": 0.339629709924146, "blind_refined": 0.02734375, "blind_alpha": 0.95, "blind_E": 0.6074469738905559, "oracle_refined": 0.078125, "curve": [0.0117, 0.0117, 0.043, 0.0781, 0.0078, 0.0352, 0.0156, 0.0156, 0.0039, 0.0078, 0.0273, 0.0234]} +{"trial": 0, "cost": "tlbperm", "mode": "warm", "E_truth": 0.339629709924146, "blind_refined": 0.03125, "blind_alpha": 1.0, "blind_E": 0.7038378812415739, "oracle_refined": 0.04296875, "curve": [0.0117, 0.043, 0.0234, 0.0234, 0.043, 0.0391, 0.0391, 0.0391, 0.0312, 0.0312, 0.0312, 0.0312]} diff --git a/logs/fgw_alpha_lean.log b/logs/fgw_alpha_lean.log new file mode 100644 index 0000000..0787289 --- /dev/null +++ b/logs/fgw_alpha_lean.log @@ -0,0 +1,71 @@ +# omit-size trial0 E_truth=0.3396 + cold a=0.0 raw=0.0625 ref=0.0781 E=0.6734 (8.4s) + cold a=0.2 raw=0.0117 ref=0.0117 E=0.6397 (9.5s) + cold a=0.4 raw=0.0703 ref=0.0742 E=0.6363 (2.8s) + cold a=0.6 raw=0.0586 ref=0.0664 E=0.6366 (1.2s) + cold a=0.75 raw=0.0508 ref=0.0430 E=0.6678 (5.7s) + cold a=0.85 raw=0.0312 ref=0.0234 E=0.6678 (1.6s) + cold a=0.88 raw=0.0586 ref=0.0508 E=0.6401 (2.4s) + cold a=0.9 raw=0.0000 ref=0.0000 E=0.6654 (2.2s) + cold a=0.92 raw=0.0352 ref=0.0312 E=0.7672 (0.5s) + cold a=0.94 raw=0.1680 ref=0.1875 E=0.6087 (0.5s) + cold a=0.95 raw=0.0625 ref=0.0547 E=0.6642 (2.1s) + cold a=0.96 raw=0.0664 ref=0.0664 E=0.6551 (0.3s) + cold a=0.97 raw=0.8203 ref=0.8203 E=0.3396 (0.4s) + cold a=0.98 raw=0.0039 ref=0.0039 E=0.6459 (1.0s) + cold a=0.99 raw=0.0156 ref=0.0156 E=0.6409 (1.0s) + cold a=1.0 raw=0.0195 ref=0.0234 E=0.6229 (0.7s) + == cold: blind-by-E refined=0.8203 (alpha=0.97), oracle-best=0.8203 + warm a=0.0 raw=0.0625 ref=0.0781 E=0.6734 (0.6s) + warm a=0.2 raw=0.8242 ref=0.8320 E=0.3396 (1.2s) + warm a=0.4 raw=0.8281 ref=0.8320 E=0.3396 (0.4s) + warm a=0.6 raw=0.8320 ref=0.8320 E=0.3396 (0.3s) + warm a=0.75 raw=0.8320 ref=0.8320 E=0.3396 (0.2s) + warm a=0.85 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.88 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.9 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.92 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.94 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.95 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.3s) + warm a=0.96 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.97 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.98 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.3s) + warm a=0.99 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=1.0 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.0s) + == warm: blind-by-E refined=0.8320 (alpha=0.2), oracle-best=0.8320 +# synth-full trial0 E_truth=0.1430 + cold a=0.0 raw=0.2578 ref=0.8633 E=0.2199 (2.1s) + cold a=0.2 raw=0.9531 ref=0.9609 E=0.1425 (2.6s) + cold a=0.4 raw=0.9531 ref=0.9688 E=0.1425 (1.7s) + cold a=0.6 raw=0.9570 ref=0.9766 E=0.1425 (2.2s) + cold a=0.75 raw=0.9375 ref=0.9609 E=0.1425 (3.0s) + cold a=0.85 raw=0.9531 ref=0.9609 E=0.1425 (4.3s) + cold a=0.88 raw=0.9453 ref=0.9609 E=0.1425 (3.8s) + cold a=0.9 raw=0.9766 ref=0.9844 E=0.1425 (3.5s) + cold a=0.92 raw=0.9453 ref=0.9531 E=0.1425 (0.9s) + cold a=0.94 raw=0.9531 ref=0.9766 E=0.1425 (2.0s) + cold a=0.95 raw=0.9453 ref=0.9531 E=0.1425 (1.1s) + cold a=0.96 raw=0.9453 ref=0.9531 E=0.1425 (0.8s) + cold a=0.97 raw=0.9375 ref=0.9531 E=0.1425 (6.2s) + cold a=0.98 raw=0.3789 ref=0.3789 E=0.5018 (7.1s) + cold a=0.99 raw=0.4219 ref=0.4141 E=0.4734 (10.8s) + cold a=1.0 raw=0.3633 ref=0.3906 E=0.5082 (12.6s) + == cold: blind-by-E refined=0.9609 (alpha=0.2), oracle-best=0.9844 + warm a=0.0 raw=0.2578 ref=0.8633 E=0.2199 (4.5s) + warm a=0.2 raw=0.9375 ref=0.9453 E=0.1425 (1.9s) + warm a=0.4 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (2.5s) + warm a=0.6 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.3s) + warm a=0.75 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (1.4s) + warm a=0.85 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.9s) + warm a=0.88 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.6s) + warm a=0.9 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.8s) + warm a=0.92 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (1.2s) + warm a=0.94 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.8s) + warm a=0.95 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.4s) + warm a=0.96 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.6s) + warm a=0.97 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.7s) + warm a=0.98 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.4s) + warm a=0.99 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.1s) + warm a=1.0 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.2s) + == warm: blind-by-E refined=0.9453 (alpha=0.2), oracle-best=0.9453 +{"done": true} diff --git a/logs/fgw_alpha_mech.log b/logs/fgw_alpha_mech.log new file mode 100644 index 0000000..0c3f85f --- /dev/null +++ b/logs/fgw_alpha_mech.log @@ -0,0 +1,12 @@ +--- trial 0 E_truth=0.339630 + tlb-vertex start=0.0586 | 1 power step ->0.0312 (E=0.7114) | 5 steps ->0.1719 (E=0.6409) | warm FGW ->0.8555 (E=0.3396) + random-vertex start=0.0000 | 1 power step ->0.0312 (E=0.7279) | 5 steps ->0.0781 (E=0.6819) | warm FGW ->0.8125 (E=0.3396) + grampa-vertex start=0.0156 | 1 power step ->0.0938 (E=0.7821) | 5 steps ->0.0469 (E=0.7049) | warm FGW ->0.8203 (E=0.3396) +--- trial 1 E_truth=0.339630 + tlb-vertex start=0.0508 | 1 power step ->0.0352 (E=0.7114) | 5 steps ->0.1836 (E=0.6409) | warm FGW ->0.8164 (E=0.3396) + random-vertex start=0.0039 | 1 power step ->0.0156 (E=0.7051) | 5 steps ->0.0195 (E=0.6692) | warm FGW ->0.0195 (E=0.6554) + grampa-vertex start=0.0195 | 1 power step ->0.0742 (E=0.7821) | 5 steps ->0.0547 (E=0.7049) | warm FGW ->0.8164 (E=0.3396) +--- trial 2 E_truth=0.339630 + tlb-vertex start=0.0430 | 1 power step ->0.0312 (E=0.7114) | 5 steps ->0.1797 (E=0.6409) | warm FGW ->0.8359 (E=0.3396) + random-vertex start=0.0000 | 1 power step ->0.0195 (E=0.7149) | 5 steps ->0.0117 (E=0.6700) | warm FGW ->0.0391 (E=0.6387) + grampa-vertex start=0.0156 | 1 power step ->0.0820 (E=0.7821) | 5 steps ->0.0508 (E=0.7049) | warm FGW ->0.8203 (E=0.3396) diff --git a/logs/lean.log b/logs/lean.log new file mode 100644 index 0000000..0787289 --- /dev/null +++ b/logs/lean.log @@ -0,0 +1,71 @@ +# omit-size trial0 E_truth=0.3396 + cold a=0.0 raw=0.0625 ref=0.0781 E=0.6734 (8.4s) + cold a=0.2 raw=0.0117 ref=0.0117 E=0.6397 (9.5s) + cold a=0.4 raw=0.0703 ref=0.0742 E=0.6363 (2.8s) + cold a=0.6 raw=0.0586 ref=0.0664 E=0.6366 (1.2s) + cold a=0.75 raw=0.0508 ref=0.0430 E=0.6678 (5.7s) + cold a=0.85 raw=0.0312 ref=0.0234 E=0.6678 (1.6s) + cold a=0.88 raw=0.0586 ref=0.0508 E=0.6401 (2.4s) + cold a=0.9 raw=0.0000 ref=0.0000 E=0.6654 (2.2s) + cold a=0.92 raw=0.0352 ref=0.0312 E=0.7672 (0.5s) + cold a=0.94 raw=0.1680 ref=0.1875 E=0.6087 (0.5s) + cold a=0.95 raw=0.0625 ref=0.0547 E=0.6642 (2.1s) + cold a=0.96 raw=0.0664 ref=0.0664 E=0.6551 (0.3s) + cold a=0.97 raw=0.8203 ref=0.8203 E=0.3396 (0.4s) + cold a=0.98 raw=0.0039 ref=0.0039 E=0.6459 (1.0s) + cold a=0.99 raw=0.0156 ref=0.0156 E=0.6409 (1.0s) + cold a=1.0 raw=0.0195 ref=0.0234 E=0.6229 (0.7s) + == cold: blind-by-E refined=0.8203 (alpha=0.97), oracle-best=0.8203 + warm a=0.0 raw=0.0625 ref=0.0781 E=0.6734 (0.6s) + warm a=0.2 raw=0.8242 ref=0.8320 E=0.3396 (1.2s) + warm a=0.4 raw=0.8281 ref=0.8320 E=0.3396 (0.4s) + warm a=0.6 raw=0.8320 ref=0.8320 E=0.3396 (0.3s) + warm a=0.75 raw=0.8320 ref=0.8320 E=0.3396 (0.2s) + warm a=0.85 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.88 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.9 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.92 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.94 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.95 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.3s) + warm a=0.96 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.2s) + warm a=0.97 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=0.98 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.3s) + warm a=0.99 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.1s) + warm a=1.0 raw=0.8320 ref=0.8320 E=0.3396 FROZEN (0.0s) + == warm: blind-by-E refined=0.8320 (alpha=0.2), oracle-best=0.8320 +# synth-full trial0 E_truth=0.1430 + cold a=0.0 raw=0.2578 ref=0.8633 E=0.2199 (2.1s) + cold a=0.2 raw=0.9531 ref=0.9609 E=0.1425 (2.6s) + cold a=0.4 raw=0.9531 ref=0.9688 E=0.1425 (1.7s) + cold a=0.6 raw=0.9570 ref=0.9766 E=0.1425 (2.2s) + cold a=0.75 raw=0.9375 ref=0.9609 E=0.1425 (3.0s) + cold a=0.85 raw=0.9531 ref=0.9609 E=0.1425 (4.3s) + cold a=0.88 raw=0.9453 ref=0.9609 E=0.1425 (3.8s) + cold a=0.9 raw=0.9766 ref=0.9844 E=0.1425 (3.5s) + cold a=0.92 raw=0.9453 ref=0.9531 E=0.1425 (0.9s) + cold a=0.94 raw=0.9531 ref=0.9766 E=0.1425 (2.0s) + cold a=0.95 raw=0.9453 ref=0.9531 E=0.1425 (1.1s) + cold a=0.96 raw=0.9453 ref=0.9531 E=0.1425 (0.8s) + cold a=0.97 raw=0.9375 ref=0.9531 E=0.1425 (6.2s) + cold a=0.98 raw=0.3789 ref=0.3789 E=0.5018 (7.1s) + cold a=0.99 raw=0.4219 ref=0.4141 E=0.4734 (10.8s) + cold a=1.0 raw=0.3633 ref=0.3906 E=0.5082 (12.6s) + == cold: blind-by-E refined=0.9609 (alpha=0.2), oracle-best=0.9844 + warm a=0.0 raw=0.2578 ref=0.8633 E=0.2199 (4.5s) + warm a=0.2 raw=0.9375 ref=0.9453 E=0.1425 (1.9s) + warm a=0.4 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (2.5s) + warm a=0.6 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.3s) + warm a=0.75 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (1.4s) + warm a=0.85 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.9s) + warm a=0.88 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.6s) + warm a=0.9 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.8s) + warm a=0.92 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (1.2s) + warm a=0.94 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.8s) + warm a=0.95 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.4s) + warm a=0.96 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.6s) + warm a=0.97 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.7s) + warm a=0.98 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.4s) + warm a=0.99 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.1s) + warm a=1.0 raw=0.9375 ref=0.9453 E=0.1425 FROZEN (0.2s) + == warm: blind-by-E refined=0.9453 (alpha=0.2), oracle-best=0.9453 +{"done": true} diff --git a/logs/lean_probe.py b/logs/lean_probe.py new file mode 100644 index 0000000..f3b394c --- /dev/null +++ b/logs/lean_probe.py @@ -0,0 +1,84 @@ +"""Lean version: 15-alpha grid, TLB cost, cold vs warm, omit-size then synth-full.""" +from __future__ import annotations +import json, sys, time +sys.path.insert(0, "/home/yurenh2/emm") +import numpy as np, ot, torch +from scipy.optimize import linear_sum_assignment +from worldalign.synth_fast_gate import ClosedFormEnergy, fast_pair_descent +from worldalign.synth_triangle_gate import standardized + +DEV = torch.device("cuda:1") +ALPHAS = [0.0, 0.2, 0.4, 0.6, 0.75, 0.85, 0.88, 0.90, 0.92, 0.94, 0.95, + 0.96, 0.97, 0.98, 0.99, 1.0] + + +def standardise(m): + mask = ~np.eye(len(m), dtype=bool); v = m[mask] + o = (m - v.mean()) / v.std(); np.fill_diagonal(o, 0.0); return o + + +def tlb(V, S): + n = len(V); a = np.sort(V, 1); b = np.sort(S, 1) + return (a * a).sum(1)[:, None] / n + (b * b).sum(1)[None, :] / n - 2 * (a @ b.T) / n + + +def run(path, label, trial=0): + st = torch.load(path, map_location="cpu", weights_only=False) + V = standardise(st["visual_field"].double().numpy()) + T = standardise(st["text_field"].double().numpy()) + n = len(V); w = ot.unif(n) + rng = np.random.default_rng(trial); hid = rng.permutation(n) + S = T[np.ix_(hid, hid)] + Vg = standardized(torch.from_numpy(V).to(DEV)).double() + Sg = standardized(torch.from_numpy(S).to(DEV)).double() + en = ClosedFormEnergy(Sg, Vg, 1.0, 0.0, 64) + M = tlb(V, S); M = M / M.max() + + def acc(c): return float((hid[c] == np.arange(n)).mean()) + + def refine(c): + f = fast_pair_descent(Sg, Vg, torch.from_numpy(np.ascontiguousarray(c)).to(DEV), + 600).cpu().numpy() + e = float(en.energy(torch.from_numpy(np.ascontiguousarray(f)).to(DEV)[None])[0]) + return e, acc(f) + + et = float(en.energy(torch.from_numpy(np.ascontiguousarray(np.argsort(hid))).to(DEV)[None])[0]) + print(f"# {label} trial{trial} E_truth={et:.4f}", flush=True) + + out = {"label": label, "E_truth": et, "cold": [], "warm": []} + for mode in ("cold", "warm"): + G = np.outer(w, w) + prev = None + for a in ALPHAS: + t0 = time.time() + g0 = None if mode == "cold" else G + G2 = ot.gromov.fused_gromov_wasserstein(M, V, S, w, w, "square_loss", + alpha=a, G0=g0, max_iter=200, + tol_rel=1e-9) + if mode == "warm": + G = G2 + _, cols = linear_sum_assignment(-G2) + e, ar = refine(cols) + same = (prev is not None and np.array_equal(cols, prev)) + prev = cols + rec = {"alpha": a, "raw": acc(cols), "refined": ar, "E": e, + "sec": round(time.time() - t0, 1), "same_as_prev": bool(same)} + out[mode].append(rec) + print(f" {mode} a={a:<5} raw={rec['raw']:.4f} ref={ar:.4f} E={e:.4f}" + f" {'FROZEN' if same else ''} ({rec['sec']}s)", flush=True) + good = out[mode] + best = min(good, key=lambda x: x["E"]) + print(f" == {mode}: blind-by-E refined={best['refined']:.4f} (alpha={best['alpha']})," + f" oracle-best={max(x['refined'] for x in good):.4f}", flush=True) + return out + + +if __name__ == "__main__": + base = "/home/yurenh2/emm/artifacts/synth_v1/" + res = [] + for p, l in [("omit_size.pt", "omit-size"), ("omit_none.pt", "synth-full")]: + res.append(run(base + p, l)) + scratch = ("/tmp/claude-1273071/-home-yurenh2-emm/" + "bb97201d-d9b5-46e6-82fc-981d37e7ab98/scratchpad/") + json.dump(res, open(scratch + "lean.json", "w")) + print(json.dumps({"done": True})) diff --git a/logs/mech.log b/logs/mech.log new file mode 100644 index 0000000..817501a --- /dev/null +++ b/logs/mech.log @@ -0,0 +1,11 @@ +--- trial 0 E_truth=0.339630 + tlb-vertex start=0.0586 | 1 power step ->0.0312 (E=0.7114) | 5 steps ->0.1719 (E=0.6409) | warm FGW ->0.8555 (E=0.3396) + random-vertex start=0.0000 | 1 power step ->0.0312 (E=0.7279) | 5 steps ->0.0781 (E=0.6819) | warm FGW ->0.8125 (E=0.3396) + grampa-vertex start=0.0156 | 1 power step ->0.0938 (E=0.7821) | 5 steps ->0.0469 (E=0.7049) | warm FGW ->0.8203 (E=0.3396) +--- trial 1 E_truth=0.339630 + tlb-vertex start=0.0508 | 1 power step ->0.0352 (E=0.7114) | 5 steps ->0.1836 (E=0.6409) | warm FGW ->0.8164 (E=0.3396) + random-vertex start=0.0039 | 1 power step ->0.0156 (E=0.7051) | 5 steps ->0.0195 (E=0.6692) | warm FGW ->0.0195 (E=0.6554) + grampa-vertex start=0.0195 | 1 power step ->0.0742 (E=0.7821) | 5 steps ->0.0547 (E=0.7049) | warm FGW ->0.8164 (E=0.3396) +--- trial 2 E_truth=0.339630 + tlb-vertex start=0.0430 | 1 power step ->0.0312 (E=0.7114) | 5 steps ->0.1797 (E=0.6409) | warm FGW ->0.8359 (E=0.3396) + random-vertex start=0.0000 | 1 power step ->0.0195 (E=0.7149) | 5 steps ->0.0117 (E=0.6700) | warm FGW ->0.0391 (E=0.6387) diff --git a/logs/mech.py b/logs/mech.py new file mode 100644 index 0000000..6762777 --- /dev/null +++ b/logs/mech.py @@ -0,0 +1,51 @@ +import sys, numpy as np, ot, torch +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from worldalign.synth_fast_gate import ClosedFormEnergy, fast_pair_descent +from worldalign.synth_triangle_gate import standardized +DEV=torch.device("cuda:1") +def standardise(m): + mask=~np.eye(len(m),dtype=bool); v=m[mask] + o=(m-v.mean())/v.std(); np.fill_diagonal(o,0.0); return o +def tlb(V,S): + n=len(V); a=np.sort(V,1); b=np.sort(S,1) + return (a*a).sum(1)[:,None]/n+(b*b).sum(1)[None,:]/n-2*(a@b.T)/n +st=torch.load("/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt",map_location='cpu',weights_only=False) +V=standardise(st["visual_field"].double().numpy()); T=standardise(st["text_field"].double().numpy()) +n=len(V); w=ot.unif(n) +for trial in range(3): + rng=np.random.default_rng(trial); hid=rng.permutation(n); S=T[np.ix_(hid,hid)] + Vg=standardized(torch.from_numpy(V).to(DEV)).double(); Sg=standardized(torch.from_numpy(S).to(DEV)).double() + en=ClosedFormEnergy(Sg,Vg,1.0,0.0,64) + Et=float(en.energy(torch.from_numpy(np.ascontiguousarray(np.argsort(hid))).to(DEV)[None])[0]) + acc=lambda c: float((hid[c]==np.arange(n)).mean()) + def refine(c,it=2000): + f=fast_pair_descent(Sg,Vg,torch.from_numpy(np.ascontiguousarray(c)).to(DEV),it).cpu().numpy() + return float(en.energy(torch.from_numpy(np.ascontiguousarray(f)).to(DEV)[None])[0]), acc(f) + M=tlb(V,S); M=M/M.max() + _,tlbcols=linear_sum_assignment(M) + def vert(cols): + P=np.zeros((n,n)); P[np.arange(n),cols]=1.0/n; return P + starts={"tlb-vertex":tlbcols, + "random-vertex":np.random.default_rng(50+trial).permutation(n), + "grampa-vertex":None} + from worldalign.spectral_match import grampa + starts["grampa-vertex"]=grampa(V,S,1.0) + print(f"--- trial {trial} E_truth={Et:.6f}") + for name,cols in starts.items(): + # (i) pure power step: LSAP on gradient V P S (anchor-propagation, one round) + P=vert(cols); grad=V@P@S + _,c1=linear_sum_assignment(-grad) + e1,a1=refine(c1) + # (ii) 5 rounds of power step + c=cols.copy() + for _ in range(5): + P=vert(c); _,c=linear_sum_assignment(-(V@P@S)) + e2,a2=refine(c) + # (iii) warm FGW at alpha=0.1 then 0.2 from this vertex + G=vert(cols) + for a in (0.1,0.2,0.5): + G=ot.gromov.fused_gromov_wasserstein(M,V,S,w,w,"square_loss",alpha=a,G0=G,max_iter=200,tol_rel=1e-9) + _,c3=linear_sum_assignment(-G); e3,a3=refine(c3) + print(" %-14s start=%.4f | 1 power step ->%.4f (E=%.4f) | 5 steps ->%.4f (E=%.4f) | warm FGW ->%.4f (E=%.4f)" + %(name,acc(cols),a1,e1,a2,e2,a3,e3), flush=True) diff --git a/logs/verify.py b/logs/verify.py new file mode 100644 index 0000000..17fd643 --- /dev/null +++ b/logs/verify.py @@ -0,0 +1,59 @@ +import sys, numpy as np, ot, torch +sys.path.insert(0,'/home/yurenh2/emm') +from scipy.optimize import linear_sum_assignment +from worldalign.synth_fast_gate import ClosedFormEnergy, fast_pair_descent +from worldalign.synth_triangle_gate import standardized +DEV=torch.device("cuda:1") +def standardise(m): + mask=~np.eye(len(m),dtype=bool); v=m[mask] + o=(m-v.mean())/v.std(); np.fill_diagonal(o,0.0); return o +def tlb(V,S): + n=len(V); a=np.sort(V,1); b=np.sort(S,1) + return (a*a).sum(1)[:,None]/n+(b*b).sum(1)[None,:]/n-2*(a@b.T)/n +st=torch.load("/home/yurenh2/emm/artifacts/synth_v1/omit_size.pt",map_location='cpu',weights_only=False) +V=standardise(st["visual_field"].double().numpy()); T=standardise(st["text_field"].double().numpy()) +n=len(V); w=ot.unif(n) +rng=np.random.default_rng(0); hid=rng.permutation(n); S=T[np.ix_(hid,hid)] +Vg=standardized(torch.from_numpy(V).to(DEV)).double(); Sg=standardized(torch.from_numpy(S).to(DEV)).double() +en=ClosedFormEnergy(Sg,Vg,1.0,0.0,64) +M=tlb(V,S); M=M/M.max() +G=np.outer(w,w) +for a in (0.0,0.2): + G=ot.gromov.fused_gromov_wasserstein(M,V,S,w,w,"square_loss",alpha=a,G0=G,max_iter=200,tol_rel=1e-9) +_,cols=linear_sum_assignment(-G) +fin=fast_pair_descent(Sg,Vg,torch.from_numpy(np.ascontiguousarray(cols)).to(DEV),2000).cpu().numpy() +truth=np.argsort(hid) +Ef=float(en.energy(torch.from_numpy(np.ascontiguousarray(fin)).to(DEV)[None])[0]) +Et=float(en.energy(torch.from_numpy(np.ascontiguousarray(truth)).to(DEV)[None])[0]) +acc=float((hid[fin]==np.arange(n)).mean()) +print("accuracy=%.4f E_found=%.10f E_truth=%.10f diff=%.3e"%(acc,Ef,Et,Ef-Et)) +wrong=np.where(hid[fin]!=np.arange(n))[0] +print("n_wrong=%d"%len(wrong)) +# for each wrong scene i: correlation between S-row of assigned node and S-row of true node +Sal=S # in shuffled index space +corr=[] +vcorr=[] +for i in wrong: + a_idx=fin[i]; t_idx=truth[i] + x=np.delete(Sal[a_idx],[a_idx,t_idx]); y=np.delete(Sal[t_idx],[a_idx,t_idx]) + corr.append(np.corrcoef(x,y)[0,1]) + u=np.delete(V[i],[i]); + # visual side: is the visually-assigned scene similar to i? compare V rows of i and of the scene truly at a_idx + j=hid[a_idx] + v1=np.delete(V[i],[i,j]); v2=np.delete(V[j],[i,j]) + vcorr.append(np.corrcoef(v1,v2)[0,1]) +rng2=np.random.default_rng(7); pairs=rng2.integers(0,n,(len(wrong),2)) +base=[np.corrcoef(np.delete(Sal[p],[p,q]),np.delete(Sal[q],[p,q]))[0,1] for p,q in pairs if p!=q] +print("text-row corr of confused partners: median=%.3f mean=%.3f | random pair baseline median=%.3f"%(np.median(corr),np.mean(corr),np.median(base))) +print("visual-row corr of confused partners: median=%.3f"%np.median(vcorr)) +# cycle structure of the error +perm=hid[fin] +seen=set(); cyc=[] +for i in range(n): + if i in seen or perm[i]==i: continue + c=0; j=i + while j not in seen: + seen.add(j); j=perm[j]; c+=1 + cyc.append(c) +import collections +print("error cycle lengths:", dict(collections.Counter(cyc))) |