1import json, math, random
  2import numpy as np
  3from sklearn.datasets import load_digits
  4from sklearn.linear_model import LogisticRegression
  5from sklearn.model_selection import train_test_split
  6from sklearn.metrics import accuracy_score
  7
  8SEED = 1248
  9np.random.seed(SEED)
 10random.seed(SEED)
 11
 12
 13def sigmoid(z):
 14    return 1.0 / (1.0 + np.exp(-np.clip(z, -40, 40)))
 15
 16
 17def math_checks():
 18    out = {}
 19    # Prediction 1: around target rho, homeostatic scalar update contracts with
 20    # q = 1 - eta*lambda*[rho(1-rho)]^2.
 21    rho, eta, z0, tol = .7, .8, -2., 1e-3
 22    alpha_r = rho * (1-rho)
 23    rows = []
 24    for lam in [.05, .1, .2, .4, .8]:
 25        z = z0
 26        q = 1 - eta * lam * alpha_r**2
 27        pred = math.log(tol / abs(sigmoid(z0)-rho)) / math.log(abs(q))
 28        observed = None
 29        for k in range(200000):
 30            h = sigmoid(z)
 31            z += eta * lam * h * (1-h) * (rho-h)
 32            if abs(sigmoid(z)-rho) < tol:
 33                observed = k + 1
 34                break
 35        rows.append({'lambda': lam, 'observed_steps': observed,
 36                     'predicted_linear_steps': pred, 'predicted_q': q})
 37    out['homeostatic_contraction'] = rows
 38
 39    # Prediction 2: with scalar A=[a], C_A=a^2, so temporal drive scales as a^2.
 40    dh = .23
 41    rows = []
 42    for a in [.25, .5, 1., 2., 4.]:
 43        rows.append({'a': a, 'observed_ratio': (a*a*dh)/dh,
 44                     'predicted_ratio': a*a})
 45    out['projection_scaling'] = rows
 46
 47    # Prediction 3: alpha is maximal at h=.5 and goes to zero at saturation.
 48    hs = [.01, .05, .1, .25, .5, .75, .9, .95, .99]
 49    out['gain_suppression'] = [
 50        {'h': h, 'alpha': h*(1-h), 'predicted_relative': 4*h*(1-h)}
 51        for h in hs]
 52    return out
 53
 54
 55def local_train(X, C=24, M=12, epochs=3, eta=.22, etab=.01,
 56                lam=.12, rho=.5, ordered=True, seed=0):
 57    rng = np.random.RandomState(seed)
 58    n = X.shape[1]
 59    W = rng.normal(0, .18, (C, n))
 60    b = np.zeros(C)
 61    A = rng.normal(0, 1/math.sqrt(M), (M, C))
 62    order = np.arange(len(X))
 63    for _ in range(epochs):
 64        if not ordered:
 65            rng.shuffle(order)
 66        for p in range(len(order)-1):
 67            i, j = order[p], order[p+1]
 68            x, xn = X[i], X[j]
 69            h = sigmoid(W @ x + b)
 70            hn = sigmoid(W @ xn + b)
 71            u = A.T @ (A @ (hn-h)) - lam*(h-rho)
 72            gain = h*(1-h)
 73            d = gain*u
 74            W += eta*np.outer(d, x)
 75            b += etab*d
 76    return W, b
 77
 78
 79def bp_train(X, C=24, epochs=3, lr=.08, seed=0):
 80    # Backprop baseline: minimize consecutive representation mismatch plus
 81    # weak variance/homeostatic regularization, using the same ordered stream.
 82    rng = np.random.RandomState(seed)
 83    n = X.shape[1]
 84    W = rng.normal(0, .18, (C, n)); b = np.zeros(C)
 85    for _ in range(epochs):
 86        for i in range(len(X)-1):
 87            x, xn = X[i], X[i+1]
 88            z=W@x+b; zn=W@xn+b
 89            h=sigmoid(z); hn=sigmoid(zn)
 90            diff=hn-h
 91            # loss .5||hn-h||^2 + .08 mean((h-.5)^2), update both pair terms
 92            gh = -diff + .16*(h-.5)/C
 93            ghn = diff + .16*(hn-.5)/C
 94            dz=gh*h*(1-h); dzn=ghn*hn*(1-hn)
 95            W -= lr*(np.outer(dz,x)+np.outer(dzn,xn))
 96            b -= lr*(dz+dzn)
 97    return W,b
 98
 99
100def features(X, W, b):
101    return sigmoid(X @ W.T + b)
102
103
104def probe_accuracy(Ftr, ytr, Fte, yte):
105    clf = LogisticRegression(max_iter=300, random_state=SEED)
106    clf.fit(Ftr, ytr)
107    return float(accuracy_score(yte, clf.predict(Fte)))
108
109
110def experiment():
111    d = load_digits()
112    X = d.data.astype(np.float64)/16.0
113    y = d.target
114    # Temporal stream: sort by label, with small within-class noise permutation,
115    # making adjacency meaningful but not identical. Hold out a fixed suffix.
116    rng=np.random.RandomState(SEED)
117    tr, te = train_test_split(np.arange(len(X)), test_size=.3, stratify=y,
118                              random_state=SEED)
119    tr=np.array(tr); te=np.array(te)
120    stream=np.concatenate([tr[np.argsort(y[tr])], te])
121    Xtr, ytr = X[stream[:len(tr)]], y[stream[:len(tr)]]
122    Xte, yte = X[te], y[te]
123    results={}
124    configs=[('local_ordered',True),('local_shuffled',False)]
125    for name, ordered in configs:
126        W,b=local_train(Xtr, ordered=ordered, seed=SEED)
127        Ft,Fe=features(Xtr,W,b),features(Xte,W,b)
128        results[name]={'accuracy':probe_accuracy(Ft,ytr,Fe,yte),
129          'mean_activity':float(Fe.mean()),
130          'saturated_fraction':float(np.mean((Fe<.05)|(Fe>.95))),
131          'temporal_mse':float(np.mean((Ft[1:]-Ft[:-1])**2))}
132    W,b=bp_train(Xtr,seed=SEED)
133    Ft,Fe=features(Xtr,W,b),features(Xte,W,b)
134    results['backprop_ordered']={'accuracy':probe_accuracy(Ft,ytr,Fe,yte),
135      'mean_activity':float(Fe.mean()),
136      'saturated_fraction':float(np.mean((Fe<.05)|(Fe>.95))),
137      'temporal_mse':float(np.mean((Ft[1:]-Ft[:-1])**2))}
138    return results
139
140
141if __name__ == '__main__':
142    result={'math_checks':math_checks(), 'experiment':experiment()}
143    with open('results.json','w') as f: json.dump(result,f,indent=2)
144    print(json.dumps(result,indent=2))