Skip to content

Instantly share code, notes, and snippets.

@giuseppebonaccorso
Last active October 22, 2017 16:16
Show Gist options
  • Save giuseppebonaccorso/ee57c79334164dcfd6d9e5c422229e19 to your computer and use it in GitHub Desktop.
Save giuseppebonaccorso/ee57c79334164dcfd6d9e5c422229e19 to your computer and use it in GitHub Desktop.
Example of Self-Organizing Map (Kohonen Network) based on the Olivetti faces dataset (Cupy-based)
import cupy as cp
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import fetch_olivetti_faces
# Set random seed for reproducibility
np.random.seed(1000)
cp.random.seed(1000)
# Retrieve the dataset
faces = fetch_olivetti_faces()
nb_iterations = 5000
nb_startup_iterations = 100
nb_patterns = 100
pattern_length = faces['data'].shape[1]
pattern_width = pattern_height = int(np.sqrt(pattern_length))
eta0 = 1.0
sigma0 = 10
tau = 400.0
matrix_side = 20
precomputed_distances = np.zeros((matrix_side, matrix_side, matrix_side, matrix_side))
# Create the training set
X = faces['data']
X_train = cp.array(X[0:nb_patterns])
# Initialize the weights
W = cp.random.uniform(cp.min(X_train), cp.max(X_train), size=(matrix_side, matrix_side, pattern_length))
def winning_unit(xt):
distances = cp.linalg.norm(W - xt, ord=2, axis=2)
max_activation_unit = cp.argmax(distances)
return int(cp.floor(max_activation_unit / matrix_side)), max_activation_unit % matrix_side
def eta(t):
return eta0 * cp.exp(-float(t) / tau)
def sigma(t):
return float(sigma0) * cp.exp(-float(t) / tau)
# Precompute the distances
for i in range(matrix_side):
for j in range(matrix_side):
for k in range(matrix_side):
for t in range(matrix_side):
precomputed_distances[i, j, k, t] = np.power(float(i) - float(k), 2) + np.power(float(j) - float(t), 2)
precomputed_distances = cp.array(precomputed_distances)
def distance_matrix(xt, yt, sigmat):
dm = precomputed_distances[xt, yt, :, :]
de = 2.0 * cp.power(sigmat, 2)
return cp.exp(-dm / de)
if __name__ == '__main__':
# Main cycle
sequence = np.arange(0, nb_patterns)
for t in range(nb_iterations):
np.random.shuffle(sequence)
if t < nb_startup_iterations:
etat = eta(t)
sigmat = sigma(t)
else:
etat = 0.1
sigmat = 0.01
for n in sequence:
x_sample = X_train[n]
xw, yw = winning_unit(x_sample)
dm = distance_matrix(xw, yw, sigmat)
dW = etat * cp.expand_dims(dm, axis=2) * (x_sample - W)
W += dW
# Show the weight matrix
matrix_w = cp.zeros((matrix_side * pattern_height, matrix_side * pattern_width))
for i in range(matrix_side):
for j in range(matrix_side):
matrix_w[i * pattern_height:i * pattern_height + pattern_height,
j * pattern_height:j * pattern_height + pattern_width] = W[i, j].reshape((pattern_height, pattern_width)) * 255.0
fig, ax = plt.subplots(figsize=(12, 12))
ax.matshow(matrix_w.tolist(), cmap='gray')
ax.set_xticks([])
ax.set_yticks([])
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment