import numpy as np from scipy.io.wavfile import read from ran_array import ( RANArray, ran_add, ran_sub, ran_mul, const_ran, real_ran, zero_ran, ran_matmul, ran_sum0, ) # ------------------------------------------------------------ # Why this version is fast but still genuinely uses RAN: # # The original implementation wrapped every single weight, activation, # and gradient in its own scalar RAN Python object, inside triple-nested # Python loops (784 x 100 x 10 elements, per sample, per batch, per # epoch). That's billions of individual object allocations - the actual # cause of the slowness, not a math bug. # # RANArray (ran_array.py) implements the exact same a/b + c algebra - # addition, subtraction, multiplication, collapse - but stores a, b, c # as numpy arrays and applies the formulas elementwise/batched via # numpy, so an entire weight matrix updates in one vectorized call # instead of tens of thousands of Python-level RAN objects. Collapse # still only happens where the original collapsed (relu's sign check, # softmax, and the final loss/accuracy read-out) - the weights # themselves stay in raw (a, b, c) form across the whole run, so this # still gets you delayed collapse's actual benefit: the tiny per-step # SGD updates (lr * grad) accumulate into the weights without being # rounded together with the weight's current magnitude at every single # step, the way plain float += would. # ------------------------------------------------------------ # ------------------------------------------------------------ # Load data (unchanged from original) # ------------------------------------------------------------ X_train = read('../X_train.wav')[1].reshape(-1, 784).astype(np.float64) y_train = (read('../y_train.wav')[1] * 9).astype(int) X_test = read('../X_test.wav')[1].reshape(-1, 784).astype(np.float64) y_test = (read('../y_test.wav')[1] * 9).astype(int) X20 = X_test[:1000] yt20 = y_test[:1000] # ------------------------------------------------------------ # MLP Classifier - all parameters are RAN values (a/b + c), stored as # vectorized RANArrays instead of grids of scalar RAN objects. # ------------------------------------------------------------ class MLPClassifierRAN: def __init__(self, input_size, hidden_size, output_size, learning_rate=0.01, seed=None): rng = np.random.default_rng(seed) # weights as RAN(value, 1, 0) -- same convention as RAN(random.gauss(...), 1, 0) self.W1 = const_ran(rng.normal(0, 0.01, size=(input_size, hidden_size))) self.b1 = zero_ran((hidden_size,)) self.W2 = const_ran(rng.normal(0, 0.01, size=(hidden_size, output_size))) self.b2 = zero_ran((output_size,)) learning_rate = np.random.rand(4) self.lr0 = const_ran(learning_rate[0]) # RAN(learning_rate, 1, 0) self.lr1 = const_ran(learning_rate[1]) # RAN(learning_rate, 1, 0) self.lr2 = const_ran(learning_rate[2]) # RAN(learning_rate, 1, 0) self.lr3 = const_ran(learning_rate[3]) # RAN(learning_rate, 1, 0) def relu(self, z): # z: RANArray. Keep RAN structure where z.collapse() > 0, # else RAN(0,1,0) -- same rule as the original relu(). vals = z.collapse() mask = vals > 0 a = np.where(mask, z.a, 0.0) b = np.where(mask, z.b, 1.0) c = np.where(mask, z.c, 0.0) return RANArray(a, b, c), mask def softmax(self, z_ran): # softmax needs real floats for exp/sum, so we collapse here -- # exactly the same boundary the original collapsed at. vals = z_ran.collapse() shifted = vals - np.max(vals, axis=1, keepdims=True) exp_vals = np.exp(shifted) return exp_vals / np.sum(exp_vals, axis=1, keepdims=True) def forward(self, X): Xr = real_ran(X) # matches float*RAN auto-conversion RAN(0,1,x) z1 = ran_add(ran_matmul(Xr, self.W1), self.b1) a1, relu_mask = self.relu(z1) z2 = ran_add(ran_matmul(a1, self.W2), self.b2) probs = self.softmax(z2) return probs, a1, relu_mask def compute_loss(self, y_true, y_pred): return -np.mean(np.sum(y_true * np.log(y_pred + 1e-9), axis=1)) def backward(self, X, y_true, probs, a1, relu_mask): N = X.shape[0] # dz2 = RAN(probs - y_true, 1, 0), averaged over the batch # (real mini-batch gradient, vs. the original's per-sample loop) dz2 = const_ran((probs - y_true) / N) dW2 = ran_matmul(a1.T, dz2) db2 = ran_sum0(dz2) da1 = ran_matmul(dz2, self.W2.T) deriv = const_ran(relu_mask.astype(np.float64)) # RAN(relu_derivative, 1, 0) dz1 = ran_mul(da1, deriv) Xr = real_ran(X) dW1 = ran_matmul(Xr.T, dz1) db1 = ran_sum0(dz1) return dW1, db1, dW2, db2 def update(self, X, y_true): probs, a1, relu_mask = self.forward(X) dW1, db1, dW2, db2 = self.backward(X, y_true, probs, a1, relu_mask) # W -= lr * dW, done entirely in RAN space (no collapse) self.W1 = ran_sub(self.W1, ran_mul(self.lr0, dW1)) self.b1 = ran_sub(self.b1, ran_mul(self.lr1, db1)) self.W2 = ran_sub(self.W2, ran_mul(self.lr2, dW2)) self.b2 = ran_sub(self.b2, ran_mul(self.lr3, db2)) return self.compute_loss(y_true, probs) def predict(self, X): probs, _, _ = self.forward(X) return np.argmax(probs, axis=1) def score(self, X, y_true): preds = self.predict(X) return np.mean(preds == y_true) # ------------------------------------------------------------ # Training loop (real mini-batch gradient descent, vectorized) # ------------------------------------------------------------ if __name__ == "__main__": mlp = MLPClassifierRAN(input_size=784, hidden_size=256, output_size=10, learning_rate=0.01) epochs = 50 batch_size = 100 n_samples = X_train.shape[0] for epoch in range(epochs): idx = np.random.permutation(n_samples) X_shuffled = X_train[idx] y_shuffled = y_train[idx] total_loss = 0.0 n_batches = 0 for i in range(0, n_samples, batch_size//10): X_batch = X_shuffled[i:i + batch_size] y_batch = y_shuffled[i:i + batch_size] y_onehot = np.eye(10)[y_batch] loss = mlp.update(X_batch, y_onehot) total_loss += loss n_batches += 1 if i%100==0: acc = mlp.score(X20, yt20) print(f"Epoch {epoch+1}: average loss = {total_loss/n_batches:.6f}, test accuracy = {acc:.4f}")