リザバーコンピューティングの種類を図解

import numpy as np

# 状態 X（各行が時刻ごとの状態ベクトル）と教師 d から読み出し重みを1回で求める
def train_readout(X, d, beta=1e-6):
    """X:(T,N)  d:(T,)  beta:正則化係数"""
    N = X.shape[1]
    W_out = np.linalg.solve(X.T @ X + beta * np.eye(N), X.T @ d)   # リッジ回帰の閉形式
    return W_out                                                  # 予測は X @ W_out

import numpy as np

# ---------- 1. データ：Mackey-Glass 時系列 ----------
def mackey_glass(n, tau=17):
    x = np.full(n + tau, 1.2)
    for t in range(tau, n + tau - 1):
        x[t + 1] = x[t] + 0.2 * x[t - tau] / (1 + x[t - tau] ** 10) - 0.1 * x[t]
    return x[tau:]

data = mackey_glass(4000)
n_train = 3000
mu, sd = data[:n_train].mean(), data[:n_train].std()
data = (data - mu) / sd                       # 学習区間の統計だけで標準化
u, d = data[:-1], data[1:]                    # 入力＝現在値、教師＝1ステップ先

# ---------- 2. リング（環）を作る：これだけ ----------
N       = 200        # ノード数
r       = 0.9        # 環に沿った1つの結合強度（＝スペクトル半径そのもの）
v       = 0.1        # 入力の大きさ（全ノード共通）
beta    = 1e-8       # リッジ回帰の正則化係数
washout = 200

pi_digits = (
    "14159265358979323846264338327950288419716939937510"
    "58209749445923078164062862089986280348253421170679"
    "82148086513282306647093844609550582231725359408128"
    "48111745028410270193852110555964462294895493038196"
    "44288109756659334461284756482337867831652712019091"
    "45648566923460348610454326648213393607260249141273"
    "72458700660631558817488152092096282925409171536436"
    "78925903600113305305488204665213841469519415116094"
)
signs = np.array([1.0 if c >= "5" else -1.0 for c in pi_digits[:N]])  # πの小数桁で決定
W_in  = v * signs    # 入力重み：大きさは共通、符号パターンのみ

# ---------- 3. 状態を集める（環に沿って前のノードの値を1つ受け取る） ----------
T = len(u)
X = np.zeros((T, N + 1))                       # 先頭列はバイアス
x = np.zeros(N)
for t in range(T):
    x = np.tanh(W_in * u[t] + r * np.roll(x, 1))  # np.roll(x,1) が「環の1つ隣」
    X[t, 0] = 1.0
    X[t, 1:] = x

# ---------- 4. 読み出し層をリッジ回帰で学習（他方式と全く同じ） ----------
Xtr, dtr = X[washout:n_train], d[washout:n_train]
W_out = np.linalg.solve(Xtr.T @ Xtr + beta * np.eye(N + 1), Xtr.T @ dtr)

# ---------- 5. テスト区間で評価 ----------
y = X[n_train:] @ W_out
nrmse = np.sqrt(np.mean((y - d[n_train:]) ** 2)) / d[n_train:].std()
print(f"リングリザバー テストNRMSE: {nrmse:.4f}")   # 0.0003 前後

import numpy as np
rng = np.random.default_rng(0)

# ---------- 設定 ----------
N        = 200          # ニューロン数
n_class  = 3            # クラス数（3種類の周波数）
n_trials = 600          # 試行数
L        = 80           # 1試行の長さ（ステップ）
noise    = 0.25         # 入力ノイズ
freqs    = np.linspace(1/28, 1/10, n_class)   # 各クラス＝ある周波数

# ---------- リザバー（作ったら固定）：多様な時定数がカギ ----------
tau_m = rng.uniform(5, 16, N)                 # 膜時定数（ニューロンごとに違える）
a_s   = np.exp(-1.0 / rng.uniform(3, 25, N))  # シナプス痕跡の減衰（記憶の長さを決める）
th    = rng.uniform(0.9, 1.1, N)              # 発火閾値
bias  = rng.uniform(0.75, 1.0, N) * 1.15      # 閾値近くに置くバイアス電流
W_in  = rng.uniform(-1, 1, N) * 2.2
W = rng.normal(0, 1, (N, N)); W[rng.random((N, N)) > 0.12] = 0.0
inh = rng.random(N) < 0.2                      # 2割を抑制性に
W[:, inh] = -np.abs(W[:, inh]); W[:, ~inh] = np.abs(W[:, ~inh])
W *= 0.9 / max(abs(np.linalg.eigvals(W)))      # 再帰結合のスケール調整

# ---------- 1試行を流して「液体状態」を取り出す ----------
def run_trial(sig, tref=2):
    v = np.zeros(N); refr = np.zeros(N, int); trace = np.zeros(N); spk = np.zeros(N)
    for t in range(len(sig)):
        I = bias + W_in * sig[t] + W @ spk     # 入力電流
        active = refr <= 0
        v = v + active * (1.0 / tau_m) * (-v + I)   # 漏れ積分（LIF）
        spk = np.zeros(N); fire = active & (v >= th); spk[fire] = 1.0
        v[fire] = 0.0; refr[fire] = tref; refr = np.maximum(refr - 1, 0)  # 発火→リセット→不応期
        trace = trace * a_s + spk              # スパイクを指数フィルタ＝液体状態
    return trace                               # 試行末の状態ベクトル

# ---------- データ生成＋状態収集 ----------
feats, labels = [], []
for _ in range(n_trials):
    cls = rng.integers(n_class); ph = rng.uniform(0, 2*np.pi)
    sig = np.sin(2*np.pi*freqs[cls]*np.arange(L) + ph) + noise*rng.standard_normal(L)
    feats.append(run_trial(sig)); labels.append(cls)
X = np.array(feats); y = np.array(labels)

# ---------- 線形読み出し（リッジ回帰）で分類：他方式と同じ型 ----------
ntr = int(0.7 * n_trials)
F  = np.concatenate([np.ones((ntr, 1)), X[:ntr]], axis=1)
Ft = np.concatenate([np.ones((n_trials - ntr, 1)), X[ntr:]], axis=1)
Y  = np.eye(n_class)[y[:ntr]]                  # 教師をワンホットに
W_out = np.linalg.solve(F.T @ F + 1e-1*np.eye(F.shape[1]), F.T @ Y)
pred = np.argmax(Ft @ W_out, axis=1)
acc = np.mean(pred == y[ntr:])
print(f"スパイキングリザバー テスト精度: {acc*100:.1f}%  (3クラス・チャンス33%)")  # 96%前後

import numpy as np

def mackey_glass(n, tau=17):
    x = np.full(n + tau, 1.2)
    for t in range(tau, n + tau - 1):
        x[t + 1] = x[t] + 0.2 * x[t - tau] / (1 + x[t - tau] ** 10) - 0.1 * x[t]
    return x[tau:]

data = mackey_glass(4000); n_train = 3000
mu, sd = data[:n_train].mean(), data[:n_train].std()
data = (data - mu) / sd
u, d = data[:-1], data[1:]

# ---------- 振動子リザバー（作ったら固定） ----------
N       = 100      # 振動子の数
K       = 1.0      # 結合の強さ（同期のしやすさ）
eps     = 0.6      # 入力の強さ
a       = 1.0      # 各振動子の利得（振幅を生む）
dt      = 0.1      # 積分刻み
nsub    = 5        # 1入力あたりの内部積分ステップ
beta    = 1e-6
washout = 200

rng   = np.random.default_rng(0)
omega = rng.normal(0, 0.5, N)          # 自然周波数（振動子ごとにばらつく）
c_in  = rng.uniform(-1, 1, N)          # 入力の配り方
z     = 0.1 * (rng.standard_normal(N) + 1j * rng.standard_normal(N))  # 複素振幅

# ---------- 状態を集める ----------
T = len(u); X = np.zeros((T, 2 * N + 1))
for t in range(T):
    drive = eps * c_in * u[t]
    for _ in range(nsub):
        zbar = z.mean()                                  # 平均場（結合の中身）
        dz = (a + 1j * omega - np.abs(z) ** 2) * z + K * (zbar - z) + drive
        z = z + dt * dz                                  # −|z|²z が振幅を安定化
    X[t, 0] = 1.0
    X[t, 1:N+1]  = z.real                                 # 状態＝各振動子の実部・虚部
    X[t, N+1:]   = z.imag

# ---------- 読み出し層をリッジ回帰で学習：他方式と同じ型 ----------
Xtr, dtr = X[washout:n_train], d[washout:n_train]
M = X.shape[1]
W_out = np.linalg.solve(Xtr.T @ Xtr + beta * np.eye(M), Xtr.T @ dtr)
y = X[n_train:] @ W_out
nrmse = np.sqrt(np.mean((y - d[n_train:]) ** 2)) / d[n_train:].std()
print(f"結合振動子リザバー テストNRMSE: {nrmse:.4f}")   # 0.012 前後

import numpy as np
from itertools import combinations_with_replacement

# ---------- 1. データ：ローレンツ63（RK4で数値積分） ----------
def lorenz(n, dt=0.025, s=10., r=28., b=8/3., warm=1000):
    X = np.zeros((n + warm, 3)); X[0] = [1., 1., 1.]
    f = lambda v: np.array([s*(v[1]-v[0]), v[0]*(r-v[2])-v[1], v[0]*v[1]-b*v[2]])
    for t in range(n + warm - 1):
        v = X[t]; k1=f(v); k2=f(v+dt/2*k1); k3=f(v+dt/2*k2); k4=f(v+dt*k3)
        X[t+1] = v + dt/6*(k1 + 2*k2 + 2*k3 + k4)
    return X[warm:]

dt = 0.025; n = 6000; n_train = 4000; warm = 100; beta = 1e-3
k, s = 2, 1                                    # 時間遅れ：2タップ、間隔1
data = lorenz(n, dt=dt)
mu, sd = data[:n_train].mean(0), data[:n_train].std(0)
D = (data - mu) / sd

# ---------- 2. 特徴ベクトル：定数 + 線形（遅れ） + 2次の積 ----------
def features(taps):                             # taps: 形 (k, 3)
    lin = taps.flatten()                        # 線形部（3k 次元）
    quad = np.array([lin[i]*lin[j]              # 2次の単項式（重複ありの組合せ）
                     for i, j in combinations_with_replacement(range(len(lin)), 2)])
    return np.concatenate([[1.0], lin, quad])   # ← ここが「リザバーの代わり」

span = (k - 1) * s
F = np.array([features(D[t-span:t+1:s]) for t in range(span, n-1)])
Y = np.array([D[t+1] - D[t] for t in range(span, n-1)])   # 教師＝1ステップ差分

# ---------- 3. リッジ回帰1回で学習：他方式と同じ型 ----------
Ftr, Ytr = F[warm:n_train], Y[warm:n_train]
M = Ftr.shape[1]                                # 特徴次元＝28（=1+6+21）
W_out = np.linalg.solve(Ftr.T @ Ftr + beta*np.eye(M), Ftr.T @ Ytr)
pred = F[n_train:] @ W_out
nrmse = np.sqrt(np.mean((pred - Y[n_train:])**2)) / Y[n_train:].std()
print(f"1ステップ NRMSE: {nrmse:.5f}   特徴次元: {M}")     # 0.0017 前後

# ---------- 4. 自律走行：予測だけでアトラクタを再現 ----------
hist = list(D[n_train-span:n_train+1:s]); cur = D[n_train].copy(); gen = [cur.copy()]
for _ in range(1999):
    cur = cur + features(np.array(hist[-k:])) @ W_out
    gen.append(cur.copy()); hist.append(cur.copy())
gen = np.array(gen); truth = D[n_train:n_train+len(gen)]
m = min(len(gen), len(truth))
err = np.linalg.norm(gen[:m]-truth[:m], axis=1) / np.sqrt(3)
vt = int(np.argmax(err > 0.3)) if np.any(err > 0.3) else m
print(f"自律走行の有効予測: 約 {vt*dt:.1f} 時間単位（ローレンツのリアプノフ時間の数倍）")