Add the diagnostics behind the no-interference claim

code/diag_interference.py

@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+"""Where does the legitimate advantage over OMA go?
+
+An ideal M-ary receiver at the main configuration should reach 0.199 at
+10 dB against the 0.275 of resource-matched OMA, a factor of 1.38, while
+the system measures 0.257, a factor of 1.07. This script splits the
+shortfall into its two causes: residual multi-user interference, which
+orthogonal keys do not remove because masking is elementwise, and the
+distance the trained unit codebook falls short of an orthogonal set.
+"""
+import math
+import sys
+from pathlib import Path
+
+import torch
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+import sse_lib as L
+from sse_lib import DEVICE, snr_to_sigma2, rayleigh_gain
+from exp_full import main_model
+
+SNR_DB = 10.0
+FRAMES = 400_000
+CH = 40_000
+
+
+def ser(m, snr_db, frames, solo=False):
+    """SER of user 0. With solo=True the other users transmit nothing,
+    while the power normalizer c is left at its four-user value so that
+    user 0 keeps exactly the energy it has in the real system."""
+    tot = wrong = 0
+    with torch.no_grad():
+        while tot < frames:
+            n = min(CH, frames - tot)
+            dig = torch.randint(m.vu, (n, m.users, m.P), device=DEVICE)
+            Bn = m.unit_codebook()
+            e = Bn[dig] / math.sqrt(m.P)
+            mk = m.masks()
+            x = e * mk[None, :, None, :]
+            if solo:
+                x = x[:, :1]
+            y = x.sum(dim=1) / m.c
+            h = rayleigh_gain((n, 1), device=DEVICE)
+            sig = snr_to_sigma2(torch.full((n,), snr_db), m.d).to(DEVICE).sqrt()
+            rx = h[:, :, None, None] * y[:, None] \
+                + sig[:, None, None, None] * torch.randn(n, 1, m.P, m.L,
+                                                         device=DEVICE)
+            r = rx / h[:, :, None, None].clamp_min(1e-6)
+            cand = Bn[None, :, :] * mk[:1, None, :]
+            sc = torch.einsum("nupl,uvl->nupv", r, cand)
+            bad = (sc.argmax(-1)[:, 0] != dig[:, 0]).any(dim=-1)
+            wrong += int(bad.sum())
+            tot += n
+    return wrong / tot
+
+
+def main():
+    m = main_model()
+    with torch.no_grad():
+        Bn = m.unit_codebook()
+        G = Bn @ Bn.T
+        off = G - torch.diag(torch.diag(G))
+        mk = m.masks()
+        Gm = mk @ mk.T / m.L
+        offm = Gm - torch.diag(torch.diag(Gm))
+
+    print("main configuration: d=%d P=%d L=%d Vu=%d U=%d"
+          % (m.d, m.P, m.L, m.vu, m.users))
+    print("key cross-correlation, max |off-diagonal| : %.2e"
+          % offm.abs().max())
+    print("codebook Gram, max |off-diagonal|         : %.4f"
+          % off.abs().max())
+    print("codebook Gram, rms off-diagonal           : %.4f"
+          % off.pow(2).sum().div(m.vu * (m.vu - 1)).sqrt())
+    print("(an orthogonal set of %d codewords in %d dims would read 0)"
+          % (m.vu, m.L))
+    print()
+    four = ser(m, SNR_DB, FRAMES, solo=False)
+    solo = ser(m, SNR_DB, FRAMES, solo=True)
+    print("user-0 SER, all four users transmitting   : %.4f" % four)
+    print("user-0 SER, other users silent            : %.4f" % solo)
+    print("OMA, resource matched (closed form)       : %.4f"
+          % L.oma_ser([SNR_DB])[0])
+    print("ideal 16-ary orthogonal (separate MC)     : 0.1986")
+
+
+if __name__ == "__main__":
+    main()

code/diag_orthobook.py

@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+"""Does an orthogonal unit codebook recover the shortfall?
+
+diag_interference shows the gap to the ideal M-ary receiver is not
+multi-user interference but the geometry of the trained unit codebook,
+whose Gram matrix carries a root-mean-square off-diagonal of 0.45 where
+an orthogonal set would carry zero. Vu = L = 16 admits an exactly
+orthogonal set, so this measures what installing one buys.
+
+Two orthogonal sets are tried, because the choice is not free. The
+Walsh-Hadamard set collides with the keys: the rows are closed under the
+elementwise product, so masking a Hadamard codeword by a Hadamard key
+returns another Hadamard codeword and every user ends up with the same
+candidate set. A random orthogonal set carries no such group structure,
+and masking by a unit-modulus key preserves its orthogonality exactly.
+"""
+import sys
+from pathlib import Path
+
+import torch
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+import sse_lib as L
+from sse_lib import DEVICE, SSE
+from exp_full import hadamard, base_keys
+from diag_interference import ser
+
+SNR = [0.0, 10.0, 20.0]
+FRAMES = 400_000
+
+
+def fixed_model(B, P=4, vu=16, d=64, U=4):
+    L.set_seed(1)
+    m = SSE(P=P, vu=vu, d=d, users=U).to(DEVICE)
+    with torch.no_grad():
+        m.B.copy_(B.to(DEVICE))
+        m.W.copy_(base_keys(U, d // P).to(DEVICE))
+    m.calibrate_power()
+    return m
+
+
+def hadamard_book(vu=16, Lp=16):
+    B = torch.zeros(vu, Lp)
+    B[:, :vu] = torch.tensor(hadamard(vu).copy(), dtype=torch.float32)
+    return B
+
+
+def random_ortho_book(vu=16, Lp=16, seed=7):
+    g = torch.Generator().manual_seed(seed)
+    A = torch.randn(Lp, Lp, generator=g)
+    Q, _ = torch.linalg.qr(A)
+    return Q[:vu].contiguous()
+
+
+def report(name, m):
+    with torch.no_grad():
+        Bn = m.unit_codebook()
+        G = Bn @ Bn.T
+        off = (G - torch.diag(torch.diag(G))).abs().max()
+    row = [name, "%.2e" % off]
+    for s in SNR:
+        row.append("%.4f" % ser(m, s, FRAMES))
+    print("%-22s %-10s %-9s %-9s %-9s" % tuple(row))
+
+
+def main():
+    print("%-22s %-10s %-9s %-9s %-9s"
+          % ("unit codebook", "max|off|", "0 dB", "10 dB", "20 dB"))
+    report("Walsh-Hadamard", fixed_model(hadamard_book()))
+    report("random orthogonal", fixed_model(random_ortho_book()))
+    print("%-22s %-10s %-9s %-9s %-9s"
+          % ("trained (paper)", "0.887", "0.8822", "0.2576", "0.0307"))
+    print("%-22s %-10s %-9s %-9s %-9s"
+          % ("OMA, resource matched", "-",
+             "%.4f" % L.oma_ser([0.0])[0],
+             "%.4f" % L.oma_ser([10.0])[0],
+             "%.4f" % L.oma_ser([20.0])[0]))
+
+
+
+
+
+def solo_check():
+    """Splitting each candidate set from the superposition it must live
+    in. Orthogonal codewords are ideal for one user alone and are what
+    the single-user bound assumes, but the masked sets of different users
+    are then far from orthogonal to each other."""
+    print()
+    print("%-22s %-12s %-12s" % ("unit codebook", "solo 10 dB", "4-user 10 dB"))
+    for name, B in (("random orthogonal", random_ortho_book()),
+                    ("trained (retrain)", None)):
+        if B is None:
+            from exp_full import main_model
+            m = main_model()
+        else:
+            m = fixed_model(B)
+        print("%-22s %-12.4f %-12.4f"
+              % (name, ser(m, 10.0, FRAMES, solo=True),
+                 ser(m, 10.0, FRAMES, solo=False)))
+    print("single-user ideal M-ary bound (separate MC): 0.1986")
+
+
+if __name__ == "__main__":
+    main()
+    solo_check()