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

# 1. Dependencies are pinned in notebooks/requirements.txt.  Nothing is installed when the
#    pinned stack is already present (local runs, CI); a fresh Colab or Binder kernel installs
#    it once.  On Colab, run from a clone of the repository so that notebooks/_shared/ is
#    available (repository URL: TODO(confirm), spec section 13 item 3).
_needed = ('mne', 'scipy', 'matplotlib', '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)

# 2. Shared helpers, located relative to the working directory -- notebooks/<level>/ or
#    notebooks/ -- never through an absolute path.
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
                if (d / "_shared" / "helpers_l4.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L4/ (or notebooks/) so that "
                            "_shared/helpers_l4.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l1
import helpers_l4 as L4
from scipy import stats as sstats

# 3. Plotting: Jupyter's default inline backend renders static PNGs through Agg (no windows,
#    nothing blocks); outside Jupyter the helpers select Agg.  Every MNE figure is requested
#    with show=False and each figure cell ends with plt.show().
import matplotlib.pyplot as plt
import numpy as np
import mne

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

# 4. Quiet the downloader.  pooch, which MNE uses to fetch datasets, logs
#    "Downloading file '...' from '...' to '<cache directory>'" at INFO, and that last field is an
#    ABSOLUTE PATH from whichever machine executed the notebook.  Absolute paths are not allowed in a
#    stored notebook (scripts/scrub-notebooks.py is a CI gate) and re-executing would put them straight
#    back, so the message is suppressed at the source rather than cleaned up afterwards.  Nothing is
#    hidden by this: every cell below prints the file NAMES it fetched and helpers_l4.Downloads prints
#    free disk before and after.  Please do not delete this as noise.
try:
    import pooch

    pooch.get_logger().setLevel("WARNING")
except Exception:  # pooch absent or its API moved: the scrub script is the backstop
    pass
print(f"MNE {mne.__version__}; helpers_l4 imported from notebooks/_shared")
print(f"downloads go to {helpers_l1.download_dir().name}/ (resolved relative to the working directory, "
      "or $EEG_COURSE_DOWNLOADS) and are deleted at the end of this notebook")

SUBJECTS = list(L4.SUBSET_DEFAULT)
CHANNELS = ["C3", "C4"]

dl = L4.Downloads("nb-4-7").start()
maps, meta = {}, []
for sid in SUBJECTS:
    ep, nfo = L4.load_imagery_epochs(sid, L4.MI_RUNS, dl)
    per = {}
    for cond in ("T1", "T2"):
        X = ep[cond].get_data(picks=CHANNELS) * 1e6
        r = L4.morlet_power(X, sfreq=float(ep.info["sfreq"]), times=ep.times)
        P = r["power"].mean(0)                                   # (2, n_freqs, n_times) raw power
        b = (r["times"] >= L4.TF_BASELINE[0]) & (r["times"] <= L4.TF_BASELINE[1])
        per[cond] = np.stack([100 * (P[k] - P[k][:, b].mean(1, keepdims=True))
                              / P[k][:, b].mean(1, keepdims=True) for k in range(len(CHANNELS))])
    maps[sid] = per
    freqs, times = r["freqs"], r["times"]
    meta.append({"subject": sid, "T1 trials": nfo["n_epochs"].get("T1", 0),
                 "T2 trials": nfo["n_epochs"].get("T2", 0),
                 "max |sample| (uV)": max(nfo["max_abs_uv"].values())})
print(L4.fmt_table(meta, floatfmt="{:.0f}"))
edge = L4.edge_seconds(freqs, L4.TF_CYCLES)
safe = (times >= times[0] + edge.max()) & (times <= times[-1] - edge.max())
print()
print(f"{len(SUBJECTS)} subjects; map {len(freqs)} frequencies x {len(times)} time points; the edge-safe "
      f"window is {times[safe][0]:+.3f}..{times[safe][-1]:+.3f} s "
      f"({safe.sum()} of {len(times)} points, {len(freqs) * safe.sum()} cells enter every test below)")
print(f"Excluded subjects, and why: {', '.join(helpers.EEGBCI_EXCLUDE_DEFAULT)} (inconsistent event "
      f"timestamps, data/directory.yaml) always; {', '.join(helpers.EEGBCI_EXCLUDE_OPTIONAL)} by default. "
      f"None of them is in this subset.")

i_c3, i_c4 = CHANNELS.index("C3"), CHANNELS.index("C4")
lat = {sid: {c: maps[sid][c][i_c3] - maps[sid][c][i_c4] for c in ("T1", "T2")} for sid in SUBJECTS}
X_diff = np.stack([(lat[s]["T1"] - lat[s]["T2"])[:, safe] for s in SUBJECTS])     # (n_subj, n_freq, n_time)
print(f"contrast array for the test: {X_diff.shape} (subjects x frequencies x times), "
      f"units: percentage points of (C3 - C4) lateralisation, left-hand minus right-hand imagery")

fig, axes = plt.subplots(1, 3, figsize=(17.5, 4.0))
for ax, (name, M) in zip(axes, (
        ("T1, left-hand imagery", np.mean([lat[s]["T1"] for s in SUBJECTS], axis=0)),
        ("T2, right-hand imagery", np.mean([lat[s]["T2"] for s in SUBJECTS], axis=0)),
        ("T1 minus T2 (the contrast)", np.mean([lat[s]["T1"] - lat[s]["T2"] for s in SUBJECTS], axis=0)))):
    L4.plot_tfr(M, freqs, times, ax=ax, vlim=(-40, 40), edge_s=edge, baseline=L4.TF_BASELINE,
                cbar_label="C3 - C4 lateralisation (percentage points)",
                title=f"{name}\ngrand average over {len(SUBJECTS)} subjects (percentage points)")
fig.suptitle("Lateralisation of the mu/beta desynchronisation, (C3 - C4) in percent-change units", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

t_obs = sstats.ttest_1samp(X_diff, 0, axis=0).statistic
p_unc = sstats.ttest_1samp(X_diff, 0, axis=0).pvalue
n_cells = t_obs.size
print(f"Uncorrected paired t-tests, one per cell: {n_cells} cells "
      f"({len(freqs)} frequencies x {safe.sum()} time points), n = {len(SUBJECTS)} subjects")
print(f"    cells with p < .05 : {(p_unc < 0.05).sum()} ({100 * (p_unc < 0.05).mean():.2f} %)")
print(f"    expected by chance : {0.05 * n_cells:.0f} (5 % of {n_cells})")
print(f"    Bonferroni threshold for {n_cells} tests: p < {0.05 / n_cells:.3g}; cells surviving it: "
      f"{(p_unc < 0.05 / n_cells).sum()}")

rng = np.random.default_rng(20260918)
null_counts = []
for _ in range(200):
    flip = rng.choice([-1.0, 1.0], size=(len(SUBJECTS), 1, 1))
    null_counts.append(int((sstats.ttest_1samp(X_diff * flip, 0, axis=0).pvalue < 0.05).sum()))
null_counts = np.array(null_counts)
print(f"    the same count with the signs flipped at random (the null is true by construction), "
      f"200 draws: median {np.median(null_counts):.0f}, 95th percentile {np.percentile(null_counts, 95):.0f}, "
      f"max {null_counts.max()}")
print(f"    the largest connected blob of p < .05 cells under that null: "
      f"see the cluster null distribution in section 3 -- counting cells is not the right summary, "
      f"which is the point.")

fig, axes = plt.subplots(1, 2, figsize=(14, 4.0))
L4.plot_tfr(np.where(p_unc < 0.05, t_obs, np.nan), freqs, times[safe], ax=axes[0], vlim=(-6, 6),
            cbar_label="t (paired, across subjects)",
            title=f"Uncorrected p < .05 only ({(p_unc < 0.05).sum()} of {n_cells} cells)\n"
                  f"{len(SUBJECTS)} subjects -- this picture is not a result")
axes[1].hist(null_counts, bins=25, color="0.7")
axes[1].axvline((p_unc < 0.05).sum(), color="tab:red", lw=2,
                label=f"observed {(p_unc < 0.05).sum()} cells")
axes[1].axvline(0.05 * n_cells, color="tab:blue", lw=1.2, ls="--", label=f"5 % of cells = {0.05 * n_cells:.0f}")
axes[1].set(xlabel="Cells with uncorrected p < .05", ylabel="Sign-flip draws",
            title="How many cells pass uncorrected when there is nothing there\n(200 sign-flip draws)")
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

N_PERM = 1024          # 2**10 sign flips: exact for 10 subjects
adjacency = mne.stats.combine_adjacency(len(freqs), int(safe.sum()))
rows = []
results = {}
for thresh_p in (0.05, 0.01, 0.001):
    thr = sstats.t.ppf(1 - thresh_p / 2, len(SUBJECTS) - 1)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        T, clusters, cluster_p, H0 = mne.stats.permutation_cluster_1samp_test(
            X_diff, threshold=thr, n_permutations=N_PERM, tail=0, adjacency=adjacency,
            out_type="mask", seed=20260918, verbose=False)
    results[thresh_p] = (T, clusters, cluster_p, H0, thr)
    # Report the LARGEST cluster at each threshold whether or not it is significant: a null result with the
    # extent hidden is not a report, and the extent is what moves with the threshold.
    order = sorted(range(len(clusters)), key=lambda i: -abs(T[clusters[i]].sum()))
    for rank, i in enumerate(order[:2]):
        m = clusters[i]
        f_i, t_i = np.where(m)
        rows.append({"cluster-forming p": thresh_p, "|t| threshold": thr,
                     "rank by |sum t|": rank + 1, "p": cluster_p[i], "cells": int(m.sum()),
                     "freq (Hz)": f"{freqs[f_i.min()]:.0f}-{freqs[f_i.max()]:.0f}",
                     "time (s)": f"{times[safe][t_i.min()]:+.2f}..{times[safe][t_i.max()]:+.2f}",
                     "sum t": float(T[m].sum()),
                     "p < .05": bool(cluster_p[i] < 0.05)})
    if not len(clusters):
        rows.append({"cluster-forming p": thresh_p, "|t| threshold": thr, "rank by |sum t|": "-",
                     "p": float("nan"), "cells": 0, "freq (Hz)": "-", "time (s)": "-",
                     "sum t": float("nan"), "p < .05": False})
print(f"Cluster-permutation test, {N_PERM} sign flips (exact for {len(SUBJECTS)} subjects), two-tailed, "
      f"lattice adjacency over {len(freqs)} x {int(safe.sum())} cells")
print(L4.fmt_table(rows, floatfmt="{:.4f}"))

fig, axes = plt.subplots(1, 3, figsize=(17.5, 4.0))
for ax, thresh_p in zip(axes, (0.05, 0.01, 0.001)):
    T, clusters, cluster_p, H0, thr = results[thresh_p]
    L4.plot_tf_clusters(T, clusters, cluster_p, freqs, times[safe], ax=ax,
                        title=f"cluster-forming p < {thresh_p:g} (|t| > {thr:.2f})\n"
                              f"{sum(p < 0.05 for p in cluster_p)} significant cluster(s); smallest p = "
                              f"{min(cluster_p) if len(cluster_p) else float('nan'):.4f}")
fig.suptitle(f"The same data, three cluster-forming thresholds (t map, {len(SUBJECTS)} subjects, "
             f"left minus right imagery lateralisation)", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    T_tfce, cl_tfce, p_tfce, H0_tfce = mne.stats.permutation_cluster_1samp_test(
        X_diff, threshold=dict(start=0, step=0.2), n_permutations=N_PERM, tail=0, adjacency=adjacency,
        out_type="mask", seed=20260918, verbose=False)
p_map = p_tfce.reshape(T_tfce.shape) if p_tfce.ndim == 1 and p_tfce.size == T_tfce.size else None
print(f"Threshold-free cluster enhancement (start 0, step 0.2), {N_PERM} sign flips:")
if p_map is not None:
    print(f"    cells with corrected p < .05: {(p_map < 0.05).sum()} of {p_map.size} "
          f"({100 * (p_map < 0.05).mean():.2f} %); smallest corrected p = {p_map.min():.4f}")
else:
    print(f"    {sum(p < 0.05 for p in p_tfce)} of {len(p_tfce)} TFCE clusters at p < .05; smallest p = "
          f"{min(p_tfce):.4f}")
fig, ax = plt.subplots(figsize=(8.6, 4.2))
if p_map is not None:
    L4.plot_tfr(np.where(p_map < 0.05, T_tfce, np.nan), freqs, times[safe], ax=ax, vlim=(-6, 6),
                cbar_label="t (paired, across subjects)",
                title=f"TFCE, corrected p < .05 ({(p_map < 0.05).sum()} cells of {p_map.size})\n"
                      f"{len(SUBJECTS)} subjects, left minus right imagery lateralisation")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print("TFCE needs no cluster-forming threshold, which removes one analyst choice; it replaces it with the "
      "start and step of its integration, which are two more. What it buys is a corrected p per CELL rather "
      "than per cluster, so the map can be read pointwise -- but the inference is still 'the null is "
      "rejected somewhere', not 'the effect is here'.")

ROI_F = L4.MU_BAND
ROI_T = L4.TF_ACTIVE
fm = (freqs >= ROI_F[0]) & (freqs <= ROI_F[1])
tm = (times[safe] >= ROI_T[0]) & (times[safe] <= ROI_T[1])
roi = X_diff[:, fm][:, :, tm].mean(axis=(1, 2))
t_roi, p_roi = sstats.ttest_1samp(roi, 0)
print(f"PLAN A -- one a-priori ROI: mu {ROI_F[0]:g}-{ROI_F[1]:g} Hz over {ROI_T[0]:g}..{ROI_T[1]:g} s, "
      f"(C3 - C4) left minus right, {len(SUBJECTS)} subjects")
print(f"    mean {roi.mean():+.2f} percentage points, SD {roi.std(ddof=1):.2f}, "
      f"t({len(SUBJECTS) - 1}) = {t_roi:+.3f}, p = {p_roi:.4f}, dz = {roi.mean() / roi.std(ddof=1):+.3f}")
print(f"    tests run: 1. No correction needed, and nothing was chosen after seeing the data.")
T05, cl05, p05, H005, thr05 = results[0.05]
print()
print(f"PLAN B -- cluster-permutation over the whole map (4-40 Hz, {times[safe][0]:+.2f}..{times[safe][-1]:+.2f} s)")
print(f"    {len(cl05)} candidate clusters; {sum(p < 0.05 for p in p05)} at p < .05; smallest p = "
      f"{min(p05) if len(p05) else float('nan'):.4f}")
print(f"    tests run: 1 (the family-wise null). Analyst choices still open: the cluster-forming threshold, "
      f"the tail, the adjacency, the frequency range, the time range, the baseline, the normalisation and "
      f"the averaging order.")
print()
print(f"PLAN C -- uncorrected cell-by-cell tests, then describe whatever is significant")
print(f"    {(p_unc < 0.05).sum()} cells at p < .05, of which {0.05 * n_cells:.0f} are expected when nothing "
      f"is there. Analyst choices: all of Plan B's, plus which blob to describe and how to word it.")
print()
print(f"ANSWER -- ranked by degrees of freedom left to the analyst, fewest first: A, B, C.")
print(f"    A is not the most sensitive plan and it is not supposed to be. It is the one whose p-value means "
      f"what it says, because every choice it makes was made before the data.")

fig, ax = plt.subplots(figsize=(7.6, 4.2))
ax.bar(range(len(SUBJECTS)), roi, color=["tab:blue" if v > 0 else "tab:orange" for v in roi])
ax.axhline(0, color="k", lw=0.8)
ax.axhline(roi.mean(), color="tab:red", lw=1.4,
           label=f"mean {roi.mean():+.2f} pp (t = {t_roi:+.2f}, p = {p_roi:.4f})")
ax.set_xticks(range(len(SUBJECTS)))
ax.set_xticklabels(SUBJECTS, rotation=60, fontsize=7)
ax.set(xlabel="Subject", ylabel="(C3 - C4) left minus right (percentage points)",
       title=f"The a-priori ROI, one number per subject (percentage points)\n"
             f"mu {ROI_F[0]:g}-{ROI_F[1]:g} Hz, {ROI_T[0]:g}..{ROI_T[1]:g} s, {len(SUBJECTS)} subjects")
ax.legend(fontsize=8)
ax.grid(alpha=0.3, axis="y")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print(f"Subjects with the expected sign (positive = more C3-relative desynchronisation for RIGHT-hand "
      f"imagery, i.e. contralateral): {(roi > 0).sum()} of {len(SUBJECTS)}.")

try:
    print("nb-4-7-tf-cluster -- L4.7 numbers (draft; TODO(confirm) at author review)")
    print(f"Data: ds-eegbci, {len(SUBJECTS)} subjects ({SUBJECTS[0]}..{SUBJECTS[-1]}), runs "
          f"{'+'.join(f'R{r:02d}' for r in L4.MI_RUNS)} (PhysioNet DOI "
          f"{L4.DATASETS_L4['ds-eegbci']['dataset_doi']}; ODC-By 1.0). "
          f"helpers_l4.MI_PIPELINE unchanged; percent change from "
          f"{L4.TF_BASELINE[0]:g}..{L4.TF_BASELINE[1]:g} s, band-first order.")
    print(f"Contrast: (C3 - C4) lateralisation, left-hand imagery minus right-hand imagery, "
          f"{len(freqs)} frequencies x {int(safe.sum())} edge-safe time points = {n_cells} cells.")
    print()
    print("L4.7's exercise is a RANK exercise (three analysis plans by degrees of freedom left to the "
          "analyst), so there is no numeric key. What this notebook supplies:")
    print("    ANSWER -- fewest degrees of freedom first: (A) a-priori ROI, (B) cluster-permutation over "
          "the whole map, (C) uncorrected cell-by-cell tests described after the fact.")
    print(f"    A: 1 test, 0 post-hoc choices. t({len(SUBJECTS) - 1}) = {t_roi:+.3f}, p = {p_roi:.4f}, "
          f"mean {roi.mean():+.2f} percentage points, dz {roi.mean() / roi.std(ddof=1):+.3f}.")
    print(f"    B: 1 family-wise test over {n_cells} cells; {sum(p < 0.05 for p in p05)} cluster(s) at "
          f"p < .05; at least 8 choices remain open (cluster-forming threshold, tail, adjacency, frequency "
          f"range, time range, baseline window, normalisation, averaging order).")
    print(f"    C: {(p_unc < 0.05).sum()} cells at p < .05 where {0.05 * n_cells:.0f} are expected under the "
          f"null; sign-flipped draws give a median of {np.median(null_counts):.0f} and as many as "
          f"{null_counts.max()} with nothing there at all.")
    print()
    print("Supporting -- the cluster test, three cluster-forming thresholds, LARGEST cluster of each "
          "(significant or not):")
    for r in rows:
        print(f"    forming p {r['cluster-forming p']:<6g} (|t| > {r['|t| threshold']:.2f}) rank "
              f"{r['rank by |sum t|']}: p = {r['p']:.4f}{'  (p < .05)' if r['p < .05'] else ''}, "
              f"{r['cells']:5d} cells, {r['freq (Hz)']:>7s} Hz, {r['time (s)']:>15s} s, "
              f"sum t = {r['sum t']:+.1f}")
    n_sig_any = sum(sum(p < 0.05 for p in results[tp][2]) for tp in results)
    print(f"    RESULT: {'no cluster reaches p < .05 at any of the three thresholds' if n_sig_any == 0 else f'{n_sig_any} cluster(s) reach p < .05'}. "
          f"That is the honest report for this contrast on {len(SUBJECTS)} subjects, and the lesson does not "
          f"depend on the result being positive.")
    print(f"    The cluster-forming threshold moves the smallest p from "
          + ", ".join(f"{min(results[tp][2]) if len(results[tp][2]) else float('nan'):.4f} at forming p {tp:g}"
                      for tp in (0.05, 0.01, 0.001))
          + " and the largest cluster's extent from "
          + " to ".join(f"{r['cells']} cells" for r in rows if r["rank by |sum t|"] == 1)
          + ", on the same data. A threshold chosen after seeing which one works is "
          "pf-cluster-inference-misread in its most direct form.")
    if p_map is not None:
        print(f"    TFCE: {(p_map < 0.05).sum()} of {p_map.size} cells at corrected p < .05, smallest "
              f"corrected p = {p_map.min():.4f}.")
    print()
    print(f"Supporting -- Bonferroni over {n_cells} cells would need p < {0.05 / n_cells:.3g}; "
          f"{(p_unc < 0.05 / n_cells).sum()} cells reach it. That is the price of not using the map's "
          f"structure, and it is why cluster methods exist.")
    print()
    print(f"Supporting -- per-subject ROI values (percentage points): "
          + ", ".join(f"{s} {v:+.1f}" for s, v in zip(SUBJECTS, roi)))
    print(f"    {(roi > 0).sum()} of {len(SUBJECTS)} subjects show the contralateral sign.")
    print()
    print("Pitfalls: pf-uncorrected-timepoint-tests, pf-cluster-inference-misread. "
          "Widget: w-cluster-permutation-viz (mode tf). Statistics thread: L3.7 -> L4.7 -> L6.1 -> L6.4.")
    print("TODO(confirm): the catalog carries no published lateralisation effect size for this dataset, so "
          "nothing above is compared with a literature value.")
finally:
    dl.finish()