Buckets:
| #!/usr/bin/env python3 | |
| """build_l2_decay — 给 L2 补模型无关的衰减尺度（不拟合任何曲线族）。 | |
| 两个量，分别作用在两条不同的曲线上，不能混用： | |
| rho_half_ms / rho_zero_ms | |
| 作用在"跨单元共享相关曲线" rho 上（rho_curve 的输出）。 | |
| 定义：从 lag=1 的值出发，首次降到 50% / 0 的滞后（线性插值），单位 ms。 | |
| 尖峰的合适不变量；对 EEG 无效（通道间已被白化）。 | |
| acf1_half_ms / acf1_zero_ms | |
| 作用在"平均单通道自相关"上（rho_curve 同一次 FFT 的对角项 dia 归一化）。 | |
| 定义同上。不受跨通道白化影响 —— 白化只改通道间关系，不改单通道时间自相关。 | |
| 两者都只要求曲线在 max_lag 内降到阈值；降不到就标 censored（记录下界）， | |
| 不外推、不假设函数族。 | |
| 数值对拍：acf1 用 40 个单元的逐条 ACF 平均验证（应与 FFT 路径一致）。 | |
| 输出: /domains/framework/l2_decay.parquet，并左连接回 l2_invariants.parquet | |
| """ | |
| from __future__ import annotations | |
| import glob | |
| import sys | |
| import time | |
| from pathlib import Path | |
| import numpy as np | |
| import pandas as pd | |
| ROOT = Path("/domains") | |
| OUT = ROOT / "framework" | |
| MAX_LAG_SPIKE = 400 | |
| MAX_LAG_EEG = 100 | |
| N_EEG_TRIALS = 100 | |
| N_4D_TRIALS = 50 | |
| FRAC = 0.5 | |
| def n2(n): | |
| k = 1 | |
| while k < n: | |
| k <<= 1 | |
| return k | |
| def curve_parts(x, max_lag): | |
| """与 rho_curve 相同的 FFT 路径；额外返回归一化对角项（平均单通道 ACF）。""" | |
| T, d = x.shape | |
| xc = x - x.mean(axis=0, keepdims=True) | |
| N = n2(2 * T) | |
| F = np.fft.rfft(xc, n=N, axis=0) | |
| Fs = F.sum(axis=1) | |
| tot = np.fft.irfft(Fs * np.conj(Fs), n=N)[:max_lag + 1] | |
| dia = np.fft.irfft((np.abs(F) ** 2).sum(axis=1), n=N)[:max_lag + 1] | |
| nrm = np.arange(T, T - max_lag - 1, -1) | |
| tot = tot / np.maximum(nrm, 1) | |
| dia = dia / np.maximum(nrm, 1) | |
| t0 = max(tot[0], 1e-12) | |
| rho = (tot - dia) / t0 | |
| acf1 = dia / max(dia[0], 1e-12) | |
| return rho, acf1 | |
| def cross_time(curve, bin_ms, frac): | |
| """首次降到 frac*curve[1] 的滞后（线性插值）。frac=0 → 过零。""" | |
| if len(curve) < 3: | |
| return None, "too_short" | |
| y1 = float(curve[1]) | |
| if not np.isfinite(y1) or y1 <= 0: | |
| return None, "bad_start" | |
| thr = frac * y1 | |
| idx = np.flatnonzero(curve[2:] <= thr) | |
| if len(idx) == 0: | |
| return None, "censored" | |
| l = int(idx[0]) + 2 | |
| y0, yb = float(curve[l - 1]), float(curve[l]) | |
| f = (y0 - thr) / (y0 - yb) if y0 > yb else 0.0 | |
| return float((l - 1 + f) * bin_ms), "ok" | |
| rows = [] | |
| t00 = time.time() | |
| # ---------- 数值对拍 ---------- | |
| f0 = sorted(glob.glob(str(ROOT / "steinmetz/spikes_session_*.npz")))[0] | |
| z0 = np.load(f0, mmap_mode="r") | |
| x0 = np.round(np.asarray(z0["spikes"], dtype=np.float64) * 8).T | |
| ml0 = min(MAX_LAG_SPIKE, x0.shape[0] - 2) | |
| rho_f, acf_f = curve_parts(x0.astype(np.float32), ml0) | |
| sub = np.random.default_rng(0).choice(x0.shape[1], size=min(40, x0.shape[1]), replace=False) | |
| acc, rpos = [], [] | |
| for j in sub: | |
| v = x0[:, j] - x0[:, j].mean() | |
| ac = np.correlate(v, v, mode="full")[len(v) - 1:][:ml0 + 1] | |
| acc.append(ac / max(ac[0], 1e-12)) | |
| print(f"[check] acf1 max|fast-brute| = {np.abs(acf_f - np.mean(acc, axis=0)).max():.3e}", | |
| flush=True) | |
| del x0, z0 | |
| def push(member, slice_id, rho, acf1, bin_ms, curv_kind, note=""): | |
| """curv_kind: 'spikes'（rho 与 acf1 都用）/ 'eeg'（只用 acf1）。""" | |
| rec = dict(member=member, slice=slice_id, note=note, | |
| rho_half_ms=None, rho_zero_ms=None, | |
| rho_half_status=None, rho_zero_status=None, | |
| acf1_lag1=float(acf1[1]) if len(acf1) > 1 else None, | |
| acf1_half_ms=None, acf1_zero_ms=None, | |
| acf1_half_status=None, acf1_zero_status=None, | |
| max_lag_ms=len(acf1) * bin_ms) | |
| rec["acf1_half_ms"], rec["acf1_half_status"] = cross_time(acf1, bin_ms, FRAC) | |
| rec["acf1_zero_ms"], rec["acf1_zero_status"] = cross_time(acf1, bin_ms, 0.0) | |
| if curv_kind == "spikes" and rho is not None: | |
| rec["rho_half_ms"], rec["rho_half_status"] = cross_time(rho, bin_ms, FRAC) | |
| rec["rho_zero_ms"], rec["rho_zero_status"] = cross_time(rho, bin_ms, 0.0) | |
| rows.append(rec) | |
| # ---------- 1. 原生尖峰 ---------- | |
| for member, pat, bin_ms in [("steinmetz", "steinmetz/spikes_session_*.npz", 50.0), | |
| ("primate_macaque_nlb", "primate/macaque_nlb/**/*.npz", 20.0)]: | |
| for f in sorted(glob.glob(str(ROOT / pat), recursive=True)): | |
| z = np.load(f, mmap_mode="r") | |
| k = "spikes" if "spikes" in z.keys() else list(z.keys())[0] | |
| a = np.asarray(z[k], dtype=np.float32) | |
| x = a.T if a.shape[0] < a.shape[1] else a | |
| if float((np.abs(x.astype(np.float64) * 8 | |
| - np.round(x.astype(np.float64) * 8)) < 1e-5).mean()) > 0.999: | |
| x = np.round(x.astype(np.float64) * 8).astype(np.float32) | |
| ml = min(MAX_LAG_SPIKE, x.shape[0] - 2) | |
| rho, acf1 = curve_parts(x, ml) | |
| push(member, Path(f).stem, rho, acf1, bin_ms, "spikes") | |
| print(f" {member} {Path(f).stem} [{time.time()-t00:.0f}s]", flush=True) | |
| del a, x | |
| # ---------- 2. 人类 EEG：逐 trial 的通道 ACF 平均 ---------- | |
| for f in sorted(glob.glob(str(ROOT / "primate/human_eeg/**/*.npz"), recursive=True)): | |
| if not f.endswith(".npz"): | |
| continue | |
| z = np.load(f, mmap_mode="r") | |
| k = "eeg" if "eeg" in z.keys() else list(z.keys())[0] | |
| a = np.asarray(z[k], dtype=np.float32) | |
| nt, nch, T = a.shape | |
| ml = min(MAX_LAG_EEG, T - 2) | |
| sel = np.random.default_rng(0).choice(nt, size=min(N_EEG_TRIALS, nt), replace=False) | |
| acc = np.mean([curve_parts(a[i].T, ml)[1] for i in sel], axis=0) | |
| push("primate_human_eeg", Path(f).stem, None, acc, 4.0, "eeg", | |
| note=f"single_channel_acf_over_{len(sel)}_trials") | |
| print(f" human_eeg {Path(f).stem} [{time.time()-t00:.0f}s]", flush=True) | |
| del a | |
| # ---------- 3. 四轴 EEG ---------- | |
| for member, pat in [("omnibrain4d_eeg_test", "omnibrain4d/eeg/**/*.npz"), | |
| ("omnibrain4d_eeg_train", "omnibrain4d/eeg_train/**/*.npz")]: | |
| for f in sorted(glob.glob(str(ROOT / pat), recursive=True)): | |
| z = np.load(f, mmap_mode="r") | |
| if "eeg" not in z.keys(): | |
| continue | |
| a = np.asarray(z["eeg"], dtype=np.float32) | |
| if a.ndim != 4: | |
| continue | |
| a = a[:50] | |
| ni, nr, nch, T = a.shape | |
| ml = min(MAX_LAG_EEG, T - 2) | |
| flat = a.transpose(0, 1, 3, 2).reshape(ni * nr, T, nch) | |
| sel = np.random.default_rng(0).choice(len(flat), size=min(N_4D_TRIALS, len(flat)), | |
| replace=False) | |
| acc = np.mean([curve_parts(flat[i], ml)[1] for i in sel], axis=0) | |
| push(member, Path(f).stem, None, acc, 4.0, "eeg", | |
| note=f"single_channel_acf_over_{len(sel)}_trials") | |
| print(f" {member} {Path(f).stem} {a.shape} [{time.time()-t00:.0f}s]", flush=True) | |
| del a, flat | |
| tdf = pd.DataFrame(rows) | |
| tdf.to_parquet(OUT / "l2_decay.parquet", index=False) | |
| P = OUT / "l2_invariants.parquet" | |
| df = pd.read_parquet(P) | |
| newcols = [c for c in tdf.columns if c not in ("member", "slice", "note")] | |
| for c in newcols: | |
| if c in df.columns: | |
| df = df.drop(columns=[c]) | |
| df = df.merge(tdf[["member", "slice"] + newcols], on=["member", "slice"], how="left") | |
| df.to_parquet(P, index=False) | |
| df.to_csv(OUT / "l2_invariants.csv", index=False) | |
| print(f"\nrows={len(tdf)} [{time.time()-t00:.0f}s]") | |
| for c in ["rho_half_status", "rho_zero_status", "acf1_half_status", "acf1_zero_status"]: | |
| print(f"\n{c}: {tdf[c].value_counts().to_dict()}") | |
| print("\nrho_half_ms (ok):") | |
| print(tdf[tdf.rho_half_status == "ok"].groupby("member").rho_half_ms | |
| .agg(["size", "median", lambda s: s.quantile(.1), | |
| lambda s: s.quantile(.9)]).to_string()) | |
| print("\nacf1_half_ms (ok):") | |
| print(tdf[tdf.acf1_half_status == "ok"].groupby("member").acf1_half_ms | |
| .agg(["size", "median", lambda s: s.quantile(.1), | |
| lambda s: s.quantile(.9)]).to_string()) | |
| print("\nacf1_lag1 median:") | |
| print(tdf.groupby("member").acf1_lag1.median().to_string()) | |
Xet Storage Details
- Size:
- 8.19 kB
- Xet hash:
- 5be6f2116dbbb46eee19e80ecc26307e80132d9f84e5a939fb30591b15d37011
·
Xet efficiently stores files, intelligently splitting them into unique chunks and accelerating uploads and downloads. More info.