| import numpy as np, json |
| from scoobdoob.validation import msm_build as mb, tpt, bridge |
| from scoobdoob.validation import mb_ground_truth as gt |
| from scoobdoob.validation import mb_data as V |
| import make_mb_panels as MP |
|
|
| LAG, N, TL, ALPHA = MP.LAG_OPT, MP.N_OPT, MP.TICA_LAG, 1.5 |
| CLUST_SEEDS = list(range(20)) |
| PATH_SEEDS = [0, 1, 2] |
| NPATH = 2000 |
| sad = gt.saddle_points() |
| S1 = sad[np.argmin(sad[:, 0])] |
| trajs = V.load_trajs(str(V.DEFAULT_DATA)) |
|
|
| raw_m, fk_m, within = [], [], [] |
| for cs in CLUST_SEEDS: |
| model = mb.build_full_msm( |
| trajs, trajs, lag=LAG, tica_lag=TL, k=V.K, n_comp=2, seed=cs, use_tica=True |
| ) |
| A, B, _ = mb.q2_sign_endpoints(model) |
| start, end = V.orient_endpoints(model, A, B) |
| xy = model.state_xy |
| _, q, _ = tpt.committor_agreement(model, start, end) |
| rf = tpt.reactive_flux(model, start, end) |
| thru = np.asarray(rf.gross_flux).sum(1) |
| share = thru / thru.sum() |
| on = share > 0.5 / model.n_active |
| on[list(start)] = True |
| on[list(end)] = True |
| Axy = model.active_to_xy(start).mean(0) |
| Bxy = model.active_to_xy(end).mean(0) |
| land = np.vstack([Axy, Bxy, sad]) |
| lo = land.min(0) - 0.35 |
| hi = land.max(0) + 0.35 |
| on = on & ( |
| (xy[:, 0] >= lo[0]) |
| & (xy[:, 0] <= hi[0]) |
| & (xy[:, 1] >= lo[1]) |
| & (xy[:, 1] <= hi[1]) |
| ) |
| labels = np.concatenate(model.dtrajs) |
| counts = np.bincount(labels, minlength=model.n_active)[: model.n_active].astype( |
| float |
| ) |
| rho = np.maximum(counts, 1) / counts.mean() |
| react = (q > 0.05) & (q < 0.95) |
| wfk = np.exp(-LAG * ALPHA * np.maximum(rho, np.percentile(rho[react], 25))) |
| wfk[~on] = 0.0 |
| wraw = np.ones_like(wfk) |
| wraw[~on] = 0.0 |
| n1 = np.linalg.norm(xy - S1, axis=1) < 0.25 |
| Kraw = bridge.doob_kernels(model.P, end, N, w=wraw) |
| Kfk = bridge.doob_kernels(model.P, end, N, w=wfk) |
| r_ps, f_ps = [], [] |
| for psd in PATH_SEEDS: |
| pr = bridge.sample_bridge_paths(Kraw, start, n_paths=NPATH, seed=psd) |
| pf = bridge.sample_bridge_paths(Kfk, start, n_paths=NPATH, seed=psd) |
| r_ps.append(np.mean([any(n1[s] for s in p) for p in pr]) * 100) |
| f_ps.append(np.mean([any(n1[s] for s in p) for p in pf]) * 100) |
| raw_m.append(np.mean(r_ps)) |
| fk_m.append(np.mean(f_ps)) |
| within.append((np.std(r_ps), np.std(f_ps))) |
| print( |
| f" [clust seed {cs:2d}] raw {raw_m[-1]:5.1f}% fk {fk_m[-1]:5.1f}% " |
| f"diff {fk_m[-1] - raw_m[-1]:+5.1f} (path-seed sd {within[-1][0]:.1f}/{within[-1][1]:.1f})" |
| ) |
|
|
| raw = np.array(raw_m) |
| fk = np.array(fk_m) |
| d = fk - raw |
| n = len(d) |
| print( |
| f"\n--- variance decomposition (n={n} clustering seeds, {len(PATH_SEEDS)} path seeds, " |
| f"{NPATH} paths each) ---" |
| ) |
| wr = np.mean([w[0] for w in within]) |
| wf = np.mean([w[1] for w in within]) |
| print(f"between-clustering sd : raw {raw.std(ddof=1):.1f} fk {fk.std(ddof=1):.1f}") |
| print(f"within (path-sampling) sd : raw {wr:.1f} fk {wf:.1f}") |
| print( |
| f"-> clustering dominates by {fk.std(ddof=1) / max(wf, 1e-9):.0f}x; Monte-Carlo error is negligible" |
| ) |
|
|
| print(f"\n--- paired comparison, FK vs raw ---") |
| print(f"raw {raw.mean():.1f} +- {raw.std(ddof=1):.1f}") |
| print(f"fk {fk.mean():.1f} +- {fk.std(ddof=1):.1f}") |
| se = d.std(ddof=1) / np.sqrt(n) |
| print(f"paired diff {d.mean():+.1f} (sd {d.std(ddof=1):.1f}, se {se:.1f})") |
| print( |
| f" min diff over seeds {d.min():+.1f}, max {d.max():+.1f} " |
| f"(sign consistent: {'yes' if (d > 0).all() else 'NO'})" |
| ) |
| print(f"Cohen's d_z = {d.mean() / d.std(ddof=1):.2f}") |
| try: |
| from scipy import stats |
|
|
| t, pt = stats.ttest_rel(fk, raw) |
| w, pw = stats.wilcoxon(fk, raw) |
| ci = stats.t.interval(0.95, n - 1, loc=d.mean(), scale=se) |
| print(f"paired t: t={t:.1f}, p={pt:.2e}") |
| print( |
| f"Wilcoxon signed-rank: W={w:.0f}, p={pw:.2e} " |
| f"(floor for n={n} is {2 / 2**n:.1e})" |
| ) |
| print(f"95% CI of the mean difference: [{ci[0]:.1f}, {ci[1]:.1f}]") |
| except ImportError: |
| print("scipy unavailable; bootstrap only") |
| rng = np.random.default_rng(0) |
| bs = np.array([rng.choice(d, n, replace=True).mean() for _ in range(20000)]) |
| print( |
| f"bootstrap 95% CI of mean diff: [{np.percentile(bs, 2.5):.1f}, {np.percentile(bs, 97.5):.1f}]" |
| ) |
| json.dump( |
| { |
| "raw": raw.tolist(), |
| "fk": fk.tolist(), |
| "clust_seeds": CLUST_SEEDS, |
| "path_seeds": PATH_SEEDS, |
| "n_paths": NPATH, |
| "tica_lag": TL, |
| "lag": LAG, |
| "N": N, |
| "alpha": ALPHA, |
| }, |
| open("data/mb_stats.json", "w"), |
| indent=2, |
| ) |
| print("[done]") |
|
|
