# Setup: dependencies, the shared helpers, non-interactive plotting.
import importlib.util
import subprocess
import sys
import time
import warnings
from pathlib import Path

_needed = ("mne", "scipy", "matplotlib", "pandas", "pooch")
_missing = [p for p in _needed if importlib.util.find_spec(p) is None]
if _missing:
    _req = next((d / "requirements.txt" for d in (Path.cwd(), *Path.cwd().parents)
                 if (d / "requirements.txt").exists()), None)
    _cmd = [sys.executable, "-m", "pip", "install", "-q"]
    _cmd += ["-r", str(_req)] if _req else ["mne==1.10.2", "pooch>=1.8"]
    subprocess.check_call(_cmd)

_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
                if (d / "_shared" / "helpers_l6.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L6/ (or notebooks/) so that _shared/helpers_l6.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l1 as L1
import helpers_l3 as L3
import helpers_l6 as L6

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import mne
from scipy import stats

mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72

# Iteration counts.  A power simulation is linear in its simulation count and every download is
# linear in the cohort size, so both are budgeted here and FULL_RUN says what a thorough run costs.
FULL_RUN = False
N_SIMS = 2000 if not FULL_RUN else 20000        # simulated studies per grid cell
N_SIMS_BETWEEN = 4000 if not FULL_RUN else 20000
N_PILOT = 10 if not FULL_RUN else 49            # ds-iowapd subjects PER GROUP (49 controls is the cohort cap)
N_POOL = 36 if not FULL_RUN else 120            # ds-dortmund subjects in the subsampling pool
N_BOOT = 400 if not FULL_RUN else 4000          # bootstrap replicates of a pilot
SEED = L6.SEED
ALPHA = L6.ALPHA
TARGET_POWER = 0.80

print(f"MNE {mne.__version__}; helpers_l6 imported from notebooks/_shared")
print(f"FULL_RUN = {FULL_RUN}: {N_SIMS} simulated studies per power cell, {N_PILOT} ds-iowapd subjects per "
      f"group, {N_POOL} ds-dortmund subjects in the pool, {N_BOOT} bootstrap replicates")
print(f"Downloads: ~{N_PILOT * 2 * 36} MB of ds-iowapd and ~{N_POOL * 24} MB of ds-dortmund, fetched one "
      f"subject at a time and deleted before the next is fetched, so at most one recording is ever on disk")

SUBJECTS = list(L6.SUBSET_DEFAULT)
L6.disk_report("before any download")
t0 = time.time()
cohort = L6.load_cohort(SUBJECTS)
print(f"{len(cohort)} subjects in {time.time() - t0:.0f} s")
L6.disk_report("after the ERP CORE loop")

i_prespec = L6.path_index(L6.PRESPECIFIED)
MIN_TRIALS = int(L6.FIXED_CHOICES["min_trials_per_condition"])
rows = []
for s in sorted(cohort):
    p = cohort[s]
    keep = p["keep"][i_prespec]
    t_scores = p["amp"][i_prespec][keep & p["is_target"]].astype(float)
    s_scores = p["amp"][i_prespec][keep & ~p["is_target"]].astype(float)
    if min(len(t_scores), len(s_scores)) < MIN_TRIALS:
        print(f"EXCLUDED {s}: {len(t_scores)} target / {len(s_scores)} standard trials survive "
              f"(stated minimum {MIN_TRIALS})")
        continue
    rows.append({"subject": s, "n_target": len(t_scores), "n_standard": len(s_scores),
                 "diff_uv": t_scores.mean() - s_scores.mean(),
                 "sd_target": t_scores.std(ddof=1), "sd_standard": s_scores.std(ddof=1),
                 "sme_uv": np.sqrt(t_scores.var(ddof=1) / len(t_scores)
                                   + s_scores.var(ddof=1) / len(s_scores))})
pilot = pd.DataFrame(rows)
mu = float(pilot.diff_uv.mean())
sd_obs = float(pilot.diff_uv.std(ddof=1))
mean_meas_var = float((pilot.sme_uv ** 2).mean())
tau2 = sd_obs ** 2 - mean_meas_var
tau = float(np.sqrt(max(tau2, 0.0)))
sigma = float(np.sqrt(np.average(np.r_[pilot.sd_target ** 2, pilot.sd_standard ** 2],
                                 weights=np.r_[pilot.n_target - 1, pilot.n_standard - 1])))
print(pilot.round(3).to_string(index=False))
print(f"\nVARIANCE DECOMPOSITION ({len(pilot)} subjects, P3 mean amplitude 300-600 ms at Pz, "
      f"target minus standard):")
print(f"   mu    group effect                              {mu:+.4f} uV")
print(f"   observed SD of the subject scores              {sd_obs:.4f} uV")
print(f"   mean measurement variance (SME^2)              {mean_meas_var:.4f} uV^2  "
      f"(mean SME {pilot.sme_uv.mean():.3f} uV)")
print(f"   tau   between-subject SD of the TRUE effect    {tau:.4f} uV   "
      f"= sqrt({sd_obs:.4f}^2 - {mean_meas_var:.4f})")
print(f"   sigma trial-to-trial SD (pooled within-subject) {sigma:.4f} uV")
print()
print(f"   measurement noise accounts for {100 * mean_meas_var / sd_obs ** 2:.1f} % of the observed "
      f"between-subject variance; the remaining {100 * tau2 / sd_obs ** 2:.1f} % is real "
      f"between-subject difference that no number of trials will remove.")
print(f"   Cohen's dz on the observed scores: {mu / sd_obs:.3f}.  With an infinite number of trials per "
      f"subject it would rise to mu/tau = {mu / tau:.3f}, which is the ceiling this design has.")

SUBJECT_COUNTS = [6, 8, 10, 12, 15, 20, 25, 30, 40]
TRIAL_COUNTS = [5, 10, 15, 20, 30, 40, 60, 100]
t0 = time.time()
grid = L6.power_grid(effect_uv=mu, between_sd_uv=tau, trial_sd_uv=sigma,
                     subject_counts=SUBJECT_COUNTS, trial_counts=TRIAL_COUNTS,
                     trial_ratio=4.0, n_sims=N_SIMS, alpha=ALPHA, seed=SEED)
_dt = time.time() - t0
print(f"{len(SUBJECT_COUNTS)} x {len(TRIAL_COUNTS)} = {grid.size} cells x {N_SIMS} simulated studies "
      f"= {grid.size * N_SIMS:,} t-tests in {_dt:.2f} s")
print(f"(a full run at 20,000 simulations per cell would take about {_dt * 20000 / N_SIMS:.1f} s -- the "
      f"simulation is the cheap part of this notebook; the downloads in sections 3 and 4 are not.  The "
      f"Monte-Carlo standard error on a power near 0.8 is sqrt(.8 x .2 / n) = "
      f"{np.sqrt(0.8 * 0.2 / N_SIMS):.4f} at {N_SIMS} and {np.sqrt(0.8 * 0.2 / 20000):.4f} at 20,000)")

small = L6.smallest_design(grid, SUBJECT_COUNTS, TRIAL_COUNTS, target=TARGET_POWER)
print(f"\nex-6-3 GRID READ-OFF -- smallest subject count reaching {TARGET_POWER:.0%} power, "
      f"per target-trial count:")
print(f"{'target trials':>14s} {'standards':>10s} {'subjects':>9s} {'power there':>12s}")
for r in small:
    print(f"{r['n_trials_rare']:14d} {r['n_trials_rare'] * 4:10d} "
          f"{str(r['n_subjects']) if r['reached'] else '> ' + str(SUBJECT_COUNTS[-1]):>9s} "
          f"{r['power']:12.3f}")
print()
print(f"Read that table with one grid step of slack.  A cell whose power lands near {TARGET_POWER:.2f} can "
      f"cross the line either way between runs: with FULL_RUN = True (20,000 simulations instead of "
      f"{N_SIMS}) one row of this table moves on this machine -- 30 target trials asks for 10 subjects "
      f"rather than 8, because the 8-subject cell sits at 0.802.  Every other row is unchanged.  The "
      f"answer to 'how many subjects' is a region of this grid, not a cell.")

fig, ax = L6.plot_power_grid(grid, SUBJECT_COUNTS, TRIAL_COUNTS, target=TARGET_POWER,
                             title=f"Simulated power, P3 target − standard at Pz "
                                   f"(µ = {mu:+.2f} µV, τ = {tau:.2f} µV, σ = {sigma:.2f} µV, "
                                   f"{N_SIMS} studies per cell)")
plt.show()   # render the static figure(s) of this cell inline

print("Where the two sample sizes stop helping:")
j40 = TRIAL_COUNTS.index(40)
for i, ns in enumerate(SUBJECT_COUNTS):
    print(f"   {ns:3d} subjects: power {grid[i, 0]:.3f} at {TRIAL_COUNTS[0]:3d} target trials -> "
          f"{grid[i, j40]:.3f} at 40 -> {grid[i, -1]:.3f} at {TRIAL_COUNTS[-1]}")
print()
print(f"Read the rows, not the cells.  Adding trials moves power a long way while the measurement noise is "
      f"still large compared with tau, and then stops: past about "
      f"{TRIAL_COUNTS[int(np.argmax(grid[len(SUBJECT_COUNTS) // 2] > 0.95 * grid[len(SUBJECT_COUNTS) // 2].max()))]} "
      f"target trials the remaining variance is between-subject variance, which only subjects can buy.")
print(f"The ceiling for an infinite number of trials, at each subject count:")
inf = [L6.simulate_power(effect_uv=mu, between_sd_uv=tau, trial_sd_uv=0.0, n_subjects=ns,
                         n_trials_rare=1, n_sims=N_SIMS, alpha=ALPHA,
                         rng=np.random.default_rng(SEED + 7)) for ns in SUBJECT_COUNTS]
print("   " + "  ".join(f"{ns}:{p:.2f}" for ns, p in zip(SUBJECT_COUNTS, inf)))
print(f"\nThis dataset's own design ({int(pilot.n_target.median())} median surviving target trials, "
      f"{len(pilot)} subjects) sits at simulated power "
      f"{L6.simulate_power(effect_uv=mu, between_sd_uv=tau, trial_sd_uv=sigma, n_subjects=len(pilot), n_trials_rare=int(pilot.n_target.median()), n_sims=N_SIMS, alpha=ALPHA, rng=np.random.default_rng(SEED + 11)):.3f} "
      f"for an effect of {mu:+.2f} uV.")

print(f"Expected half-width of the 95 % CI (uV), by design, for tau = {tau:.3f} and sigma = {sigma:.3f}:")
print(f"{'subjects':>9s} " + "".join(f"{t:>9d}" for t in TRIAL_COUNTS) + "   <- target trials")
for ns in SUBJECT_COUNTS:
    widths = []
    for nt in TRIAL_COUNTS:
        sd_score = np.sqrt(tau ** 2 + sigma ** 2 * (1 / nt + 1 / (4 * nt)))
        widths.append(stats.t.ppf(1 - ALPHA / 2, ns - 1) * sd_score / np.sqrt(ns))
    print(f"{ns:9d} " + "".join(f"{w:9.3f}" for w in widths))
print()
print(f"At {len(pilot)} subjects and {int(pilot.n_target.median())} target trials the expected half-width is "
      f"about {stats.t.ppf(1 - ALPHA / 2, len(pilot) - 1) * np.sqrt(tau ** 2 + sigma ** 2 * (1 / pilot.n_target.median() + 1 / pilot.n_standard.median())) / np.sqrt(len(pilot)):.2f} uV "
      f"on an effect of {mu:.2f} uV -- so the study can say the effect is there and can barely say how big.")

for k, v in L6.RESTING_MEASURE.items():
    print(f"   {k:14s} : {v}")
print()
print("ds-iowapd, from data/directory.yaml:")
for k in ("device", "sfreq", "n_channels", "population", "paradigm", "reference", "online_filters",
          "mains_hz", "license", "access"):
    print(f"   {k:14s} : {L6.DATASETS_L6['ds-iowapd'][k]}")
print("   caveats        :")
for c in L6.DATASETS_L6["ds-iowapd"]["caveats"]:
    print(f"                    - {c}")

parts = L6.iowapd_participants()
print(f"participants.tsv: {len(parts)} rows, columns {list(parts.columns)}")
print(parts.GROUP.value_counts().to_string())
print(f"\nage: PD {parts.loc[parts.GROUP == 'PD', 'AGE'].mean():.1f} +- "
      f"{parts.loc[parts.GROUP == 'PD', 'AGE'].std():.1f}, "
      f"Control {parts.loc[parts.GROUP == 'Control', 'AGE'].mean():.1f} +- "
      f"{parts.loc[parts.GROUP == 'Control', 'AGE'].std():.1f}")
# Selection rule, stated before any recording is read: the first N_PILOT of each group in file order.
pd_ids = parts.loc[parts.GROUP == "PD", "participant_id"].tolist()[:N_PILOT]
hc_ids = parts.loc[parts.GROUP == "Control", "participant_id"].tolist()[:N_PILOT]
print(f"\nselection rule (stated before any recording is read): the first {N_PILOT} PD and the first "
      f"{N_PILOT} Control participants in participants.tsv order.  No recording is inspected first, and "
      f"no subject is swapped for another afterwards.")
print(f"   PD      : {', '.join(pd_ids)}")
print(f"   Control : {', '.join(hc_ids)}")

t0 = time.time()
free0 = L6.disk_report("before the ds-iowapd loop",
                       folders={"course downloads": L1.download_dir()})["free_gb"]
iowa = L6.resting_cohort("ds-iowapd", pd_ids + hc_ids, progress=True)
iowa["group"] = ["PD"] * len(pd_ids) + ["Control"] * len(hc_ids)
print(f"\n{int(iowa.ok.sum())} of {len(iowa)} recordings measured in {time.time() - t0:.0f} s")
L6.disk_report("after the ds-iowapd loop", folders={"course downloads": L1.download_dir()})
if (~iowa.ok).any():
    print("failures:")
    print(iowa.loc[~iowa.ok, ["subject", "reason"]].to_string(index=False))

ok = iowa[iowa.ok]
a = ok.loc[ok.group == "PD", "value"].to_numpy()
b = ok.loc[ok.group == "Control", "value"].to_numpy()
g = L6.hedges_g(a, b)
t_ind, p_ind = stats.ttest_ind(a, b, equal_var=False)
print(f"PILOT RESULT -- {L6.RESTING_MEASURE['name']} ({L6.RESTING_MEASURE['definition']}), eyes-open rest:")
print(f"   PD      n = {len(a):2d}, mean {a.mean():.4f}, SD {a.std(ddof=1):.4f}")
print(f"   Control n = {len(b):2d}, mean {b.mean():.4f}, SD {b.std(ddof=1):.4f}")
print(f"   difference {a.mean() - b.mean():+.4f} (PD minus Control)")
print(f"   Welch t({stats.ttest_ind(a, b, equal_var=False).df:.1f}) = {t_ind:.3f}, p = {p_ind:.4f}")
print(f"   Hedges' g = {g['g']:+.4f}, 95% CI [{g['ci'][0]:+.4f}, {g['ci'][1]:+.4f}] "
      f"(Cohen's d {g['d']:+.4f}, correction factor {g['correction_j']:.4f})")
print()
print(f"   The confidence interval is what the sample size question actually depends on, and it spans "
      f"{g['ci'][1] - g['ci'][0]:.2f} standardized units.")

t0 = time.time()
g_small = min(g["ci"], key=abs)            # the CI end nearer zero: the smallest effect still compatible
g_large = max(g["ci"], key=abs)            # the CI end further from zero: the largest
n_point = L6.required_n_between(g["g"], power=TARGET_POWER, alpha=ALPHA, n_sims=N_SIMS_BETWEEN, seed=SEED)
n_small = L6.required_n_between(g_small, power=TARGET_POWER, alpha=ALPHA, n_sims=N_SIMS_BETWEEN, seed=SEED)
n_large = L6.required_n_between(g_large, power=TARGET_POWER, alpha=ALPHA, n_sims=N_SIMS_BETWEEN, seed=SEED)
print(f"REQUIRED N per group for {TARGET_POWER:.0%} power at alpha = {ALPHA}, by simulation "
      f"({N_SIMS_BETWEEN} studies per candidate n, bisection; {time.time() - t0:.1f} s):")
print(f"   from the point estimate       g = {g['g']:+.3f} : "
      f"{n_point if n_point else '> 4000'} per group")
print(f"   from the CI end NEARER zero   g = {g_small:+.3f} : "
      f"{n_small if n_small else '> 4000'} per group   <- the smallest effect the pilot still allows, so "
      f"the LARGEST study it implies")
print(f"   from the CI end FURTHER away  g = {g_large:+.3f} : "
      f"{n_large if n_large else '> 4000'} per group   <- the most optimistic reading")
print()
print(f"   A pilot of {len(a)} and {len(b)} does not give you a sample size.  It gives you a RANGE of sample "
      f"sizes, and the range here is wide enough that the point estimate is close to useless on its own.")

rng = np.random.default_rng(SEED)
boot_g, boot_n = [], []
t0 = time.time()
for _ in range(N_BOOT):
    ga = L6.hedges_g(rng.choice(a, len(a), replace=True), rng.choice(b, len(b), replace=True))["g"]
    boot_g.append(ga)
    boot_n.append(L6.required_n_between(ga, power=TARGET_POWER, alpha=ALPHA, n_sims=600, seed=SEED))
boot_g = np.array(boot_g)
finite = np.array([n for n in boot_n if n is not None], float)
print(f"{N_BOOT} bootstrap pilots of the same size in {time.time() - t0:.0f} s "
      f"(600 simulations per required-N search, so each n is coarse; the spread, not the digits, is the point)")
print(f"   bootstrap g      : median {np.median(boot_g):+.3f}, 2.5-97.5 % "
      f"[{np.percentile(boot_g, 2.5):+.3f}, {np.percentile(boot_g, 97.5):+.3f}]")
print(f"   sign flips in    : {100 * np.mean(np.sign(boot_g) != np.sign(g['g'])):.1f} % of bootstrap pilots")
print(f"   required N       : median {np.median(finite):.0f}, 2.5-97.5 % "
      f"[{np.percentile(finite, 2.5):.0f}, {np.percentile(finite, 97.5):.0f}] per group"
      f" ({len(boot_n) - len(finite)} of {N_BOOT} pilots implied more than 4000 per group)")
print(f"   ratio of the 97.5th to the 2.5th percentile: "
      f"{np.percentile(finite, 97.5) / max(np.percentile(finite, 2.5), 1):.0f}x")

for k in ("device", "sfreq", "n_channels", "population", "paradigm", "reference", "online_filters",
          "mains_hz", "license", "license_note", "access"):
    print(f"   {k:14s} : {L6.DATASETS_L6['ds-dortmund'][k]}")
print("   caveats        :")
for c in L6.DATASETS_L6["ds-dortmund"]["caveats"]:
    print(f"                    - {c}")

dparts = L6.dortmund_participants()
print(f"\nparticipants.tsv: {len(dparts)} rows, columns {list(dparts.columns)}")
print(f"age {dparts.age.min()}-{dparts.age.max()}, median {dparts.age.median():.0f}; "
      f"sex {dparts.sex.value_counts().to_dict()}")
# Selection rule, stated before any recording is read.
pool_ids = dparts.loc[dparts.session1 == "yes", "participant_id"].tolist()[:N_POOL]
ages = dparts.set_index("participant_id").loc[pool_ids, "age"]
cut = float(ages.median())
print(f"\nselection rule (stated before any recording is read): the first {N_POOL} participants with a "
      f"ses-1 recording, split at the pool's median age ({cut:.0f} years) into 'younger' and 'older'.  "
      f"The split is a median split for the arithmetic's sake, not a claim about ageing.")
print(f"   pool ages: {ages.min()}-{ages.max()}, median {cut:.0f}")

t0 = time.time()
L6.disk_report("before the ds-dortmund loop", folders={"course downloads": L1.download_dir()})
dort = L6.resting_cohort("ds-dortmund", pool_ids, task="EyesOpen", acq="pre", progress=True)
dort["age"] = ages.to_numpy()
dort["group"] = np.where(dort.age <= cut, "younger", "older")
print(f"\n{int(dort.ok.sum())} of {len(dort)} recordings measured in {time.time() - t0:.0f} s")
L6.disk_report("after the ds-dortmund loop", folders={"course downloads": L1.download_dir()})
if (~dort.ok).any():
    print(dort.loc[~dort.ok, ["subject", "reason"]].to_string(index=False))

dok = dort[dort.ok]
ya = dok.loc[dok.group == "younger", "value"].to_numpy()
oa = dok.loc[dok.group == "older", "value"].to_numpy()
g_pool = L6.hedges_g(oa, ya)
t_pool, p_pool = stats.ttest_ind(oa, ya, equal_var=False)
n_pool = L6.required_n_between(g_pool["g"], power=TARGET_POWER, alpha=ALPHA, n_sims=N_SIMS_BETWEEN, seed=SEED)
print(f"THE POOL ({len(oa)} older, {len(ya)} younger), same measure:")
print(f"   older   mean {oa.mean():.4f}, SD {oa.std(ddof=1):.4f}")
print(f"   younger mean {ya.mean():.4f}, SD {ya.std(ddof=1):.4f}")
print(f"   Welch t = {t_pool:.3f}, p = {p_pool:.4f}; Hedges' g = {g_pool['g']:+.4f} "
      f"[{g_pool['ci'][0]:+.4f}, {g_pool['ci'][1]:+.4f}]")
print(f"   required N from the pool's own estimate: "
      f"{n_pool if n_pool else '> 4000'} per group")
print()
print(f"   The pool is {len(dok)} people, which is itself a small study -- it is the best available truth "
      f"here, not the truth.  Everything below compares pilots with IT, not with reality.")

rng = np.random.default_rng(SEED + 3)
# A pilot size must be strictly smaller than the pool's group size, or every "independent pilot"
# is the same pilot and the spread collapses to zero.
PILOT_SIZES = [k for k in ([4, 8, 12] if not FULL_RUN else [5, 10, 20, 40])
               if k < min(len(ya), len(oa))]
N_PILOTS = 200 if not FULL_RUN else 1000
t0 = time.time()
sub_rows = []
for k in PILOT_SIZES:
    gs, ns = [], []
    for _ in range(N_PILOTS):
        ia = rng.choice(len(oa), k, replace=False)
        ib = rng.choice(len(ya), k, replace=False)
        gg = L6.hedges_g(oa[ia], ya[ib])["g"]
        gs.append(gg)
        ns.append(L6.required_n_between(gg, power=TARGET_POWER, alpha=ALPHA, n_sims=400, seed=SEED))
    gs = np.array(gs)
    fin = np.array([v for v in ns if v is not None], float)
    sub_rows.append({"n_per_group": k, "g_median": float(np.median(gs)),
                     "g_lo": float(np.percentile(gs, 2.5)), "g_hi": float(np.percentile(gs, 97.5)),
                     "sign_flip": float(np.mean(np.sign(gs) != np.sign(g_pool["g"]))),
                     "n_median": float(np.median(fin)) if fin.size else np.nan,
                     "n_lo": float(np.percentile(fin, 2.5)) if fin.size else np.nan,
                     "n_hi": float(np.percentile(fin, 97.5)) if fin.size else np.nan,
                     "n_unbounded": int(len(ns) - len(fin)), "n_pilots": len(ns),
                     "gs": gs, "ns": fin})
print(f"{len(sub_rows)} pilot sizes x {N_PILOTS} pilots drawn without replacement from "
      f"{len(oa)} older and {len(ya)} younger in {time.time() - t0:.0f} s")
print(f"(the pilots overlap, because the pool is finite; the spread below is therefore a LOWER bound on "
      f"what genuinely independent pilots would show)")
print(f"\npool: g = {g_pool['g']:+.3f}, required N = {n_pool if n_pool else '> 4000'} per group\n")
print(f"{'pilot n/group':>14s} {'median g':>9s} {'g 2.5-97.5%':>20s} {'sign flip':>10s} "
      f"{'median req N':>13s} {'req N 2.5-97.5%':>18s} {'no answer':>10s}")
for r in sub_rows:
    print(f"{r['n_per_group']:14d} {r['g_median']:+9.3f} "
          f"{r['g_lo']:+8.3f} to {r['g_hi']:+7.3f} {r['sign_flip']:10.1%} "
          f"{r['n_median']:13.0f} {r['n_lo']:8.0f} to {r['n_hi']:7.0f} "
          f"{r['n_unbounded']:4d}/{r['n_pilots']:<5d}")
print()
print(f"'no answer' counts pilots whose estimated effect was so near zero that no n up to 4000 per group "
      f"reached {TARGET_POWER:.0%} power.  Those are not failures of the method: they are what a pilot that "
      f"sampled badly actually tells you, and dropping them is how the published required-N estimates get "
      f"optimistic.")

fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.3))
for r in sub_rows:
    axes[0].hist(r["gs"], bins=30, alpha=0.5, label=f"pilot n = {r['n_per_group']}/group", density=True)
axes[0].axvline(g_pool["g"], color="k", lw=2, label=f"pool g = {g_pool['g']:+.3f}")
axes[0].axvline(0, color="gray", lw=0.8)
axes[0].set(xlabel="Hedges' g estimated by the pilot (dimensionless)", ylabel="Density",
            title=f"{N_PILOTS} independent pilots per size, drawn from a {len(dok)}-person pool")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3)

pos = np.arange(len(sub_rows))
axes[1].boxplot([np.clip(r["ns"], 1, 4000) for r in sub_rows], positions=pos, widths=0.55,
                showfliers=False, patch_artist=True,
                boxprops=dict(facecolor="tab:blue", alpha=0.45))
if n_pool:
    axes[1].axhline(n_pool, color="k", lw=2, label=f"from the pool: {n_pool}/group")
axes[1].set_yscale("log")
axes[1].set_xticks(pos, [f"{r['n_per_group']}/group" for r in sub_rows])
axes[1].set(xlabel="Pilot size", ylabel="Required N per group for 80 % power (count, log scale)",
            title="What each pilot would have told you to plan for")
axes[1].legend(fontsize=8)
axes[1].grid(alpha=0.3, axis="y")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print("ds-bonn, from data/directory.yaml -- discussed, never derived from:")
for k in ("population", "structure", "reference", "online_filters", "license", "license_note",
          "paper_doi", "source"):
    print(f"   {k:14s} : {L6.DATASETS_L6['ds-bonn'][k]}")
print("   caveats        :")
for c in L6.DATASETS_L6["ds-bonn"]["caveats"]:
    print(f"                    - {c}")
print()
print(L6.BONN_POLICY)
print()
print("The arithmetic of the mistake, with no data required:")
for n_people, n_seg in ((10, 500), (5, 100)):
    se_seg = 1 / np.sqrt(n_seg)
    se_ppl = 1 / np.sqrt(n_people)
    print(f"   {n_seg} segments from {n_people} people: a standard error computed on the segments is "
          f"{se_seg:.3f} sigma, on the people {se_ppl:.3f} sigma -- too small by a factor of "
          f"{se_ppl / se_seg:.1f} if the within-person correlation is perfect.")
print("   The truth lies between the two and depends on the intraclass correlation, which is why the fix is "
      "to model the grouping (nb-6-2-lmm) rather than to pick one of the two numbers.")

print("nb-6-3-power-sim -- L6.3 numbers (draft; TODO(confirm) at author review)")
print()
print(f"1. VARIANCE DECOMPOSITION -- ds-erpcore P3, {len(pilot)} subjects, mean amplitude 300-600 ms at Pz "
      f"under helpers_l6.PRESPECIFIED:")
print(f"     mu (group effect)                    {mu:+.4f} uV")
print(f"     tau (between-subject SD, true)       {tau:.4f} uV")
print(f"     sigma (trial-to-trial SD)            {sigma:.4f} uV")
print(f"     observed SD of subject scores        {sd_obs:.4f} uV; mean SME {pilot.sme_uv.mean():.4f} uV")
print(f"     dz observed {mu / sd_obs:.3f}; ceiling with infinite trials mu/tau = {mu / tau:.3f}")
print()
print(f"2. ex-6-3 POWER GRID ({N_SIMS} simulated studies per cell, seed {SEED}, 4 standards per target, "
      f"alpha = {ALPHA}) -- smallest subject count for {TARGET_POWER:.0%} power:")
for r in small:
    print(f"     {r['n_trials_rare']:3d} target trials -> "
          f"{str(r['n_subjects']) if r['reached'] else '> ' + str(SUBJECT_COUNTS[-1])} subjects "
          f"(power {r['power']:.3f})")
print()
print(f"3. ds-iowapd PILOT ({len(a)} PD, {len(b)} Control, first-in-file-order rule), "
      f"{L6.RESTING_MEASURE['name']}:")
print(f"     PD {a.mean():.4f} +- {a.std(ddof=1):.4f}; Control {b.mean():.4f} +- {b.std(ddof=1):.4f}")
print(f"     Hedges' g {g['g']:+.4f} [{g['ci'][0]:+.4f}, {g['ci'][1]:+.4f}]; Welch p = {p_ind:.4f}")
print(f"     required N per group for {TARGET_POWER:.0%} power: {n_point if n_point else '> 4000'} from the "
      f"point estimate; {n_large if n_large else '> 4000'} (optimistic CI end) to "
      f"{n_small if n_small else '> 4000'} (pessimistic CI end)")
print(f"     bootstrap over {N_BOOT} resampled pilots: median required N {np.median(finite):.0f}, "
      f"2.5-97.5 % [{np.percentile(finite, 2.5):.0f}, {np.percentile(finite, 97.5):.0f}], "
      f"sign flip in {100 * np.mean(np.sign(boot_g) != np.sign(g['g'])):.1f} % of resamples")
print()
print(f"4. ds-dortmund SUBSAMPLING (pool of {len(dok)}, median-age split at {cut:.0f} y, "
      f"{N_PILOTS} independent pilots per size):")
print(f"     pool: g = {g_pool['g']:+.4f} [{g_pool['ci'][0]:+.4f}, {g_pool['ci'][1]:+.4f}], "
      f"required N {n_pool if n_pool else '> 4000'} per group")
for r in sub_rows:
    print(f"     pilot n = {r['n_per_group']:2d}/group: median g {r['g_median']:+.3f} "
          f"[{r['g_lo']:+.3f}, {r['g_hi']:+.3f}], sign flip {r['sign_flip']:.1%}, "
          f"median required N {r['n_median']:.0f} [{r['n_lo']:.0f}, {r['n_hi']:.0f}], "
          f"{r['n_unbounded']} of {r['n_pilots']} gave no bounded answer")
print()
print(f"5. ds-bonn: discussed, never derived from (licence informal-academic, snippets: no).  500 segments "
      f"from 10 people.")
print()
print("TODO(confirm): the resting measure here is a course decision (helpers_l6.RESTING_MEASURE), not a "
      "value taken from either paper.  Neither the ds-iowapd nor the ds-dortmund result should be read as a "
      "replication of anything those papers report -- the catalog carries no published numbers to compare "
      "with, and the pipelines are not theirs.")
print()
print("Pitfall: pf-site-device-confound.  No widget for this lesson.")