1import numpy as np
 2
 3SEED = 1149
 4
 5def stationary(T):
 6    w, v = np.linalg.eig(T.T)
 7    x = np.real(v[:, np.argmin(abs(w - 1))])
 8    if x.sum() < 0: x = -x
 9    x = np.maximum(x, 0); return x / x.sum()
10
11def make_chain():
12    # Eight states grouped into four adjacent coarse states.
13    pc = np.array([18., 12., 6., 1.]) / 37.
14    pi = np.repeat(pc / 2, 2)
15    T = np.zeros((8, 8))
16    for i in range(7):
17        # Equalize neighboring probability flows, making pi stationary.
18        T[i, i+1] = min(.22, .12 * min(1., pi[i+1]/pi[i]))
19        T[i+1, i] = min(.22, .12 * min(1., pi[i]/pi[i+1]))
20    T[np.diag_indices(8)] = 1 - T.sum(1)
21    Q = T.copy(); Q[6:8] = [.5, .5, 0, 0, 0, 0, 0, 0]
22    return T, Q, pi
23
24def coarse(T, pi, alpha=0.):
25    groups = [np.array([2*i, 2*i+1]) for i in range(4)]
26    den = np.array([pi[g].sum() for g in groups])
27    out = np.zeros((4, 4))
28    for i, gi in enumerate(groups):
29        for j, gj in enumerate(groups):
30            out[i,j] = (pi[gi,None] * T[np.ix_(gi,gj)]).sum() / max(den[i], 1e-15)
31    return out, den
32
33def hit_time(T, start_dist, sink=(6,7)):
34    keep = [i for i in range(8) if i not in sink]
35    h = np.linalg.solve(np.eye(len(keep))-T[np.ix_(keep,keep)], np.ones(len(keep)))
36    return sum(start_dist[i] * h[keep.index(i)] for i in keep if start_dist[i] > 0)
37
38def sample_counts(T, n, rg, alpha=.5):
39    s = stationary(T); starts = rg.choice(8, n, p=s)
40    ends = np.array([rg.choice(8, p=T[i]) for i in starts])
41    C = np.full((8,8), alpha); np.add.at(C, (starts, ends), 1)
42    return C / C.sum(1, keepdims=True)
43
44def main():
45    P, Q, pi = make_chain(); piq = stationary(Q)
46    _, pic = coarse(P, pi); _, piqc = coarse(Q, piq)
47    reset = np.array([.5,.5,0,0,0,0,0,0])
48    direct = hit_time(P, reset)
49    hill = 1/piqc[3] - 1
50    print('CORE_CHECK')
51    print('P_stationarity_error %.3e' % max(abs(pi@P-pi)))
52    print('Q_stationarity_error %.3e' % max(abs(piq@Q-piq)))
53    print('equilibrium_target', np.round(pic,6))
54    print('ness_target', np.round(piqc,6))
55    print('reset_process_MFPT direct %.6f hill %.6f abs_error %.3e' % (direct,hill,abs(direct-hill)))
56
57    print('PREDICTION_1 lag-invariant matched stationary error: lag, equilibrium, NESS')
58    for lag in [1,2,4,8,16,32]:
59        ep = max(abs(stationary(coarse(np.linalg.matrix_power(P,lag),pi)[0])-pic))
60        eq = max(abs(stationary(coarse(np.linalg.matrix_power(Q,lag),piq)[0])-piqc))
61        print(lag, '%.3e %.3e' % (ep,eq))
62
63    print('PREDICTION_2 reset-strength: mix, dual_MFPT_error, single_MFPT_error, sink_rate')
64    for mix in [0,.25,.5,.75,1.]:
65        reset = np.array([1-mix,mix,0,0,0,0,0,0])
66        q = P.copy(); q[6:8] = reset
67        sq = stationary(q); _, sc = coarse(q,sq)
68        truth = hit_time(P,reset); dual = 1/sc[3]-1
69        single = 1/pic[3]-1
70        print('%.2f %.6f %.6f %.6f' % (mix,abs(dual-truth),abs(single-truth),sc[3]))
71
72    print('PREDICTION_3 sample-size RMS NESS occupancy error')
73    ns=[400,1600,6400,25600]; vals=[]
74    for n in ns:
75        de=[]; se=[]
76        for r in range(30):
77            rg=np.random.default_rng(SEED+100000+n+r)
78            pe=sample_counts(P,n,rg); qe=sample_counts(Q,n,rg)
79            spe=stationary(pe); sqe=stationary(qe)
80            _, pce=coarse(pe,spe); _, qce=coarse(qe,sqe)
81            de.append(np.linalg.norm(qce-piqc)); se.append(np.linalg.norm(pce-piqc))
82        d=float(np.sqrt(np.mean(np.square(de)))); s=float(np.sqrt(np.mean(np.square(se))))
83        vals.append(d); print(n,'dual %.6f single %.6f'%(d,s))
84    print('dual_loglog_slope %.3f predicted -0.5' % np.polyfit(np.log(ns),np.log(vals),1)[0])
85
86if __name__ == '__main__': main()