"""Numerical verification of the closed forms in Sections IV and V of paper 11 (mask-as-key encryption and jamming robustness). Every claim that enters the manuscript is checked here against Monte Carlo, with a PASS/FAIL verdict and the achieved agreement level printed. Real-vector convention, dimension d, U users, codebook of V unit-norm codewords, per-user masks with zero-mean entries normalized to ||m||^2 = d (so E[m_k^2] = 1). Notation matches the tex: y = (1/c) sum_u e_{s_u} .* m_u transmit frame legit score z_{u,i} = r_u^T (e_i .* m_u), r_u = y + n_u/h_u eve score zE_{u,i} = (y_E/h_E)^T (e_i .* mtil_u) jammer adds h_J sqrt(rho) w to the victim observation Run on CPU (NumPy); no training involved, pure algebra checks. """ from __future__ import annotations import numpy as np from pathlib import Path RNG = np.random.default_rng(2026) D, U, V = 64, 4, 256 CKPT = Path(__file__).resolve().parent.parent / "data" / "model_main.pt" def cached_main_model(): """The trained main-configuration model, from a checkpoint. V8 and V9 read the trained codebook. Retraining it reproduces only on the device that trained it, so a CPU run of the released package disagreed with the shipped numbers. The checkpoint fixes the codebook, which is what both checks are about; delete it to retrain. """ import torch from exp_full import main_model m = main_model() if CKPT.exists(): m.load_state_dict(torch.load(CKPT, map_location="cpu")) else: torch.save({k: v.cpu() for k, v in m.state_dict().items()}, CKPT) return m def unit_codebook(V, d, rng): E = rng.standard_normal((V, d)) return E / np.linalg.norm(E, axis=1, keepdims=True) def masks(U, d, rng): """Zero-mean entries normalized so ||m_u||^2 = d.""" M = rng.standard_normal((U, d)) return M / np.linalg.norm(M, axis=1, keepdims=True) * np.sqrt(d) ROWS = [] # (tag, claim, emp, err, tol, verdict) def report(tag, claim, emp, tol, extra=""): err = abs(claim - emp) ok = err <= tol ROWS.append((tag, claim, emp, err, tol, "PASS" if ok else "FAIL")) print(f"[{'PASS' if ok else 'FAIL'}] {tag}: claim={claim:.5g} " f"emp={emp:.5g} |err|={err:.2g} tol={tol:g} {extra}") return ok def v1_legit_self_alignment(): """Claim: E[e_s^T diag(m^2) e_s] = 1 (signal self-correlation).""" E = unit_codebook(V, D, RNG) vals = [] for _ in range(4000): m = masks(1, D, RNG)[0] s = RNG.integers(V) vals.append(float((E[s] ** 2) @ (m ** 2))) return report("V1 legit self-alignment", 1.0, float(np.mean(vals)), 2e-2) def v2_eve_uninformed(): """Claim: with an independent substitute mask, the eavesdropper's correct-index correlation has the same mean as any wrong index, so the mean advantage is zero and the eavesdropper SER = (V-1)/V, independent of SNR.""" E = unit_codebook(V, D, RNG) # mean advantage of the true index over the wrong indices, noiseless adv = [] ser_by_snr = {} for snr_db in [0.0, 10.0, 20.0, 80.0]: # 80 dB stands in for noiseless sigma = np.sqrt(1.0 / (D * 10 ** (snr_db / 10.0))) err = 0 trials = 6000 for _ in range(trials): s = RNG.integers(V, size=U) M = masks(U, D, RNG) c = 1.0 # scale-invariant for argmax y = np.zeros(D) for u in range(U): y += E[s[u]] * M[u] y /= np.sqrt(U) # any fixed scale # eavesdropper targets user 0 with an independent wrong mask mtil = masks(1, D, RNG)[0] hE = np.sqrt(-np.log(RNG.random())) rE = y + (sigma / hE) * RNG.standard_normal(D) scores = (E * mtil) @ rE # (V,) if snr_db == 80.0: adv.append(scores[s[0]] - scores.mean()) if scores.argmax() != s[0]: err += 1 ser_by_snr[snr_db] = err / trials chance = (V - 1) / V ok1 = report("V2a eve mean advantage", 0.0, float(np.mean(adv)), 3e-3) ok2 = True for snr_db, ser in ser_by_snr.items(): tag = f"V2b eve SER @ {int(snr_db)}dB" ok2 &= report(tag, chance, ser, 1.5e-2) return ok1 and ok2 def v3_leakage_vs_correlation(): """Claim: if the substitute mask has normalized correlation rho = /d with the true mask, the eavesdropper's true-index bias grows linearly in rho; independent random masks give E|rho| = O(1/sqrt(d)); orthogonal masks give rho = 0.""" E = unit_codebook(V, D, RNG) # (a) bias vs prescribed rho slopes = [] for rho in [0.0, 0.25, 0.5, 0.75, 1.0]: bias = [] for _ in range(3000): m = masks(1, D, RNG)[0] mp = masks(1, D, RNG)[0] mp = mp - (mp @ m) / (m @ m) * m # orthogonalize mp = mp / np.linalg.norm(mp) * np.sqrt(D) mtil = rho * m + np.sqrt(1 - rho ** 2) * mp s = RNG.integers(V) # noiseless single-user useful alignment for the true index bias.append(float((E[s] ** 2) @ (m * mtil))) slopes.append((rho, float(np.mean(bias)))) # claim: bias(rho) = rho * bias(1); check linearity b1 = slopes[-1][1] lin_ok = all(abs(b - rho * b1) <= 3e-2 for rho, b in slopes) print(f"[{'PASS' if lin_ok else 'FAIL'}] V3a bias linear in rho: " + ", ".join(f"rho={r:.2f}->{b:.3f}" for r, b in slopes)) ROWS.append(("V3a bias slope in kappa", "1.0", "%.4f" % b1, "%.4f" % abs(1.0 - b1), "0.03", "PASS" if lin_ok else "FAIL")) # (b) random independent mask correlation: E|corr| = sqrt(2/(pi d)) # (the folded-normal mean of a N(0, 1/d) variable) corrs = [] for _ in range(5000): m = masks(1, D, RNG)[0] mt = masks(1, D, RNG)[0] corrs.append(abs((m @ mt) / D)) emp = float(np.mean(corrs)) claim = float(np.sqrt(2.0 / (np.pi * D))) ok_b = report("V3b random mask E|corr|", claim, emp, 0.03 * claim) return lin_ok and ok_b def v4_blind_jammer_spread(): """Claim: a mask-blind jammer (w independent of m_u) contributes a zero-mean term to every candidate score with variance (hJ^2/hu^2) rho * sum_k w_k^2 e_{i,k}^2, i.e. it is spread with no systematic bias toward any index.""" E = unit_codebook(V, D, RNG) rho = 1.0 # (a) structural claim: the mask projection g_i = (w .* m)^T e_i is # zero-mean, decoupled from the positive gain factor hJ/hu. g, var_emp, var_cl = [], [], [] for _ in range(20000): m = masks(1, D, RNG)[0] w = RNG.standard_normal(D); w /= np.linalg.norm(w) i = RNG.integers(V) gi = (w * m) @ E[i] g.append(gi) var_emp.append(gi ** 2) var_cl.append(np.sum(w ** 2 * E[i] ** 2)) ok1 = report("V4a blind jammer projection mean", 0.0, float(np.mean(g)), 3e-3) # (b) variance of the projection matches sum_k w_k^2 e_{i,k}^2; the # full contribution scales this by (hJ^2/hu^2) rho. ok2 = report("V4b blind jammer projection variance", float(np.mean(var_cl)), float(np.mean(var_emp)), 0.01 * float(np.mean(var_cl))) return ok1 and ok2 def v5_matched_jammer_concentrates(): """Claim: a mask-matched jammer aligned with the victim key for a target index t creates a bias of order one on index t, while the blind-jammer projection has zero mean and RMS of order 1/sqrt(d). The physically meaningful separation is matched bias over blind RMS, which is sqrt(d) (the sample mean of the blind bias estimates zero and is pure Monte Carlo noise, so it is NOT a valid denominator).""" E = unit_codebook(V, D, RNG) bias_matched, blind_sq = [], [] for _ in range(3000): m = masks(1, D, RNG)[0] t = RNG.integers(V) wm = E[t] * m; wm /= np.linalg.norm(wm) # matched (needs m) wb = RNG.standard_normal(D); wb /= np.linalg.norm(wb) # blind bias_matched.append(float((wm * m) @ E[t])) blind_sq.append(float(((wb * m) @ E[t]) ** 2)) bm = float(np.mean(bias_matched)) brms = float(np.sqrt(np.mean(blind_sq))) ratio = bm / brms ok = abs(ratio - np.sqrt(D)) <= 0.25 * np.sqrt(D) and bm > 0.9 print(f"[{'PASS' if ok else 'FAIL'}] V5 matched bias / blind RMS: " f"matched={bm:.3f} blind_rms={brms:.4f} ratio={ratio:.1f} " f"(claim sqrt(d)={np.sqrt(D):.1f})") ROWS.append(("V5 matched bias over blind RMS", "%.1f" % np.sqrt(D), "%.1f" % ratio, "%.2f" % abs(ratio - np.sqrt(D)), "1.0", "PASS" if ok else "FAIL")) return ok def v6_coded_oma_outage(): """The genie reference the manuscript concedes: an OMA user coding at the Rayleigh outage limit of its own allocation. The user owns L = d/U real dimensions, so d/(2U) complex uses, and carries log2(V) bits. Its per-dimension SNR is the frame SNR, the same convention snr_to_sigma2 sets. Outage is the probability that the instantaneous mutual information falls below that rate.""" import math U, d, V = 4, 256, 65536 for snr_db in (10.0,): g = 10.0 ** (snr_db / 10.0) uses = d / (2.0 * U) # complex channel uses rate = math.log2(V) / uses # bits per complex use thr = (2.0 ** rate - 1.0) / g # |h|^2 threshold pout = 1.0 - math.exp(-thr) # |h|^2 ~ Exp(1) print(f"V6 coded-OMA outage @ {snr_db:.0f} dB: {pout:.4f} " f"(rate {rate:.3f} bit/use)") ROWS.append(("V6 coded-OMA outage @ %.0f dB" % snr_db, "%.4f" % pout, "%.6f" % pout, "0", "0", "REFERENCE")) return True def v7_symbolic_identities(): """Symbolic verification of the three algebraic identities. Monte Carlo cannot check an identity, only an instance of it.""" try: import sympy as sp except ImportError: print("V7 symbolic: sympy not installed, SKIPPED") ROWS.append(("V7 symbolic identities", "-", "-", "-", "-", "SKIPPED")) return True ok = True # refresh entropy: L sign bits, an entry permutation, a user permutation L, Uu = 64, 4 ent = sp.Integer(L) + sp.log(sp.factorial(L), 2) + sp.log(sp.factorial(Uu), 2) ok &= abs(float(ent) - 364.6) < 0.05 print("V7a refresh entropy %.3f bits/block" % float(ent)) # fixed-key entropy: ordered choices of U of the L-1 non-constant rows fixed = sp.log(sp.factorial(L - 1) / sp.factorial(L - 1 - Uu), 2) ok &= abs(float(fixed) - 23.8) < 0.05 ok &= sp.factorial(L - 1) / sp.factorial(L - 1 - Uu) == 14295960 print("V7b fixed-key entropy %.3f bits, %d choices" % (float(fixed), sp.factorial(L - 1) / sp.factorial(L - 1 - Uu))) # the termwise Hadamard identity behind eq:hadamard n = 4 e = sp.Matrix(sp.symbols("e1:%d" % (n + 1))) m = sp.Matrix(sp.symbols("m1:%d" % (n + 1))) f = sp.Matrix(sp.symbols("f1:%d" % (n + 1))) lhs = sum((sp.matrix_multiply_elementwise(e, m))[k] * (sp.matrix_multiply_elementwise(f, m))[k] for k in range(n)) rhs = (e.T * sp.diag(*[m[k] ** 2 for k in range(n)]) * f)[0] ok &= sp.simplify(lhs - rhs) == 0 print("V7c (e*m).(f*m) == e^T diag(m^2) f :", sp.simplify(lhs - rhs) == 0) ROWS.append(("V7 symbolic identities", "exact", "exact", "0", "0", "PASS" if ok else "FAIL")) return ok def v8_cross_period_terms(): """Proposition 2 keeps only the diagonal of the jammer projection. The periodic key makes entries one period apart identical, so the cross-period terms do not vanish termwise. They are zero mean over the codebook, which is the claim the proof rests on.""" import math import torch m = cached_main_model() Bn = m.unit_codebook().detach().cpu() pat = m.masks().detach().cpu()[0] L, P, d = m.L, m.P, m.d # a generator of its own, seeded after the model is built: seeding the # global one first leaves the draw dependent on how main_model consumed # it, which moved this number between runs g = torch.Generator().manual_seed(7) rel = [] for _ in range(20000): w = torch.randn(d, generator=g) w /= w.norm() i = torch.randint(m.vu, (P,), generator=g) e = (Bn[i] / math.sqrt(P)).reshape(-1) a = (w * e).reshape(P, L) * pat[None, :] diag = float((a ** 2).sum()) rel.append((float((a.sum(0) ** 2).sum()) - diag) / diag) mean = sum(rel) / len(rel) ok = abs(mean) < 0.0005 print("V8 cross-period remainder, mean %+.4f of the retained term" % mean) ROWS.append(("V8 cross-period remainder", "0.0", "%.6f" % mean, "%.6f" % abs(mean), "0.0005", "PASS" if ok else "FAIL")) return ok def v9_score_variance_ratio(): """Per-digit score-variance ratio of Proposition 1's proof: the mean of sum_j e_j^4 / sum_j e_j^2 e'_j^2 over ordered codeword pairs of the trained unit codebook, quoted as 2.8 in the manuscript.""" import torch m = cached_main_model() B = m.unit_codebook().detach().cpu().double() B = B / B.norm(dim=1, keepdim=True) n = B.shape[0] num = (B ** 4).sum(1) ratios = [] for i in range(n): for j in range(n): if i != j: ratios.append(float(num[i] / (B[i] ** 2 * B[j] ** 2).sum())) mean = sum(ratios) / len(ratios) ok = abs(mean - 2.8) < 0.05 print("V9 per-digit score-variance ratio: %.3f (quoted 2.8)" % mean) ROWS.append(("V9 score-variance ratio", "2.8", "%.4f" % mean, "%.4f" % abs(mean - 2.8), "0.05", "PASS" if ok else "FAIL")) return ok def v10_format_matched_oma(): """The format-matched OMA reference of Section VI-B. The binary reference spends 16 of its 64 exclusive dimensions on antipodal bits. The same allocation spent the way the proposed scheme spends it, P=4 sixteen-ary orthogonal decisions over 16 dimensions each, is the comparison a reviewer will ask for.""" from sse_lib import oma_ser_orth from exp_full import oma_ser_keylen val = oma_ser_orth([10.0])[0] binary = oma_ser_keylen(64, 10.0) ok = abs(val - 0.055) < 0.001 print("V10 format-matched OMA at 10 dB: %.5f (binary %.5f)" % (val, binary)) ROWS.append(("V10 format-matched OMA at 10 dB", "0.055", "%.5f" % val, "%.5f" % abs(val - 0.055), "0.001", "PASS" if ok else "FAIL")) return ok def v11_oma_closed_form_vs_mc(): """The manuscript says the OMA closed form agrees with Monte Carlo to within one percent. That check had no stored artifact.""" from sse_lib import oma_ser, oma_ser_mc cf = oma_ser([16.0])[0] mc = oma_ser_mc([16.0], frames=2_000_000)[0] rel = abs(cf - mc) / mc ok = rel < 0.01 print("V11 OMA closed form %.6f vs Monte Carlo %.6f (%.2f%%)" % (cf, mc, 100 * rel)) ROWS.append(("V11 OMA closed form vs Monte Carlo", "%.6f" % mc, "%.6f" % cf, "%.4f" % rel, "0.01", "PASS" if ok else "FAIL")) return ok def main(): print(f"config d={D} U={U} V={V}\n") results = { "V1": v1_legit_self_alignment(), "V2": v2_eve_uninformed(), "V3": v3_leakage_vs_correlation(), "V4": v4_blind_jammer_spread(), "V5": v5_matched_jammer_concentrates(), "V6": v6_coded_oma_outage(), "V7": v7_symbolic_identities(), "V8": v8_cross_period_terms(), "V9": v9_score_variance_ratio(), "V10": v10_format_matched_oma(), "V11": v11_oma_closed_form_vs_mc(), } print("\nsummary:", {k: ("PASS" if v else "FAIL") for k, v in results.items()}) print("ALL PASS" if all(results.values()) else "SOME FAILED") # stored artifact so every quoted verification number has a raw file import csv as _csv from pathlib import Path as _Path data = _Path(__file__).resolve().parents[1] / "data" with open(data / "verify_math.csv", "w", newline="") as f: w = _csv.writer(f) w.writerow(["check", "claim", "empirical", "abs_err", "tol", "verdict"]) w.writerows(ROWS) print("[csv]", data / "verify_math.csv") if __name__ == "__main__": main()