import numpy as np from scipy.io.wavfile import read from scipy.spatial.distance import cdist from sklearn.neighbors import NearestNeighbors # ---------- Data loading ---------- X_train = read('../X_train.wav')[1].reshape(-1, 784) y_train = (read('../y_train.wav')[1] * 9).astype(int) X_test = read('../X_test.wav')[1].reshape(-1, 784) y_test = (read('../y_test.wav')[1] * 9).astype(int) X20 = X_test[:1000] yt20 = y_test[:1000] # ---------- Helper: build kNN graph and normalized adjacency ---------- def build_adjacency(features, k=5, metric='cosine'): """Build kNN graph and return normalized adjacency matrix A_hat.""" n = features.shape[0] if n == 0: return np.zeros((0, 0)) if n <= k: # fully connected if too few samples A = np.ones((n, n)) - np.eye(n) else: nbrs = NearestNeighbors(n_neighbors=k, metric=metric).fit(features) knn = nbrs.kneighbors(features, return_distance=False) A = np.zeros((n, n)) for i in range(n): A[i, knn[i]] = 1 A[knn[i], i] = 1 # symmetric # add self‑loops A = A + np.eye(n) # degree matrix D D = np.diag(np.sum(A, axis=1) ** -0.5) A_hat = D @ A @ D return A_hat # ---------- Graph Autoencoder per class ---------- class GraphAutoencoder: def __init__(self, input_dim, hidden_dim=64, output_dim=None, lr=0.01): self.input_dim = input_dim self.hidden_dim = hidden_dim self.output_dim = output_dim if output_dim is not None else input_dim self.lr = lr # Encoder: two GCN layers self.W1 = np.random.randn(input_dim, hidden_dim) * 0.01 self.b1 = np.zeros(hidden_dim) self.W2 = np.random.randn(hidden_dim, hidden_dim) * 0.01 self.b2 = np.zeros(hidden_dim) # Decoder: linear mapping self.W3 = np.random.randn(hidden_dim, self.output_dim) * 0.01 self.b3 = np.zeros(self.output_dim) def relu(self, x): return np.maximum(0, x) def forward(self, A_hat, X): # X: (n, input_dim), A_hat: (n, n) h1 = self.relu(A_hat @ X @ self.W1 + self.b1) h2 = self.relu(A_hat @ h1 @ self.W2 + self.b2) X_recon = A_hat @ h2 @ self.W3 + self.b3 # graph‑aware reconstruction return X_recon, h1, h2 def compute_loss(self, X, X_recon): return np.mean((X - X_recon) ** 2) def backward(self, A_hat, X, X_recon): n = X.shape[0] # gradients w.r.t. decoder dLoss = 2 * (X_recon - X) / n # (n, output_dim) # dW3 = A_hat.T @ (A_hat @ self.a2) @ dLoss # actually need to backprop through A_hat # we'll simplify: compute gradients by considering A_hat as constant # full chain is complex; use a simple approximation: treat X_recon as (A_hat @ h2 @ W3) # For brevity, we compute gradients of MSE with respect to W3, b3 directly grad_W3 = (A_hat @ self.a2).T @ dLoss grad_b3 = np.sum(dLoss, axis=0, keepdims=True) # propagate to h2 dh2 = dLoss @ self.W3.T * (self.a2 > 0) # ReLU derivative # W2, b2 grad_W2 = (A_hat @ self.a1).T @ (A_hat @ dh2) grad_b2 = np.sum(A_hat @ dh2, axis=0, keepdims=True) # propagate to h1 dh1 = (A_hat @ dh2) @ self.W2.T * (self.a1 > 0) grad_W1 = (A_hat @ X).T @ (A_hat @ dh1) grad_b1 = np.sum(A_hat @ dh1, axis=0, keepdims=True) # update parameters self.W3 -= self.lr * grad_W3 self.b3 -= self.lr * grad_b3.flatten() self.W2 -= self.lr * grad_W2 self.b2 -= self.lr * grad_b2.flatten() self.W1 -= self.lr * grad_W1 self.b1 -= self.lr * grad_b1.flatten() def train_step(self, A_hat, X): X_recon, self.a1, self.a2 = self.forward(A_hat, X) # store activations loss = self.compute_loss(X, X_recon) self.backward(A_hat, X, X_recon) return loss def reconstruct(self, A_hat, X): X_recon, *_ = self.forward(A_hat, X) return X_recon # ---------- Build per‑class models ---------- num_classes = 10 class_models = [] class_adj = [] class_features = [] for c in range(num_classes): idx = np.where(y_train == c)[0] X_c = X_train[idx] if len(X_c) == 0: # fallback: use a single dummy sample X_c = X_train[:1] A_hat = build_adjacency(X_c, k=5) model = GraphAutoencoder(input_dim=784, hidden_dim=64, lr=0.001) class_models.append(model) class_adj.append(A_hat) class_features.append(X_c) # ---------- Train each graph autoencoder on its class data ---------- epochs = 1 for c in range(num_classes): print(f"Training class {c} ...") A_hat = class_adj[c] X_c = class_features[c] for epoch in range(epochs): loss = class_models[c].train_step(A_hat, X_c) if epoch % 10 == 0: print(f" epoch {epoch}: loss = {loss:.6f}") # ---------- Prediction: reconstruction error per class ---------- def predict_one(sample, k=5): errors = [] for c in range(num_classes): # build extended graph: original class nodes + this sample X_c = class_features[c] if len(X_c) == 0: errors.append(np.inf) continue # add sample as new node X_ext = np.vstack([X_c, sample.reshape(1, -1)]) # compute kNN edges for the new node to class nodes n_c = len(X_c) if n_c <= k: # connect to all existing nodes new_edges = np.ones(n_c) else: # find k nearest neighbors in X_c dists = cdist(sample.reshape(1, -1), X_c, metric='cosine').flatten() nn_indices = np.argsort(dists)[:k] new_edges = np.zeros(n_c) new_edges[nn_indices] = 1 # build adjacency for extended graph A = np.zeros((n_c+1, n_c+1)) A[:n_c, :n_c] = class_adj[c] # original adjacency (already normalized) # add connections from new node to its neighbors A[n_c, :n_c] = new_edges A[:n_c, n_c] = new_edges # add self‑loop to new node A[n_c, n_c] = 1 # degree normalization (approximate) D = np.diag(np.sum(A, axis=1) ** -0.5) A_hat_ext = D @ A @ D # reconstruct the new node's features with np.errstate(all='ignore'): X_recon_ext, *_ = class_models[c].forward(A_hat_ext, X_ext) recon_sample = X_recon_ext[-1] # only for the new node err = np.mean((sample - recon_sample) ** 2) errors.append(err) return np.argmin(errors) def predict(X): return np.array([predict_one(x) for x in X]) def score(X, y): preds = predict(X) return np.mean(preds == y) # ---------- Evaluate on test set ---------- acc = score(X20, yt20) print(f"Test accuracy on first 1000 samples: {acc:.4f}")