"""Reproducible synthetic BPSK study. No real communication measurements."""
from pathlib import Path
import math,csv,json
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
OUT=Path(__file__).resolve().parent
SEED=240927; N=1_000_000; DBS=[-6,-3,0,3,6]
def q(x): return .5*math.erfc(x/math.sqrt(2))
def wilson(k,n):
    z=1.959963984540054;p=k/n;d=1+z*z/n;c=(p+z*z/(2*n))/d;h=z*math.sqrt(p*(1-p)/n+z*z/(4*n*n))/d
    return max(0.,c-h),min(1.,c+h)
def integral_q(x):
    # Independent numerical quadrature of the standard-normal density.
    n=20000; h=(12-x)/n
    f=lambda t:math.exp(-t*t/2)/math.sqrt(2*math.pi)
    return h/3*(f(x)+f(12)+4*sum(f(x+h*j) for j in range(1,n,2))+2*sum(f(x+h*j) for j in range(2,n,2)))
def write(name,rows):
    with (OUT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
rng=np.random.default_rng(SEED);rows=[];controls=[];checks=[]
for db in DBS:
    gamma=10**(db/10); sigma=math.sqrt(1/(2*gamma))
    bits=rng.integers(0,2,N);s=2.*bits-1
    noise=rng.normal(0,sigma,(3,N)); single=s+noise[0]
    summed=np.cumsum(noise,axis=0)
    for r in [1,2,3]:
        errors=int(np.count_nonzero((r*s+summed[r-1]>=0)!=bits))
        lo,hi=wilson(errors,N);theory=q(math.sqrt(2*r*gamma))
        rows.append(dict(per_use_energy_over_N0_dB=db,repetitions=r,bits=N,errors=errors,observed_ber=errors/N,wilson95_lower=lo,wilson95_upper=hi,theory_ber=theory,total_energy_multiplier=r,channel_uses_per_bit=r))
        z=(errors-N*theory)/math.sqrt(N*theory*(1-theory))
        checks.append(dict(db=db,repetitions=r,binomial_standardized_difference=z))
        assert abs(z)<6,'Unexpected Monte Carlo deviation; investigate rather than tune seed.'
        assert np.array_equal(r*single>=0,single>=0)
        controls.append(dict(per_use_energy_over_N0_dB=db,repetitions=r,replay_errors=int(np.count_nonzero((r*single>=0)!=bits)),identical_decisions_to_single=True))
write('repetition-monte-carlo.csv',rows);write('same-sample-control.csv',controls)
curves=[]
for db in np.linspace(-6,6,121):
    gamma=10**(db/10)
    for r in [1,2,3]:
        curves.append(dict(energy_over_N0_dB=float(db),repetitions=r,fixed_per_use_theory_ber=q(math.sqrt(2*r*gamma)),fixed_total_theory_ber=q(math.sqrt(2*gamma)),fixed_total_per_use_energy_fraction=1/r,channel_uses_per_bit=r))
write('repetition-analytic.csv',curves)
quad=max(abs(q(v)-integral_q(v)) for v in [0,.5,1,2,3,4,5])
assert quad<1e-12
# At fixed total energy E, sum mean is sqrt(r*E), noise variance r*N0/2.
# The standardized distance to zero therefore remains sqrt(2*E/N0).
assert all(abs(math.sqrt(r)*math.sqrt(2/r)-math.sqrt(2))<1e-14 for r in [1,2,3])
meta=dict(title='Repetition, noise and the energy budget',kind='Synthetic ideal coherent BPSK in independent additive white Gaussian noise',seed=SEED,rng='numpy default_rng PCG64',numpy_version=np.__version__,bits_per_condition=N,snr_grid_dB=DBS,formulas={'fixed_per_use':'BER=Q(sqrt(2*r*E_use/N0))','fixed_total':'BER=Q(sqrt(2*E_total/N0)); E_use=E_total/r','noise':'Per-use Gaussian noise variance=N0/2; unit per-use amplitude for Monte Carlo','confidence':'Pointwise 95% Wilson binomial interval; not a simultaneous confidence band'},assumptions=['Known symbol timing, known phase, ideal coherent detection and equiprobable independent bits.','Each retransmission has independent noise and the receiver adds analog-valued samples (soft combining).','At fixed per-use energy, r repeats cost r times the energy and r channel uses per information bit.','At fixed total energy, divide per-use energy by r; analytical BER is unchanged in this ideal channel.','Replaying the same stored noisy sample multiplies both signal and noise and cannot change its sign decision.'],limitations=['Not measured performance, speech intelligibility, room acoustics or a human listener.','No fading, interference, timing errors, phase errors, power limits, overhead or modern error-correcting code comparison.','Paired repetition conditions share samples; confidence intervals describe each BER, not independent between-condition differences.','Zero observed errors do not imply zero probability: retain the Wilson upper limit.'],validation={'max_Q_quadrature_abs_error':quad,'fixed_total_normalized_margin_invariant':True,'same_sample_decisions_identical':True,'monte_carlo_checks':checks})
(OUT/'study.json').write_text(json.dumps(meta,indent=2)+'\n',encoding='utf-8')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False})
fig,ax=plt.subplots(1,2,figsize=(11.5,4.8),constrained_layout=True)
colors=['#393f94','#087f8c','#b42e8d']
for r,c in zip([1,2,3],colors):
    rr=[x for x in curves if x['repetitions']==r]
    ax[0].semilogy([x['energy_over_N0_dB'] for x in rr],[x['fixed_per_use_theory_ber'] for x in rr],color=c,label=f'{r} use'+('s' if r>1 else '')+f'; {r}× total energy')
    for row in [x for x in rows if x['repetitions']==r]:
        y=row['observed_ber'];lo=row['wilson95_lower'];hi=row['wilson95_upper']
        if y>0:ax[0].errorbar(row['per_use_energy_over_N0_dB'],y,yerr=[[y-lo],[hi-y]],fmt='o',ms=4,color=c,capsize=3)
        else:ax[0].plot(row['per_use_energy_over_N0_dB'],hi,'v',color=c,ms=7)
xx=np.linspace(-6,6,121);yy=[q(math.sqrt(2*10**(x/10))) for x in xx]
ax[1].semilogy(xx,yy,color=colors[0],lw=3,label='1, 2 and 3 uses: identical ideal BER')
ax[1].text(.05,.12,'Repetition spreads the same energy\nover more channel uses.\nNo ideal AWGN energy-efficiency gain.',transform=ax[1].transAxes,bbox=dict(facecolor='white',edgecolor='#dddddd',pad=8),fontsize=10)
for a in ax:a.set_ylabel('Bit error probability');a.grid(True,which='both',alpha=.16);a.legend(fontsize=8,loc='upper right');a.set_xlim(-6.2,6.2)
ax[0].set_xlabel('Energy per channel use / N₀ (dB)');ax[0].set_title('A  |  More energy buys fewer errors',loc='left',weight='bold')
ax[1].set_xlabel('Total energy per information bit / N₀ (dB)');ax[1].set_title('B  |  Hold total energy constant',loc='left',weight='bold')
fig.suptitle('Repeating a bit: reliability has an energy and time cost',fontsize=15,weight='bold')
fig.supxlabel('Lines: analytical model. Dots: seeded simulation, 95% Wilson intervals. Triangle: upper limit when zero errors occur.',fontsize=9)
for ext in ['png','svg','pdf']:fig.savefig(OUT/f'repetition-energy.{ext}',dpi=180)
print(json.dumps({'rows':len(rows),'analytic_rows':len(curves),'max_Q_error':quad,'max_abs_mc_z':max(abs(x['binomial_standardized_difference']) for x in checks),'zero_count_cases':sum(x['errors']==0 for x in rows)}))