""" esn_baseline.py =============== Tuned Standard ESN baseline for comparison against K-R architecture. Includes: StandardESN -- Single-reservoir ESN with ridge regression DeepESN2Layer -- 2-layer DeepESN baseline (Gallicchio & Micheli, 2017) grid_search -- Held-out hyperparameter search evaluate_10seed -- 10-seed evaluation returning mean +/- std Usage: from esn_baseline import StandardESN, grid_search, evaluate_10seed best_params = grid_search(u, y, method='std') mean, std, scores = evaluate_10seed(u, y, 'std', best_params) """ import numpy as np import itertools from typing import Dict, List, Tuple, Optional # ── Reservoir construction ───────────────────────────────────── def make_reservoir(N: int, rho: float, seed: int) -> np.ndarray: np.random.seed(seed) W = np.random.randn(N, N) W *= rho / (np.max(np.abs(np.linalg.eigvals(W))) + 1e-12) return W def make_input_weights(N: int, scale: float, seed: int) -> np.ndarray: np.random.seed(seed + 1000) return np.random.randn(N) * scale # ── Standard ESN ─────────────────────────────────────────────── class StandardESN: """ Standard Echo State Network with ridge regression readout. Parameters ---------- N : int Reservoir size rho : float Spectral radius sigma : float Input scaling lam : float Ridge regularisation seed : int Random seed washout : int Washout steps """ def __init__(self, N=200, rho=0.9, sigma=0.1, lam=1e-4, seed=0, washout=200): self.N = N; self.rho = rho; self.sigma = sigma self.lam = lam; self.seed = seed; self.washout = washout self.W = make_reservoir(N, rho, seed) self.Win = make_input_weights(N, sigma, seed) self.W_out = None self.mu = None; self.sg = None def _states(self, u): x = np.zeros(self.N); Xs = [] for t in range(len(u)): x = np.tanh(self.W @ x + self.Win * u[t]) Xs.append(x.copy()) return np.array(Xs) def fit(self, u, y): X = self._states(u) self.mu = X[self.washout:].mean(0) self.sg = X[self.washout:].std(0) + 1e-8 X = (X - self.mu) / self.sg Xtr = X[self.washout:]; ytr = y[self.washout:] A = Xtr.T @ Xtr + self.lam * np.eye(self.N) self.W_out = np.linalg.solve(A, Xtr.T @ ytr) return self def evaluate(self, u, y, test_size=500): X = self._states(u) X = (X - self.mu) / self.sg Xte = X[-test_size:]; yte = y[-test_size:] yhat = Xte @ self.W_out nrmse = float(np.sqrt(np.mean((yte - yhat) ** 2)) / (np.std(yte) + 1e-12)) return nrmse, yte, yhat # ── 2-Layer DeepESN ──────────────────────────────────────────── class DeepESN2Layer: """ 2-layer DeepESN (Gallicchio et al., 2018). Layer 2 receives layer 1 states as input. Parameters ---------- N : int Neurons per layer (total = 2N features) rho1 : float Layer 1 spectral radius rho2 : float Layer 2 spectral radius lam : float Ridge regularisation seed : int Random seed """ def __init__(self, N=200, rho1=0.9, rho2=0.7, lam=1e-4, seed=0, washout=200): self.N = N; self.rho1 = rho1; self.rho2 = rho2 self.lam = lam; self.seed = seed; self.washout = washout self.W1 = make_reservoir(N, rho1, seed) self.Win1 = make_input_weights(N, 0.5, seed) self.W2 = make_reservoir(N, rho2, seed + 200) self.Win2 = make_input_weights(N, 0.5, seed + 300) self.W_out = None self.mu = None; self.sg = None def _states(self, u): x1 = np.zeros(self.N); x2 = np.zeros(self.N); Xs = [] for t in range(len(u)): x1 = np.tanh(self.W1 @ x1 + self.Win1 * u[t]) x2 = np.tanh(self.W2 @ x2 + self.Win2 @ x1) Xs.append(np.concatenate([x1, x2])) return np.array(Xs) def fit(self, u, y): X = self._states(u) self.mu = X[self.washout:].mean(0) self.sg = X[self.washout:].std(0) + 1e-8 X = (X - self.mu) / self.sg Xtr = X[self.washout:]; ytr = y[self.washout:] A = Xtr.T @ Xtr + self.lam * np.eye(X.shape[1]) self.W_out = np.linalg.solve(A, Xtr.T @ ytr) return self def evaluate(self, u, y, test_size=500): X = self._states(u) X = (X - self.mu) / self.sg Xte = X[-test_size:]; yte = y[-test_size:] yhat = Xte @ self.W_out nrmse = float(np.sqrt(np.mean((yte - yhat) ** 2)) / (np.std(yte) + 1e-12)) return nrmse, yte, yhat # ── Grid Search ──────────────────────────────────────────────── RHO_GRID = [0.5, 0.7, 0.9, 0.95, 0.99] LAM_GRID = [1e-8, 1e-6, 1e-4, 1e-2, 1.0, 10.0] SIGMA_GRID = [0.1, 0.5, 1.0, 1.5, 2.0] N_DEFAULT = 200 KF_DEFAULT = 1.5 * (200 / N_DEFAULT) ** 0.40 KS_DEFAULT = 1.0 * (200 / N_DEFAULT) ** 0.40 def grid_search(u: np.ndarray, y: np.ndarray, method: str = 'std', N: int = N_DEFAULT, M: int = 50, washout: int = 200, test_size: int = 200, val_seed: int = 99) -> Dict: """ Held-out hyperparameter grid search. Uses first half of data as validation set with val_seed. Parameters ---------- u : ndarray Full input series y : ndarray Full target series method : str 'std' | 'kr' | 'deep2' N : int Reservoir size M : int Task memory order (for K-R only) washout : int Washout steps test_size: int Validation test size val_seed : int Fixed validation seed (default 99) Returns ------- Dict with keys: rho, rho_f, rho_s, sigma, lam, method, nrmse_val """ from kr_reservoir import KRReservoir, adaptive_delays half = len(u) // 2 u_v, y_v = u[:half], y[:half] best = (np.inf, None) if method == 'std': for rho, sigma, lam in itertools.product(RHO_GRID, SIGMA_GRID, LAM_GRID): m = StandardESN(N=N, rho=rho, sigma=sigma, lam=lam, seed=val_seed, washout=washout) m.fit(u_v, y_v) n, _, _ = m.evaluate(u_v, y_v, test_size) if n < best[0]: best = (n, dict(rho=rho, sigma=sigma, lam=lam, method='std', nrmse_val=n)) elif method == 'kr': KF = 1.5 * (200 / N) ** 0.40 KS = 1.0 * (200 / N) ** 0.40 for rho_f, rho_s, lam in itertools.product(RHO_GRID, RHO_GRID, LAM_GRID): m = KRReservoir(N=N, M=M, rho_f=rho_f, rho_s=rho_s, lam=lam, seed=val_seed, washout=washout) m.fit(u_v, y_v) n, _, _ = m.evaluate(u_v, y_v, test_size) if n < best[0]: best = (n, dict(rho_f=rho_f, rho_s=rho_s, lam=lam, method='kr', nrmse_val=n)) elif method == 'deep2': for rho1, rho2, lam in itertools.product(RHO_GRID, RHO_GRID, LAM_GRID): m = DeepESN2Layer(N=N, rho1=rho1, rho2=rho2, lam=lam, seed=val_seed, washout=washout) m.fit(u_v, y_v) n, _, _ = m.evaluate(u_v, y_v, test_size) if n < best[0]: best = (n, dict(rho1=rho1, rho2=rho2, lam=lam, method='deep2', nrmse_val=n)) return best[1] def evaluate_10seed(u: np.ndarray, y: np.ndarray, method: str, params: Dict, N: int = N_DEFAULT, M: int = 50, washout: int = 200, test_size: int = 500, n_seeds: int = 10) -> Tuple[float, float, List[float]]: """ Evaluate method over n_seeds random seeds using fixed params. Returns ------- mean : float std : float scores: List[float] Per-seed NRMSE values """ from kr_reservoir import KRReservoir scores = [] for seed in range(n_seeds): if method == 'std': m = StandardESN(N=N, rho=params['rho'], sigma=params['sigma'], lam=params['lam'], seed=seed, washout=washout) elif method == 'kr': m = KRReservoir(N=N, M=M, rho_f=params['rho_f'], rho_s=params['rho_s'], lam=params['lam'], seed=seed, washout=washout) elif method == 'deep2': m = DeepESN2Layer(N=N, rho1=params['rho1'], rho2=params['rho2'], lam=params['lam'], seed=seed, washout=washout) m.fit(u, y) n, _, _ = m.evaluate(u, y, test_size) scores.append(n) return float(np.mean(scores)), float(np.std(scores)), scores