482 lines
21 KiB
Python
482 lines
21 KiB
Python
"""
|
|
GPU-accelerated Monte Carlo experiments (torch backend).
|
|
========================================================
|
|
Computes the CSV artifacts of experiments E2, E3, E4, E5, E7c, E8
|
|
of revision_sims.py with identical models and conventions, using
|
|
torch (CUDA when available) for the dense linear algebra. All
|
|
random draws come from the numpy generator with the documented seed
|
|
2026, so the sample stream is platform-independent; torch only
|
|
accelerates QR, matrix products, and linear solves in float32 /
|
|
complex64 precision. Figures are rendered separately by
|
|
replot_all.py, which reads only data/.
|
|
|
|
Run under WSL: python3 revision_sims_gpu.py E2 E3 E4 E5 E7c E8
|
|
"""
|
|
from __future__ import annotations
|
|
import csv
|
|
import math
|
|
import sys
|
|
import time
|
|
from pathlib import Path
|
|
import numpy as np
|
|
import torch
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
CSV_DIR = ROOT / "data"
|
|
|
|
SEED = 2026
|
|
rng = np.random.default_rng(SEED)
|
|
DEV = "cuda" if torch.cuda.is_available() else "cpu"
|
|
print(f"[gpu] device = {DEV}")
|
|
|
|
|
|
def write_csv(name, header, rows):
|
|
p = CSV_DIR / f"{name}.csv"
|
|
with open(p, "w", newline="") as f:
|
|
w = csv.writer(f); w.writerow(header); w.writerows(rows)
|
|
print(f"[OK] wrote {p}")
|
|
|
|
|
|
def haar_g(d):
|
|
"""Haar orthogonal on the GPU from a numpy Gaussian draw."""
|
|
G = torch.tensor(rng.standard_normal((d, d)), dtype=torch.float32,
|
|
device=DEV)
|
|
Q, R = torch.linalg.qr(G)
|
|
return Q * torch.sign(torch.diagonal(R)).unsqueeze(0)
|
|
|
|
|
|
def embed_pair(d, beta):
|
|
e1 = rng.standard_normal(d)
|
|
e1 /= np.linalg.norm(e1)
|
|
w = rng.standard_normal(d)
|
|
w -= (w @ e1) * e1
|
|
w /= np.linalg.norm(w)
|
|
e2 = beta * e1 + math.sqrt(1.0 - beta**2) * w
|
|
return (torch.tensor(e1, dtype=torch.float32, device=DEV),
|
|
torch.tensor(e2, dtype=torch.float32, device=DEV))
|
|
|
|
|
|
def rayleigh2():
|
|
h = (rng.standard_normal(2) + 1j * rng.standard_normal(2)) / math.sqrt(2)
|
|
return complex(h[0]), complex(h[1])
|
|
|
|
|
|
def cnoise_g(d):
|
|
n = (rng.standard_normal(d) + 1j * rng.standard_normal(d)) / math.sqrt(2)
|
|
return torch.tensor(n, dtype=torch.complex64, device=DEV)
|
|
|
|
|
|
def aware_g(t, Q, beta, c, nvar):
|
|
"""Batched aware Wiener demux. t: (b,d) cfloat, Q: (d,d) float,
|
|
c: python complex, nvar: (b,) tensor."""
|
|
d = Q.shape[0]
|
|
b = t.shape[0]
|
|
g = 1.0 - beta * beta
|
|
rho = g * abs(c)**2 / d + nvar
|
|
A = torch.eye(d, device=DEV, dtype=torch.complex64) \
|
|
+ beta * c * Q.to(torch.complex64)
|
|
S = (A @ A.mH / d).unsqueeze(0) \
|
|
+ rho.view(b, 1, 1) * torch.eye(d, device=DEV,
|
|
dtype=torch.complex64)
|
|
x = torch.linalg.solve(S, t.unsqueeze(-1))
|
|
return (A.mH.unsqueeze(0) @ x).squeeze(-1) / d
|
|
|
|
|
|
def abscos(a, e):
|
|
num = (a * e.to(a.dtype).conj()).sum(-1).abs()
|
|
return (num / (a.norm(dim=-1) * e.norm())).cpu().numpy()
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E2_sic(beta=0.311, d=512, ntr=400):
|
|
print("\n=== E2: receiver comparison under block-Rayleigh fading ===")
|
|
snr_db = np.arange(0, 31, 2.5)
|
|
sigs = torch.tensor(10 ** (-snr_db / 20.0), dtype=torch.float32,
|
|
device=DEV)
|
|
nb = len(snr_db)
|
|
res = {k: np.zeros(nb) for k in ("edma", "blind", "oma", "genie", "sic")}
|
|
t0 = time.time()
|
|
for tr in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
h1, h2 = rayleigh2()
|
|
n = cnoise_g(d); n2 = cnoise_g(d)
|
|
r0 = h1 * (M1 @ e1).to(torch.complex64) \
|
|
+ h2 * (M2 @ e2).to(torch.complex64)
|
|
r = r0.unsqueeze(0) + sigs.view(-1, 1) * n.unsqueeze(0)
|
|
t1 = (M1.T.to(torch.complex64) @ r.unsqueeze(-1)).squeeze(-1) / h1
|
|
t2 = (M2.T.to(torch.complex64) @ r.unsqueeze(-1)).squeeze(-1) / h2
|
|
c1, c2 = h2 / h1, h1 / h2
|
|
v1 = sigs**2 / abs(h1)**2
|
|
v2 = sigs**2 / abs(h2)**2
|
|
g1 = aware_g(t1, Q, beta, c1, v1)
|
|
g2 = aware_g(t2, Q.T, beta, c2, v2)
|
|
res["edma"] += 0.5 * (abscos(g1, e1) + abscos(g2, e2))
|
|
res["blind"] += 0.5 * (abscos(t1, e1) + abscos(t2, e2))
|
|
o1 = e1.to(torch.complex64).unsqueeze(0) \
|
|
+ math.sqrt(2) * sigs.view(-1, 1) * n.unsqueeze(0) / h1
|
|
o2 = e2.to(torch.complex64).unsqueeze(0) \
|
|
+ math.sqrt(2) * sigs.view(-1, 1) * n2.unsqueeze(0) / h2
|
|
res["oma"] += 0.5 * (abscos(o1, e1) + abscos(o2, e2))
|
|
ge1 = (M1.T.to(torch.complex64)
|
|
@ (r - h2 * (M2 @ e2).to(torch.complex64)).unsqueeze(-1)
|
|
).squeeze(-1) / h1
|
|
ge2 = (M2.T.to(torch.complex64)
|
|
@ (r - h1 * (M1 @ e1).to(torch.complex64)).unsqueeze(-1)
|
|
).squeeze(-1) / h2
|
|
res["genie"] += 0.5 * (abscos(ge1, e1) + abscos(ge2, e2))
|
|
# realizable decision-directed SIC: the stronger user is detected
|
|
# with the same aware Wiener stage (a scalar-scaled matched-filter
|
|
# decision would re-modulate to a multiple of r itself, because
|
|
# M_s M_s^T = I, and cancel nothing), its re-modulated estimate is
|
|
# subtracted, and the weaker user is read from the residual.
|
|
if abs(h1) >= abs(h2):
|
|
hs, hw, Ms, Mw, es, ew = h1, h2, M1, M2, e1, e2
|
|
Qsw, csw, vsw = Q, c1, v1
|
|
else:
|
|
hs, hw, Ms, Mw, es, ew = h2, h1, M2, M1, e2, e1
|
|
Qsw, csw, vsw = Q.T, c2, v2
|
|
t_s = (Ms.T.to(torch.complex64) @ r.unsqueeze(-1)).squeeze(-1) / hs
|
|
dec = aware_g(t_s, Qsw, beta, csw, vsw)
|
|
r_res = r - hs * (Ms.to(torch.complex64)
|
|
@ dec.unsqueeze(-1)).squeeze(-1)
|
|
d_w = (Mw.T.to(torch.complex64) @ r_res.unsqueeze(-1)).squeeze(-1) / hw
|
|
res["sic"] += 0.5 * (abscos(dec, es) + abscos(d_w, ew))
|
|
if (tr + 1) % 100 == 0:
|
|
print(f" {tr+1}/{ntr} ({time.time()-t0:.0f}s)", flush=True)
|
|
for k in res:
|
|
res[k] /= ntr
|
|
rows = [[s] + [res[k][i] for k in ("edma", "blind", "oma", "genie", "sic")]
|
|
for i, s in enumerate(snr_db)]
|
|
write_csv("sic_comparison",
|
|
["snr_db", "edma", "blind", "oma", "genie", "sic"], rows)
|
|
i20 = list(snr_db).index(20)
|
|
print(f" at 20 dB: EDMA {res['edma'][i20]:.3f}, blind "
|
|
f"{res['blind'][i20]:.3f}, SIC {res['sic'][i20]:.3f}, "
|
|
f"genie {res['genie'][i20]:.3f}, OMA {res['oma'][i20]:.3f}")
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E3_unconditional(beta=0.311, d=512, ntr=1500):
|
|
print("\n=== E3: Rayleigh unconditional MSE of the aware receiver ===")
|
|
rows = []
|
|
for s in (10, 20):
|
|
sig = 10 ** (-s / 20.0)
|
|
mses, bl = [], []
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
h1, h2 = rayleigh2()
|
|
n = cnoise_g(d)
|
|
r = h1 * (M1 @ e1).to(torch.complex64) \
|
|
+ h2 * (M2 @ e2).to(torch.complex64) + sig * n
|
|
t1 = (M1.T.to(torch.complex64) @ r) / h1
|
|
c1 = h2 / h1
|
|
nv = torch.tensor([sig**2 / abs(h1)**2], device=DEV)
|
|
g1 = aware_g(t1.unsqueeze(0), Q, beta, c1, nv)[0]
|
|
mses.append(float((g1 - e1.to(torch.complex64)).norm()**2))
|
|
lam = (1.0 / d) / (1.0 / d + abs(c1)**2 / d + sig**2 / abs(h1)**2)
|
|
b1 = lam * t1
|
|
bl.append(float((b1 - e1.to(torch.complex64)).norm()**2))
|
|
rows.append([s, float(np.mean(mses)), float(np.median(mses)),
|
|
float(np.mean(bl)), float(np.median(bl))])
|
|
print(f" {s} dB: aware mean {rows[-1][1]:.4f} "
|
|
f"(median {rows[-1][2]:.4f}) | blind mean {rows[-1][3]:.4f} "
|
|
f"(median {rows[-1][4]:.4f})")
|
|
write_csv("rayleigh_mse", ["snr_db", "aware_mean", "aware_median",
|
|
"blind_mean", "blind_median"], rows)
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E4_csi(beta=0.311, d=512, snr=30.0, ntr=400):
|
|
print("\n=== E4: imperfect CSI robustness (EDMA vs realizable SIC) ===")
|
|
sh2 = np.array([0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.3])
|
|
sig = 10 ** (-snr / 20.0)
|
|
res_e = np.zeros(len(sh2)); res_s = np.zeros(len(sh2))
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
h1, h2 = rayleigh2()
|
|
n = cnoise_g(d)
|
|
r = h1 * (M1 @ e1).to(torch.complex64) \
|
|
+ h2 * (M2 @ e2).to(torch.complex64) + sig * n
|
|
eps1, eps2 = rayleigh2()
|
|
for j, v in enumerate(sh2):
|
|
hh1 = h1 + math.sqrt(v) * eps1
|
|
hh2 = h2 + math.sqrt(v) * eps2
|
|
t1 = (M1.T.to(torch.complex64) @ r) / hh1
|
|
t2 = (M2.T.to(torch.complex64) @ r) / hh2
|
|
nv1 = torch.tensor([sig**2 / abs(hh1)**2], device=DEV)
|
|
nv2 = torch.tensor([sig**2 / abs(hh2)**2], device=DEV)
|
|
g1 = aware_g(t1.unsqueeze(0), Q, beta, hh2 / hh1, nv1)[0]
|
|
g2 = aware_g(t2.unsqueeze(0), Q.T, beta, hh1 / hh2, nv2)[0]
|
|
res_e[j] += 0.5 * (float(abscos(g1.unsqueeze(0), e1)[0])
|
|
+ float(abscos(g2.unsqueeze(0), e2)[0]))
|
|
if abs(hh1) >= abs(hh2):
|
|
hs, hw, Ms, Mw, es, ew = hh1, hh2, M1, M2, e1, e2
|
|
Qsw, csw = Q, hh2 / hh1
|
|
else:
|
|
hs, hw, Ms, Mw, es, ew = hh2, hh1, M2, M1, e2, e1
|
|
Qsw, csw = Q.T, hh1 / hh2
|
|
t_s = (Ms.T.to(torch.complex64) @ r) / hs
|
|
nvs = torch.tensor([sig**2 / abs(hs)**2], device=DEV)
|
|
dec = aware_g(t_s.unsqueeze(0), Qsw, beta, csw, nvs)[0]
|
|
r_res = r - hs * (Ms.to(torch.complex64) @ dec)
|
|
d_w = (Mw.T.to(torch.complex64) @ r_res) / hw
|
|
res_s[j] += 0.5 * (float(abscos(dec.unsqueeze(0), es)[0])
|
|
+ float(abscos(d_w.unsqueeze(0), ew)[0]))
|
|
res_e /= ntr; res_s /= ntr
|
|
print(f" EDMA: {res_e[0]:.4f} -> {res_e[-1]:.4f} "
|
|
f"(delta {100*(res_e[0]-res_e[-1]):.2f} points)")
|
|
print(f" SIC : {res_s[0]:.4f} -> {res_s[-1]:.4f} "
|
|
f"(delta {100*(res_s[0]-res_s[-1]):.2f} points)")
|
|
rows = [[v, res_e[j], res_s[j]] for j, v in enumerate(sh2)]
|
|
write_csv("csi_error", ["sigma_h2", "edma", "sic"], rows)
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E5_maskfam(beta=0.311, d=512, ntr=200):
|
|
print("\n=== E5: Walsh-Hadamard diagonal variant vs Haar ===")
|
|
snr_db = np.arange(0, 41, 5)
|
|
sigs = torch.tensor(10 ** (-snr_db / 20.0), dtype=torch.float32,
|
|
device=DEV)
|
|
nb = len(snr_db)
|
|
H = np.array([[1.0]])
|
|
while H.shape[0] < d:
|
|
H = np.block([[H, H], [H, -H]])
|
|
Ht = torch.tensor(H / math.sqrt(d), dtype=torch.float32, device=DEV)
|
|
g = 1.0 - beta**2
|
|
res = {"haar": np.zeros(nb), "wh": np.zeros(nb)}
|
|
exact = np.zeros(nb)
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
D1 = torch.tensor(np.sign(rng.standard_normal(d)),
|
|
dtype=torch.float32, device=DEV)
|
|
D2 = torch.tensor(np.sign(rng.standard_normal(d)),
|
|
dtype=torch.float32, device=DEV)
|
|
W1, W2 = Ht * D1.unsqueeze(0), Ht * D2.unsqueeze(0)
|
|
q = D1 * D2
|
|
n = cnoise_g(d)
|
|
r0h = (M1 @ e1 + M2 @ e2).to(torch.complex64)
|
|
r0w = (W1 @ e1 + W2 @ e2).to(torch.complex64)
|
|
rh = r0h.unsqueeze(0) + sigs.view(-1, 1) * n.unsqueeze(0)
|
|
rw = r0w.unsqueeze(0) + sigs.view(-1, 1) * n.unsqueeze(0)
|
|
t1 = (M1.T.to(torch.complex64) @ rh.unsqueeze(-1)).squeeze(-1)
|
|
g1 = aware_g(t1, Q, beta, 1.0, sigs**2)
|
|
res["haar"] += abscos(g1, e1)
|
|
tw = (W1.T.to(torch.complex64) @ rw.unsqueeze(-1)).squeeze(-1)
|
|
a = 1.0 + beta * q # (d,)
|
|
rho = g / d + sigs**2 # (nb,)
|
|
wdiag = a.unsqueeze(0) / (a.unsqueeze(0)**2 / d
|
|
+ rho.view(-1, 1)) # (nb,d)
|
|
w1 = (wdiag.to(torch.complex64) / d) * tw
|
|
res["wh"] += abscos(w1, e1)
|
|
exact += ((g + d * sigs.view(-1, 1)**2)
|
|
/ (a.unsqueeze(0)**2 + g + d * sigs.view(-1, 1)**2)
|
|
).mean(1).cpu().numpy()
|
|
for k in res:
|
|
res[k] /= ntr
|
|
exact /= ntr
|
|
print(f" max |WH - Haar| cosine dev: "
|
|
f"{100*np.max(np.abs(res['wh']-res['haar'])):.2f} points; "
|
|
f"40 dB WH {res['wh'][-1]:.4f} vs Haar {res['haar'][-1]:.4f}")
|
|
rows = [[s, res["haar"][i], res["wh"][i], exact[i]]
|
|
for i, s in enumerate(snr_db)]
|
|
write_csv("mask_family_rev", ["snr_db", "haar", "wh", "wh_exact_mse"],
|
|
rows)
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E7_multiuser(beta=0.311, d=512, ntr=100):
|
|
print("\n=== E7c: multi-user scaling (joint Wiener) ===")
|
|
snr_db = np.arange(0, 31, 2.5)
|
|
sigs = torch.tensor(10 ** (-snr_db / 20.0), dtype=torch.float32,
|
|
device=DEV)
|
|
nb = len(snr_db)
|
|
rows = []
|
|
for U in (2, 3, 4):
|
|
B = (1 - beta) * np.eye(U) + beta * np.ones((U, U))
|
|
A = np.linalg.cholesky(B)
|
|
Bt = torch.tensor(B, dtype=torch.float32, device=DEV)
|
|
err = np.zeros(nb)
|
|
t0 = time.time()
|
|
for _ in range(ntr):
|
|
F = rng.standard_normal((d, U))
|
|
Fq, _ = np.linalg.qr(F)
|
|
E = (Fq @ A.T).T
|
|
Et = torch.tensor(E, dtype=torch.float32, device=DEV)
|
|
masks = [haar_g(d) for _ in range(U)]
|
|
n = cnoise_g(d)
|
|
r0 = sum(masks[u] @ Et[u] for u in range(U)).to(torch.complex64)
|
|
Qs = [masks[0].T @ masks[v] for v in range(U)]
|
|
Ret = sum(Bt[0, v] * Qs[v].T for v in range(U)) / d
|
|
S0 = sum(Bt[v, w] * (Qs[v] @ Qs[w].T)
|
|
for v in range(U) for w in range(U)) / d
|
|
r = r0.unsqueeze(0) + sigs.view(-1, 1) * n.unsqueeze(0)
|
|
t = (masks[0].T.to(torch.complex64)
|
|
@ r.unsqueeze(-1)).squeeze(-1)
|
|
S = S0.to(torch.complex64).unsqueeze(0) \
|
|
+ (sigs**2).view(-1, 1, 1) \
|
|
* torch.eye(d, device=DEV, dtype=torch.complex64)
|
|
x = torch.linalg.solve(S, t.unsqueeze(-1))
|
|
eh = (Ret.to(torch.complex64).unsqueeze(0) @ x).squeeze(-1)
|
|
err += ((eh - Et[0].to(torch.complex64)).norm(dim=1)**2
|
|
).cpu().numpy()
|
|
err /= ntr
|
|
T_mc = U * np.log2(1.0 / err)
|
|
r0v = (U - 1) + d / 10 ** (snr_db / 10.0)
|
|
T_bl = U * np.log2(1.0 + 1.0 / r0v)
|
|
T_oma = U * np.log2(1 + 10 ** (snr_db / 10.0) / (U * d))
|
|
for i, s in enumerate(snr_db):
|
|
rows.append([U, s, T_mc[i], T_bl[i], T_oma[i], err[i]])
|
|
i20 = list(snr_db).index(20.0)
|
|
print(f" U={U}: at 20 dB EDMA {T_mc[i20]:.3f} vs blind "
|
|
f"{T_bl[i20]:.3f} vs OMA {T_oma[i20]:.3f} "
|
|
f"(gain {T_mc[i20]/T_oma[i20]:.2f}x), floor MSE {err[-1]:.4f}"
|
|
f" [{time.time()-t0:.0f}s]")
|
|
write_csv("multiuser_corrected",
|
|
["U", "snr_db", "edma_mc", "blind", "oma", "mse_mc"], rows)
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E8_mismatch(beta=0.311, d=512, snr=20.0, ntr=200):
|
|
print("\n=== E8: affinity mismatch of the aware receiver ===")
|
|
sig = 10 ** (-snr / 20.0)
|
|
bhs = [b for b in (beta - 0.1, beta - 0.06, beta, beta + 0.06,
|
|
beta + 0.1, beta + 2 ** -8, 0.0) if b >= 0]
|
|
accs = np.zeros(len(bhs)); msea = np.zeros(len(bhs))
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
n = cnoise_g(d)
|
|
r = (M1 @ e1 + M2 @ e2).to(torch.complex64) + sig * n
|
|
t1 = (M1.T.to(torch.complex64) @ r)
|
|
nv = torch.tensor([sig**2], device=DEV)
|
|
for j, bh in enumerate(bhs):
|
|
g1 = aware_g(t1.unsqueeze(0), Q, bh, 1.0, nv)[0]
|
|
accs[j] += float(abscos(g1.unsqueeze(0), e1)[0])
|
|
msea[j] += float((g1 - e1.to(torch.complex64)).norm()**2)
|
|
accs /= ntr; msea /= ntr
|
|
rows = []
|
|
for j, bh in enumerate(bhs):
|
|
tag = ("quant b=7" if abs(bh - beta - 2**-8) < 1e-12 else
|
|
("blind" if bh == 0.0 else f"delta={bh-beta:+.2f}"))
|
|
print(f" beta_hat={bh:.4f} ({tag}): cosine {accs[j]:.4f}, "
|
|
f"MSE {msea[j]:.4f}")
|
|
rows.append([bh, accs[j], msea[j]])
|
|
write_csv("mismatch", ["beta_hat", "cosine", "mse"], rows)
|
|
|
|
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E9_ceiling(d=512, snr=60.0, ntr=200):
|
|
print(chr(10) + '=== E9: cosine-ceiling verification (h=1) ===')
|
|
sig = 10 ** (-snr / 20.0)
|
|
rows = []
|
|
for beta in (0.311, 0.8):
|
|
g = 1.0 - beta**2
|
|
acc_a = 0.0; acc_b = 0.0
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
n = cnoise_g(d)
|
|
r = (M1 @ e1 + M2 @ e2).to(torch.complex64) + sig * n
|
|
t1 = (M1.T.to(torch.complex64) @ r)
|
|
nv = torch.tensor([sig**2], device=DEV)
|
|
g1 = aware_g(t1.unsqueeze(0), Q, beta, 1.0, nv)[0]
|
|
acc_a += float(abscos(g1.unsqueeze(0), e1)[0])
|
|
acc_b += float(abscos(t1.unsqueeze(0), e1)[0])
|
|
acc_a /= ntr; acc_b /= ntr
|
|
import math as _m
|
|
pred_a = _m.sqrt(1.0 - _m.sqrt(g) / 2.0)
|
|
pred_b = _m.sqrt(0.5)
|
|
print(f' beta={beta}: aware MC {acc_a:.4f} pred {pred_a:.4f} | '
|
|
f'blind MC {acc_b:.4f} pred {pred_b:.4f}')
|
|
rows.append([beta, acc_a, pred_a, acc_b, pred_b])
|
|
write_csv('cosine_ceiling', ['beta', 'aware_mc', 'aware_pred',
|
|
'blind_mc', 'blind_pred'], rows)
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
def E10_whpad(beta=0.311, d=768, dpad=1024, ntr=200):
|
|
"""Zero-padded WH at d=768 (padded to 1024) vs dense Haar at 768.
|
|
|
|
The embedding (768) is zero-padded to 1024, masked by H_1024 D_u,
|
|
and the exact per-coordinate Wiener uses the true prior (signal
|
|
variance 1/768 on the active support, zero on the padding, so the
|
|
padded coordinates are discarded). Reference: Haar masks at the
|
|
native d=768 with the standard aware demultiplexer. Same per-block
|
|
energy E_b = 1 and the same noise PSD; the padded block occupies
|
|
dpad channel uses, a bandwidth cost of dpad/d."""
|
|
print(chr(10) + '=== E10: zero-padded WH (768->1024) vs native Haar 768 ===')
|
|
snr_db = np.arange(0, 41, 5)
|
|
sigs = torch.tensor(10 ** (-snr_db / 20.0), dtype=torch.float32,
|
|
device=DEV)
|
|
nb = len(snr_db)
|
|
H = np.array([[1.0]])
|
|
while H.shape[0] < dpad:
|
|
H = np.block([[H, H], [H, -H]])
|
|
Ht = torch.tensor(H / math.sqrt(dpad), dtype=torch.float32, device=DEV)
|
|
g = 1.0 - beta**2
|
|
res = {'haar': np.zeros(nb), 'whpad': np.zeros(nb)}
|
|
for _ in range(ntr):
|
|
e1, e2 = embed_pair(d, beta)
|
|
# ---- native Haar at 768 ----
|
|
M1, M2 = haar_g(d), haar_g(d)
|
|
Q = M1.T @ M2
|
|
n = cnoise_g(d)
|
|
r = (M1 @ e1 + M2 @ e2).to(torch.complex64).unsqueeze(0) \
|
|
+ sigs.view(-1, 1) * n.unsqueeze(0)
|
|
t1 = (M1.T.to(torch.complex64) @ r.unsqueeze(-1)).squeeze(-1)
|
|
g1 = aware_g(t1, Q, beta, 1.0, sigs**2)
|
|
res['haar'] += abscos(g1, e1)
|
|
# ---- zero-padded WH at 1024 ----
|
|
z = torch.zeros(dpad - d, device=DEV)
|
|
e1p = torch.cat([e1, z]); e2p = torch.cat([e2, z])
|
|
D1 = torch.tensor(np.sign(rng.standard_normal(dpad)),
|
|
dtype=torch.float32, device=DEV)
|
|
D2 = torch.tensor(np.sign(rng.standard_normal(dpad)),
|
|
dtype=torch.float32, device=DEV)
|
|
W1, W2 = Ht * D1.unsqueeze(0), Ht * D2.unsqueeze(0)
|
|
npad = cnoise_g(dpad)
|
|
rp = (W1 @ e1p + W2 @ e2p).to(torch.complex64).unsqueeze(0) \
|
|
+ sigs.view(-1, 1) * npad.unsqueeze(0)
|
|
tw = (W1.T.to(torch.complex64) @ rp.unsqueeze(-1)).squeeze(-1)
|
|
q = (D1 * D2)[:d] # active coordinates only
|
|
a = 1.0 + beta * q # (d,)
|
|
s_var = 1.0 / d # true signal variance
|
|
v = g / d + (sigs**2).view(-1, 1) # interference + noise
|
|
gains = (s_var * a.unsqueeze(0)) / (a.unsqueeze(0)**2 * s_var + v)
|
|
w1 = gains.to(torch.complex64) * tw[:, :d]
|
|
res['whpad'] += abscos(w1, e1)
|
|
for k in res:
|
|
res[k] /= ntr
|
|
rows = [[s, res['haar'][i], res['whpad'][i]]
|
|
for i, s in enumerate(snr_db)]
|
|
write_csv('wh_padding', ['snr_db', 'haar768', 'whpad1024'], rows)
|
|
dev = res['haar'] - res['whpad']
|
|
print(f' cosine delta (haar - whpad): max {dev.max():.4f}, '
|
|
f'at 20 dB {dev[list(snr_db).index(20)]:.4f}, '
|
|
f'at 40 dB {dev[-1]:.4f}')
|
|
print(f' bandwidth cost: {dpad}/{d} = {dpad/d:.3f}x uses '
|
|
f'(per-use rate factor {d/dpad:.3f})')
|
|
|
|
|
|
if __name__ == "__main__":
|
|
todo = set(sys.argv[1:])
|
|
ALL = {"E2": E2_sic, "E3": E3_unconditional, "E4": E4_csi,
|
|
"E5": E5_maskfam, "E7c": E7_multiuser, "E8": E8_mismatch,
|
|
"E9": E9_ceiling, "E10": E10_whpad}
|
|
for name, fn in ALL.items():
|
|
if not todo or name in todo:
|
|
fn()
|
|
print("\nAll requested GPU simulations complete.")
|