1"""Dissipative softmax latent layer: finite-state CTMC verification."""
 2from __future__ import annotations
 3import numpy as np
 4
 5
 6def softmax(x):
 7    x = np.asarray(x, dtype=float)
 8    y = x - np.max(x)
 9    e = np.exp(y)
10    return e / e.sum()
11
12
13def conditional_rates(X, k=1.0):
14    """Reversible conditional rates with log(W_ij/W_ji)=X_j-X_i."""
15    X = np.asarray(X, dtype=float)
16    d = X[None, :] - X[:, None]
17    W = k * np.exp(d / 2)
18    np.fill_diagonal(W, 0.0)
19    return W
20
21
22def generator(X, r, k=1.0, drive=0.15):
23    """Row-generator on state 0 plus K conditional states.
24
25    Reset exchange is  i <-> 0 with rates r and r*p_i.  Thus, without
26    the driven cycle, the stationary law is P0=1/2 and P_i=p_i/2 for every r.
27    The one-way cycle carries dissipation, but its occupation perturbation is
28    O(drive/r) as r/k grows, which is the sector-separation prediction.
29    """
30    X = np.asarray(X, dtype=float)
31    p = softmax(X); K = len(X)
32    Q = np.zeros((K + 1, K + 1))
33    W = conditional_rates(X, k)
34    for i in range(K):
35        for j in range(K):
36            if i != j: Q[i+1, j+1] += W[i, j]
37        Q[i+1, 0] += r
38        Q[0, i+1] += r * p[i]
39    # Directed reset cycle 0 -> 1 -> ... -> K -> 0.
40    Q[0, 1] += drive
41    for i in range(1, K): Q[i, i+1] += drive
42    Q[K, 0] += drive
43    np.fill_diagonal(Q, -Q.sum(axis=1))
44    return Q, p
45
46
47def stationary(Q):
48    A = Q.T.copy(); A[-1] = 1.0
49    b = np.zeros(Q.shape[0]); b[-1] = 1.0
50    pi = np.linalg.solve(A, b)
51    return np.maximum(pi, 0) / np.maximum(pi, 0).sum()
52
53
54def cycle_current(pi, Q):
55    edges = [(0, 1)] + [(i, i+1) for i in range(1, Q.shape[0]-1)] + [(Q.shape[0]-1, 0)]
56    return float(np.mean([pi[i]*Q[i,j] - pi[j]*Q[j,i] for i,j in edges]))
57
58
59def rate_ratio_error(X, W):
60    X = np.asarray(X); mask = ~np.eye(len(X), dtype=bool)
61    return float(np.max(np.abs(np.log(W[mask]/W.T[mask]) - (X[None,:]-X[:,None])[mask])))
62
63
64def run_sweep():
65    X = np.array([1.0, 0.2, -0.5, -1.0]); k = 1.0; drive = 0.15
66    rows = []
67    for ratio in [0.1, 0.3, 1, 3, 10, 30, 100, 300]:
68        Q, p = generator(X, ratio*k, k, drive); pi = stationary(Q)
69        cond = pi[1:] / pi[1:].sum()
70        kl = float(np.sum(cond*np.log(np.maximum(cond,1e-300)/p)))
71        rows.append((ratio, kl, float(np.max(np.abs(cond-p))), pi[0], cycle_current(pi,Q)))
72    return X, rows, rate_ratio_error(X, conditional_rates(X,k))
73
74if __name__ == '__main__':
75    X, rows, ratio_err = run_sweep()
76    print('X=', X)
77    print('max conditional log-rate-ratio error=', ratio_err)
78    print('r_over_k, KL(cond||p), max_abs_error, P0, cycle_current')
79    for row in rows: print('%.3g %.8g %.8g %.8g %.8g' % row)
80    tail = np.array(rows[-4:])
81    print('asymptotic slope KL vs r/k (predicted -2):', np.polyfit(np.log(tail[:,0]), np.log(tail[:,1]), 1)[0])
82    print('asymptotic slope error vs r/k (predicted -1):', np.polyfit(np.log(tail[:,0]), np.log(tail[:,2]), 1)[0])