1import json, math
  2from pathlib import Path
  3import numpy as np
  4from sklearn.linear_model import LogisticRegression
  5from sklearn.model_selection import train_test_split
  6from sklearn.pipeline import make_pipeline
  7from sklearn.preprocessing import StandardScaler
  8from sklearn.metrics import accuracy_score
  9
 10SEED = 1461
 11rng = np.random.default_rng(SEED)
 12
 13
 14def laplacian_path(n):
 15    A = np.zeros((n, n))
 16    for i in range(n - 1):
 17        A[i, i + 1] = A[i + 1, i] = 1.0
 18    return np.diag(A.sum(1)) - A
 19
 20
 21def spectral_data(q):
 22    n = len(q)
 23    H = laplacian_path(n) + np.diag(q)
 24    lam, phi = np.linalg.eigh(H)
 25    S = phi ** 2
 26    sums = np.array([lam[i] + lam[j] for i in range(n) for j in range(i, n)])
 27    diffs = np.array([lam[i] - lam[j] for i in range(n) for j in range(i)])
 28    gap_sum = np.min(np.abs(sums[:, None] - sums[None, :] + np.eye(len(sums)) * 1e9)) if len(sums) > 1 else np.inf
 29    gap_diff = np.min(np.abs(diffs[:, None] - diffs[None, :] + np.eye(len(diffs)) * 1e9)) if len(diffs) > 1 else np.inf
 30    overlap = min(np.max(np.abs(phi[:, i] * phi[:, j])) for i in range(n) for j in range(i))
 31    return H, lam, phi, S, float(gap_sum), float(gap_diff), float(np.linalg.svd(S, compute_uv=False)[-1]), float(overlap)
 32
 33
 34def observations(phi, lam, u, times):
 35    c = phi.T @ u
 36    return np.array([np.abs(phi @ (np.exp(-1j * t * lam) * c)) for t in times]).ravel()
 37
 38
 39def phase_aligned_jacobian(phi, lam, u, times):
 40    # finite differences in 2n real coordinates; remove the global-phase tangent
 41    n = len(u)
 42    x = np.r_[u.real, u.imag]
 43    f0 = observations(phi, lam, u, times)
 44    J = np.empty((len(f0), 2*n))
 45    eps = 2e-6
 46    for k in range(2*n):
 47        xp = x.copy(); xp[k] += eps
 48        xm = x.copy(); xm[k] -= eps
 49        J[:, k] = (observations(phi, lam, xp[:n] + 1j*xp[n:], times) - observations(phi, lam, xm[:n] + 1j*xm[n:], times)) / (2*eps)
 50    # tangent to global phase is (-Im u, Re u), which is exactly a null direction
 51    tangent = np.r_[-u.imag, u.real]
 52    tangent /= np.linalg.norm(tangent)
 53    Q, _ = np.linalg.qr(np.column_stack([tangent, rng.normal(size=(2*n, 2*n-1))]))
 54    B = Q[:, 1:]
 55    # QR above may not preserve tangent as first column due nonorthogonal random columns;
 56    # use a robust null-space basis from SVD instead.
 57    _, _, vh = np.linalg.svd(tangent[None, :])
 58    B = vh[1:].T
 59    return np.linalg.svd(J @ B, compute_uv=False)
 60
 61
 62class SpectralPhaselessEncoder:
 63    """Fixed graph Schrodinger evolution exposed only through coordinate magnitudes."""
 64    def __init__(self, graph_laplacian, potential, dt=1.0, steps=6):
 65        self.H = np.asarray(graph_laplacian, float) + np.diag(np.asarray(potential, float))
 66        self.lam, self.phi = np.linalg.eigh(self.H)
 67        self.dt, self.steps = float(dt), int(steps)
 68        self.U = self.phi @ np.diag(np.exp(-1j * self.dt * self.lam)) @ self.phi.T
 69
 70    def transform(self, h):
 71        x = np.asarray(h, complex).copy()
 72        out = []
 73        for _ in range(self.steps):
 74            out.append(np.abs(x))
 75            x = self.U @ x
 76        return np.concatenate(out)
 77
 78    def diagnostics(self):
 79        n = len(self.lam)
 80        sums = np.array([self.lam[i]+self.lam[j] for i in range(n) for j in range(i,n)])
 81        gap = np.min(np.abs(sums[:,None]-sums[None,:] + np.eye(len(sums))*1e9))
 82        S = self.phi**2
 83        overlap = min(np.max(np.abs(self.phi[:,i]*self.phi[:,j])) for i in range(n) for j in range(i))
 84        return {"sigma_min_S": float(np.linalg.svd(S, compute_uv=False)[-1]),
 85                "pair_sum_gap": float(gap), "minimum_pair_overlap": float(overlap)}
 86
 87
 88def math_checks(n=6):
 89    rows = []
 90    for scale in [0.02, 0.1, 0.3, 1.0, 3.0]:
 91        vals = []
 92        for _ in range(40):
 93            q = scale * rng.normal(size=n)
 94            *_, gap_sum, gap_diff, smin, overlap = spectral_data(q)
 95            vals.append((gap_sum, gap_diff, smin, overlap))
 96        a = np.array(vals)
 97        rows.append({"potential_scale": scale, "median_pair_sum_gap": float(np.median(a[:,0])),
 98                     "median_difference_gap": float(np.median(a[:,1])),
 99                     "median_sigma_min_S": float(np.median(a[:,2])),
100                     "median_overlap": float(np.median(a[:,3]))})
101    q = rng.normal(size=n)
102    _, lam, phi, S, gs, gd, sm, ov = spectral_data(q)
103    u = rng.normal(size=n) + 1j*rng.normal(size=n)
104    ranks = []
105    mins = []
106    for K in range(1, 9):
107        times = np.linspace(0, 7.0, K)
108        sv = phase_aligned_jacobian(phi, lam, u, times)
109        ranks.append(int(np.sum(sv > 1e-5 * sv[0])))
110        mins.append(float(sv[-1]))
111    # Noise amplification: perturb observations and solve a local least-squares correction.
112    noise_rows = []
113    for K in [1, 2, 3, 4, 6, 8]:
114        times = np.linspace(0, 7.0, K)
115        sv = phase_aligned_jacobian(phi, lam, u, times)
116        noise_rows.append({"K": K, "jacobian_rank": int(np.sum(sv > 1e-5*sv[0])),
117                           "smallest_nonphase_singular": float(sv[-1]),
118                           "predicted_noise_gain_1_over_sigma": float(1/sv[-1])})
119    return {"selected_operator": {"eigenvalues": lam.tolist(), "pair_sum_gap": gs,
120             "difference_gap": gd, "sigma_min_S": sm, "overlap_min": ov},
121            "conditioning_sweep": rows, "time_sample_sweep": {"K": list(range(1,9)),
122             "local_rank": ranks, "smallest_singular": mins}, "noise_sweep": noise_rows, "global_phase_max_abs_error": float(max(np.max(np.abs(observations(phi, lam, u, np.linspace(0,7,6)) - observations(phi, lam, u*np.exp(1j*1.234), np.linspace(0,7,6)))), 0.0))}
123
124
125def classification(n=6, N=2400):
126    q = rng.normal(size=n)
127    _, lam, phi, *_ = spectral_data(q)
128    # Label is a phase-invariant spectral energy contrast; random global phases cannot affect it.
129    U = rng.normal(size=(N,n)) + 1j*rng.normal(size=(N,n))
130    c = U @ phi
131    label = ((np.abs(c[:,0])**2 + np.abs(c[:,1])**2) > (np.abs(c[:,2])**2 + np.abs(c[:,3])**2)).astype(int)
132    phase = rng.uniform(0, 2*np.pi, N)
133    U *= np.exp(1j*phase)[:,None]
134    def enc(K):
135        ts = np.linspace(0, 7.0, K)
136        return np.concatenate([np.abs((phi @ (np.exp(-1j*t*lam)[:,None] * (U @ phi).T)).T) for t in ts], axis=1)
137    Xidea = enc(6)
138    Xbase = np.abs(U)
139    tr, te = train_test_split(np.arange(N), test_size=.3, random_state=SEED, stratify=label)
140    def fit(X):
141        clf = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, random_state=SEED))
142        clf.fit(X[tr], label[tr]); return accuracy_score(label[te], clf.predict(X[te]))
143    # Add sensor noise at test time, matching the requested robustness check.
144    noise = .05 * rng.normal(size=Xidea[te].shape)
145    return {"clean_accuracy": {"phaseless_K6": float(fit(Xidea)), "single_time_coordinate_magnitude": float(fit(Xbase))},
146            "noisy_test_accuracy_phaseless_K6": float(accuracy_score(label[te], make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, random_state=SEED)).fit(Xidea[tr], label[tr]).predict(Xidea[te] + noise))),
147            "N": N, "train_test": [len(tr), len(te)]}
148
149
150if __name__ == "__main__":
151    out = {"seed": SEED, "math": math_checks(), "classification": classification()}
152    Path("results.json").write_text(json.dumps(out, indent=2))
153    print(json.dumps(out, indent=2))