"""Independent local likelihood check; does not import or execute supplied research."""
from pathlib import Path
import numpy as np
import csv,json,time
OUT=Path(__file__).resolve().parent
seed=16020261101; rng=np.random.default_rng(seed); trials=2000;n=512;sigma=.1;true=.137
started=time.monotonic();rows=[];errors={}
x=(np.arange(n)-(n-1)/2)/n
A=np.exp(-2*(x/2.14)**2)
noise=sigma/np.sqrt(2)*(rng.standard_normal((trials,n))+1j*rng.standard_normal((trials,n)))
noise_stats={'expected_complex_power':sigma**2,'measured_complex_power':float(np.mean(abs(noise)**2)),'real_variance':float(np.var(noise.real)),'imag_variance':float(np.var(noise.imag))}
for offset in [0.,1.]:
 f=x+offset; w=2*np.pi*f;mu=A*np.exp(-1j*w*true);y=mu[None,:]+noise
 for phase in ['known','nuisance']:
  est=np.full(trials,true+.01)
  for iteration in range(12):
   if phase=='known':
    model=A[None,:]*np.exp(-1j*est[:,None]*w[None,:]);d1=-1j*w[None,:]*model;d2=-w[None,:]**2*model;res=y-model
    score=np.real(np.sum(np.conj(d1)*res,axis=1));curv=np.real(np.sum(np.conj(d2)*res-np.conj(d1)*d1,axis=1))
   else:
    z=A[None,:]*y*np.exp(1j*est[:,None]*w[None,:]);c=z.sum(axis=1);cp=(1j*w[None,:]*z).sum(axis=1);cpp=(-w[None,:]**2*z).sum(axis=1)
    score=2*np.real(cp*np.conj(c));curv=2*np.real(cpp*np.conj(c)+cp*np.conj(cp))
   step=score/curv;est-=step
  err=est-true; weights=A*A;fbar=np.sum(weights*f)/sum(weights)
  information=2/sigma**2*np.sum(weights*(2*np.pi*(f-(fbar if phase=='nuisance' else 0)))**2)
  bound=information**-.5;old=(np.sum(weights*w*w)/sigma**2)**-.5
  row={'offset_bandwidths':offset,'phase':phase,'trials':trials,'bias':float(err.mean()),'sample_std':float(err.std(ddof=1)),'rmse':float(np.sqrt(np.mean(err**2))),'corrected_crb_std':float(bound),'rmse_over_corrected':float(np.sqrt(np.mean(err**2))/bound),'rmse_over_source_formula':float(np.sqrt(np.mean(err**2))/old),'max_final_newton_step':float(np.max(abs(step))),'failures':int(np.sum(~np.isfinite(est)))}
  rows.append(row);errors[str(offset)+'-'+phase]=err
with (OUT/'monte_carlo.csv').open('w') as f:
 writer=csv.DictWriter(f,fieldnames=rows[0]);writer.writeheader();writer.writerows(rows)
np.savez_compressed(OUT/'trial_errors.npz',**errors)
c=299792458.;N=2048;lam=2e-6*np.sin(np.deg2rad(20));be=c/(N*lam);table=[]
for B in [10e9,50e9,100e9,214e9,500e9,1000e9,2000e9]:
 freq=(np.arange(4096)-2048)*B/4096
 for target in ['point','flat_single_delay']:
  amp=np.exp(-2*(freq/be)**2) if target=='point' else np.ones(4096);weights=amp**2;fbar=np.sum(weights*freq)/sum(weights)
  old=(np.sum((2*np.pi*freq)**2*weights)/sigma**2)**-.5
  known=(2*np.sum((2*np.pi*freq)**2*weights)/sigma**2)**-.5
  nuisance=(2*np.sum((2*np.pi*(freq-fbar))**2*weights)/sigma**2)**-.5
  table.append({'bandwidth_GHz':B/1e9,'target_model':target,'source_formula_ps':old*1e12,'known_phase_ps':known*1e12,'nuisance_phase_ps':nuisance*1e12})
with (OUT/'corrected_t2.csv').open('w') as f:
 writer=csv.DictWriter(f,fieldnames=table[0]);writer.writeheader();writer.writerows(table)
result={'seed':seed,'numpy':np.__version__,'trials_per_cell':trials,'samples':n,'noise':noise_stats,'monte_carlo':rows,'runtime_seconds':time.monotonic()-started,'local_initialization':'true delay +0.01; validates local precision only'}
(OUT/'results.json').write_text(json.dumps(result,indent=2));(OUT/'RUN_LOG.txt').write_text(json.dumps(result,indent=2)+'\n');print(json.dumps(result,indent=2))