File size: 3,383 Bytes
31d1d63 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 | import numpy as np, collections
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
Ac = np.array([-200, -100, -170, 15.0])
a = np.array([-1, -1, -6.5, 0.7])
b = np.array([0, 0, 11, 0.6])
c = np.array([-10, -10, -6.5, 0.7])
X = np.array([1, 0, -0.5, -1.0])
Y = np.array([0, 0.5, 1.5, 1.0])
def E(p):
x, y = p
return float(
(Ac * np.exp(a * (x - X) ** 2 + b * (x - X) * (y - Y) + c * (y - Y) ** 2)).sum()
)
LAG, N, TL = MP.LAG_OPT, MP.N_OPT, 20
sad = gt.saddle_points()
S1 = sad[np.argmin(sad[:, 0])]
S2 = sad[np.argmax(sad[:, 0])]
trajs = V.load_trajs(str(V.DEFAULT_DATA))
AGG = collections.defaultdict(list)
DIV = collections.defaultdict(list)
for seed in [0, 1, 2, 3, 4]:
model = mb.build_full_msm(
trajs, trajs, lag=LAG, tica_lag=TL, k=V.K, n_comp=2, seed=seed, 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)
flr = np.percentile(rho[react], 25)
n1 = np.linalg.norm(xy - S1, axis=1) < 0.25
qs2 = q[np.argmin(np.linalg.norm(xy - S2, axis=1))]
for al in [0.0, 0.5, 1.0, 1.5, 2.0, 3.0]:
w = np.ones_like(rho) if al == 0 else np.exp(-LAG * al * np.maximum(rho, flr))
w = w.copy()
w[~on] = 0.0
ps = bridge.sample_bridge_paths(
bridge.doob_kernels(model.P, end, N, w=w), start, n_paths=3000, seed=0
)
rec = np.mean([any(n1[s] for s in p) for p in ps]) * 100
DIV[al].append((len({s for p in ps for s in p}), len({tuple(p) for p in ps})))
r = MP.channel(ps, q, xy, on, thru, start, end)
qr = [q[np.argmin(np.linalg.norm(xy - p, axis=1))] for p in r]
cr = [i for i in range(len(qr) - 1) if (qr[i] < qs2) != (qr[i + 1] < qs2)]
d2 = (E((r[cr[0]] + r[cr[0] + 1]) / 2) - E(S2)) / 16.7 if cr else np.nan
AGG[al].append((rec, d2))
print(f"TICA lag {TL}, count lag {LAG}, N={N}; 5 clustering seeds; kT=16.7")
print(f"{'alpha':>6} {'S1 recov':>16} {'S2 dE (kT)':>16} {'states':>10} {'routes':>10}")
for al, v in sorted(AGG.items()):
rc = np.array([x[0] for x in v])
d2 = np.array([x[1] for x in v])
st = np.array([x[0] for x in DIV[al]])
rt = np.array([x[1] for x in DIV[al]])
print(
f"{al:6.1f} {rc.mean():8.0f}%+-{rc.std():3.0f} {np.nanmean(d2):11.1f}+-{np.nanstd(d2):3.1f} "
f"{st.mean():6.0f}+-{st.std():2.0f} {rt.mean():6.0f}+-{rt.std():2.0f}"
)
|