# (c) 2026 William J. Dwyer. Free and open-source software under the MIT License (see LICENSE).
# If you use or build upon this software, please retain this notice and cite the accompanying
# methodology paper together with the archived software (Zenodo DOI: 10.5281/zenodo.20709963).
"""
CV RATIO and the recovery of parametric inference -- detailed simulation.

Premise. Under log-normality CV = sqrt(exp(sigma_e^2)-1) depends only on the log-scale SD,
so the between-group CV ratio is a one-to-one readout of the log-scale variance ratio:
  sigma_e^2 = ln(CV^2 + 1)   =>   rho := (sigma_e,dz/sigma_e,o)^2 = ln(CV_dz^2+1)/ln(CV_o^2+1).
CV ratio = 1  <=> equal log-variances <=> pooled (Student) log-t exact/UMPU;
CV ratio != 1 <=> log-scale Behrens-Fisher <=> Welch's log-t is the calibrated test.

The CV ratio's job is the VARIANCE MODEL (pooled vs Welch); it is distinct from the SCREEN,
whose job is the SCALE (raw vs log). This script characterizes the ratio's job in detail:
  A. size (type-I) of pooled / Welch / CV-ratio-adaptive across variance ratio x n-allocation;
  B. 95% CI coverage for the geometric-mean ratio (the estimation side of "recovery");
  C. power, and the power cost of defaulting to Welch when variances are actually equal;
  D. operating characteristics of the CV-ratio gate, and sensitivity to the gate width;
  E. the sample CV ratio as an estimator of the true ratio (bias/SD by n);
  F. the failure of a summaries-only reconstruction of the whole test (an honest boundary).
"""
import numpy as np, json
from scipy import stats
import warnings; warnings.filterwarnings("ignore")

rng = np.random.default_rng(20260608)
REPS = 8000
ALPHA = 0.05
SE_O = 0.833                       # control log-SD -> CV_o = 1.0 (manuscript anchor)
def cv_of(se): return float(np.sqrt(np.exp(se**2)-1.0))
CV_O = cv_of(SE_O)
LN2 = np.log(2.0)

# log-variance ratios rho and the CV ratios they induce (rho=2.62 is the 1.5+/-3.4 anchor)
RHO = [0.25, 0.50, 1.00, 1.74, 2.62, 3.69]
def se_dz_of(rho): return SE_O*np.sqrt(rho)
def cvratio_of(rho): return cv_of(se_dz_of(rho))/CV_O
ALLOC = [(10,10),(20,20),(50,50),(10,30),(30,10),(12,48),(48,12)]
GATE = 1.5

tcrit = lambda df: stats.t.ppf(0.975, df)

def analyze(go, gd, n_o, n_d, true_shift, gate=GATE):
    lo, ld = np.log(go), np.log(gd)
    mo, md = lo.mean(1), ld.mean(1)
    vo, vd = lo.var(1, ddof=1), ld.var(1, ddof=1)
    diff = md - mo                                   # estimate of ln(GMR)
    # pooled (equal-variance)
    sp2 = ((n_d-1)*vd + (n_o-1)*vo)/(n_o+n_d-2)
    se_p = np.sqrt(sp2*(1.0/n_o + 1.0/n_d)); df_p = n_o+n_d-2
    p_p = 2*stats.t.sf(np.abs(diff/se_p), df_p)
    cov_p = np.abs(diff-true_shift) <= tcrit(df_p)*se_p
    # Welch (separate-variance)
    a, b = vd/n_d, vo/n_o
    se_w = np.sqrt(a+b)
    df_w = (a+b)**2 / (a**2/(n_d-1) + b**2/(n_o-1))
    p_w = 2*stats.t.sf(np.abs(diff/se_w), df_w)
    cov_w = np.abs(diff-true_shift) <= tcrit(df_w)*se_w
    # CV-ratio adaptive gate (natural-domain sample CV ratio)
    rhat = (gd.std(1,ddof=1)/gd.mean(1)) / (go.std(1,ddof=1)/go.mean(1))
    usew = (rhat>gate) | (rhat<1.0/gate)
    p_a   = np.where(usew, p_w, p_p)
    cov_a = np.where(usew, cov_w, cov_p)
    return dict(p_p=p_p, p_w=p_w, p_a=p_a, cov_p=cov_p, cov_w=cov_w, cov_a=cov_a,
                rhat=rhat, welch_frac=float(usew.mean()))

def gen(se_dz, n_o, n_d, shift, reps=REPS):
    go = rng.lognormal(0.0,   SE_O,  (reps, n_o))
    gd = rng.lognormal(shift, se_dz, (reps, n_d))
    return go, gd

# ---------- A & B & C: size, coverage, power across the grid ----------
grid=[]
for rho in RHO:
    se_dz = se_dz_of(rho); cvr = cvratio_of(rho)
    for (n_o,n_d) in ALLOC:
        gN = gen(se_dz,n_o,n_d,0.0)       # null  -> size
        gA = gen(se_dz,n_o,n_d,LN2)       # alt   -> power, coverage
        rN = analyze(*gN, n_o,n_d, 0.0)
        rA = analyze(*gA, n_o,n_d, LN2)
        size = lambda r,k: float(np.mean(r[k]<ALPHA))
        grid.append(dict(rho=rho, cvratio=round(cvr,2), n_o=n_o, n_d=n_d,
            size_pool=round(size(rN,'p_p'),3), size_welch=round(size(rN,'p_w'),3), size_adapt=round(size(rN,'p_a'),3),
            cov_pool=round(float(rA['cov_p'].mean()),3), cov_welch=round(float(rA['cov_w'].mean()),3), cov_adapt=round(float(rA['cov_a'].mean()),3),
            pow_pool=round(size(rA,'p_p'),3), pow_welch=round(size(rA,'p_w'),3), pow_adapt=round(size(rA,'p_a'),3),
            welch_frac=round(rN['welch_frac'],2)))
json.dump(grid, open("cvratio_recover.json","w"), indent=1)

def band(key):
    v=[g[key] for g in grid]; return min(v),max(v),float(np.mean(v))

print("="*72); print("A. SIZE (type-I, nominal 0.05) across log-variance ratio x n-allocation")
print("="*72)
print(f"{'CVratio':>8}{'n_o/n_d':>10}{'pooled':>9}{'Welch':>8}{'adaptive':>10}")
for g in grid:
    if g['rho'] in (1.0, 2.62, 3.69) and (g['n_o'],g['n_d']) in [(20,20),(10,30),(30,10),(12,48),(48,12)]:
        print(f"{g['cvratio']:>8}{str(g['n_o'])+'/'+str(g['n_d']):>10}{g['size_pool']:>9}{g['size_welch']:>8}{g['size_adapt']:>10}")
print(f"  size bands over FULL grid:  pooled {band('size_pool')[0]:.3f}-{band('size_pool')[1]:.3f} | "
      f"Welch {band('size_welch')[0]:.3f}-{band('size_welch')[1]:.3f} | adaptive {band('size_adapt')[0]:.3f}-{band('size_adapt')[1]:.3f}")

print("\n"+"="*72); print("B. 95% CI COVERAGE for the geometric-mean ratio (recovery of interval estimates)")
print("="*72)
print(f"  coverage bands: pooled {band('cov_pool')[0]:.3f}-{band('cov_pool')[1]:.3f} | "
      f"Welch {band('cov_welch')[0]:.3f}-{band('cov_welch')[1]:.3f} | adaptive {band('cov_adapt')[0]:.3f}-{band('cov_adapt')[1]:.3f}")
for g in grid:
    if g['rho']==2.62 and (g['n_o'],g['n_d']) in [(30,10),(10,30),(20,20)]:
        print(f"   CVratio {g['cvratio']} n {g['n_o']}/{g['n_d']}: pooled cov {g['cov_pool']}  Welch cov {g['cov_welch']}  adaptive cov {g['cov_adapt']}")

print("\n"+"="*72); print("C. POWER cost of defaulting to Welch when variances are EQUAL (rho=1)")
print("="*72)
for g in grid:
    if g['rho']==1.0:
        print(f"   n {g['n_o']}/{g['n_d']}: pooled {g['pow_pool']}  Welch {g['pow_welch']}  (loss {g['pow_pool']-g['pow_welch']:+.3f})")

# ---------- D: gate sensitivity ----------
print("\n"+"="*72); print("D. CV-ratio GATE operating characteristics (size band over grid; power at anchor)")
print("="*72)
for gate in [1.3,1.5,2.0]:
    sizes=[];
    for rho in RHO:
        se_dz=se_dz_of(rho)
        for (n_o,n_d) in ALLOC:
            r=analyze(*gen(se_dz,n_o,n_d,0.0,reps=5000), n_o,n_d,0.0,gate=gate)
            sizes.append(float(np.mean(r['p_a']<ALPHA)))
    # power at the worst liberal anchor (rho=2.62, 30/10) and welch-usage
    ra=analyze(*gen(se_dz_of(2.62),30,10,LN2,reps=5000),30,10,LN2,gate=gate)
    print(f"   gate [1/{gate},{gate}]: adaptive size {min(sizes):.3f}-{max(sizes):.3f}  "
          f"| anchor(30/10) power {np.mean(ra['p_a']<ALPHA):.3f}  Welch used {ra['welch_frac']*100:.0f}%")

# ---------- E: sample CV-ratio estimator quality ----------
print("\n"+"="*72); print("E. SAMPLE CV ratio as estimator of the true CV ratio (mean +/- SD)")
print("="*72)
for cvr_true, rho in [(2.27,2.62),(1.0,1.0)]:
    se_dz=se_dz_of(rho)
    for n in [10,20,50,100]:
        go,gd = gen(se_dz,n,n,0.0,reps=6000)
        rh=(gd.std(1,ddof=1)/gd.mean(1))/(go.std(1,ddof=1)/go.mean(1))
        print(f"   true CVratio {cvr_true:.2f}, n={n:>3}: estimated {rh.mean():.2f} +/- {rh.std():.2f}")

# ---------- C2: small-n equal-variance power (where pooling could matter) ----------
print("\n"+"="*72); print("C2. small-n power at EQUAL log-variance (balanced): pooled vs Welch")
print("="*72)
for n in [5,6,8,10,15,20]:
    go,gd = gen(SE_O,n,n,LN2,reps=8000)
    r=analyze(go,gd,n,n,LN2)
    print(f"   n={n:>2}/group: pooled {np.mean(r['p_p']<ALPHA):.3f}  Welch {np.mean(r['p_w']<ALPHA):.3f}")

# ---------- F: summaries-only reconstruction (honest boundary) ----------
print("\n"+"="*72); print("F. summaries-only reconstruction of the test (method of moments): unstable")
print("="*72)
def recon(go,gd,n_o,n_d):
    mo,so=go.mean(1),go.std(1,ddof=1); md,sd=gd.mean(1),gd.std(1,ddof=1)
    s2o=np.log((so/mo)**2+1); s2d=np.log((sd/md)**2+1)
    diff=(np.log(md)-s2d/2)-(np.log(mo)-s2o/2)
    a,b=s2d/n_d,s2o/n_o; se=np.sqrt(a+b)
    df=(a+b)**2/(a**2/(n_d-1)+b**2/(n_o-1))
    return 2*stats.t.sf(np.abs(diff/se),df)
recs=[]
for rho in RHO:
    se_dz=se_dz_of(rho)
    for (n_o,n_d) in ALLOC:
        p=recon(*gen(se_dz,n_o,n_d,0.0,reps=5000),n_o,n_d)
        recs.append(float(np.mean(p<ALPHA)))
print(f"   reconstructed-from-summaries type-I band over grid: {min(recs):.3f}-{max(recs):.3f}")
print("\nsaved cvratio_recover.json")
