import numpy as np
import matplotlib.pyplot as plt

# Return the (g,h) index of the BMU in the grid
def find_BMU(SOM,x):
    distSq = (np.square(SOM - x)).sum(axis=2)
    return np.unravel_index(np.argmin(distSq, axis=None), distSq.shape)

# Update the weights of the SOM cells when given a single training example
# and the model parameters along with BMU coordinates as a tuple
def update_weights(SOM, train_ex, learn_rate, radius_sq,
                   BMU_coord, step=3):
    g, h = BMU_coord
    #if radius is close to zero then only BMU is changed
    if radius_sq < 1e-3:
        SOM[g,h,:] += learn_rate * (train_ex - SOM[g,h,:])
        return SOM
    # Change all cells in a small neighborhood of BMU
    for i in range(max(0, g-step), min(SOM.shape[0], g+step)):
        for j in range(max(0, h-step), min(SOM.shape[1], h+step)):
            dist_sq = np.square(i - g) + np.square(j - h)
            dist_func = np.exp(-dist_sq / 2 / radius_sq)
            SOM[i,j,:] += learn_rate * dist_func * (train_ex - SOM[i,j,:])
    return SOM

# Main routine for training an SOM. It requires an initialized SOM grid
# or a partially trained grid as parameter
def train_SOM(SOM, train_data, learn_rate = .01, radius_sq = 1,
             lr_decay = .1, radius_decay = .1, epochs = 10):
    learn_rate_0 = learn_rate
    radius_0 = radius_sq
    for epoch in np.arange(0, epochs):
        rand.shuffle(train_data)
        for train_ex in train_data:
            g, h = find_BMU(SOM, train_ex)
            SOM = update_weights(SOM, train_ex,
                                 learn_rate, radius_sq, (g,h))
        # Update learning rate and radius
        learn_rate = learn_rate_0 * np.exp(-epoch * lr_decay)
        radius_sq = radius_0 * np.exp(-epoch * radius_decay)
    return SOM

# Dimensions of the SOM grid
m = 10
n = 10
# Number of training examples
n_x = 400
rand = np.random.RandomState(42)
# Initialize the training data
train_data = rand.randint(0, 255, (n_x, 3))
# Initialize the SOM randomly
SOM = rand.randint(0, 255, (m, n, 3)).astype(float)
# Display both the training matrix and the SOM grid
fig, ax = plt.subplots(
    nrows=1, ncols=2, figsize=(12, 3.5),
    subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow(train_data.reshape(20, 20, 3))
ax[0].title.set_text('Training Data')
somMap(ax[1],SOM/255,m,n,0)
#ax[1].imshow(SOM.astype(int))
#ax[1].title.set_text('Randomly Initialized SOM Grid')

head = train_data[0:5,:]
print(head)
plt.imshow(head.reshape(5, 1, 3))

[[102 179  92]
 [ 14 106  71]
 [188  20 102]
 [121 210 214]
 [ 74 202  87]]

<matplotlib.image.AxesImage at 0x7df878854130>

fig, ax = plt.subplots(
    nrows=1, ncols=5, figsize=(18, 3.5),
    subplot_kw=dict(xticks=[], yticks=[]))
total_epochs = 0
for epochs, i in zip([1, 9, 10, 30, 50], range(0,5)):
    total_epochs += epochs
    SOM = train_SOM(SOM, train_data, epochs=epochs)
    ax[i].imshow(SOM.astype(int))
    ax[i].title.set_text('Epochs = ' + str(total_epochs))

print(SOM[1,1,:])

[ 97.81524813 150.47228211 216.97180575]

fig, ax = plt.subplots(
    nrows=1, ncols=5, figsize=(18, 3.5),
    subplot_kw=dict(xticks=[], yticks=[]))
total_epochs = 0
for epochs, i in zip([1, 9, 10, 30, 50], range(0,5)):
    total_epochs += epochs
    SOM = train_SOM(SOM, train_data, epochs=epochs)
    somMap(ax[i],SOM/255,m,n,total_epochs)
    #ax[i].imshow(SOM.astype(int))
    #ax[i].title.set_text('Epochs = ' + str(total_epochs))
plt.show()

def somMap(ax,SOM,m,n,total_epochs):
    for i in range(m):
      for j in range(n):
        ax.add_artist(plt.Circle((i, j), radius=0.5, facecolor=SOM[i, j]))
    ax.axis('image')
    # Set the limits of the axes
    ax.set_xlim(-0.5, m-0.5)
    ax.set_ylim(-0.5, n-0.5)
    ticks = range(0,n)
    #print(ticks)
    ax.set_xticks(ticks)
    ax.set_xticklabels(ticks)
    ax.set_yticks(ticks)
    ax.set_yticklabels(ticks)
    # Set the title and labels
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.title.set_text('Epochs = ' + str(total_epochs))

    # Show the plot
    #plt.show()

    return

import matplotlib.pyplot as plt
import numpy.random as rnd

m = 10
n = 10
# Generate 5 by 5 random colors
colors = rnd.rand(m,n, 3)

plt.tight_layout()
# Create a figure and axes
fig, ax = plt.subplots()

# Plot the circles
for i in range(m):
    for j in range(n):
        ax.add_artist(plt.Circle((i, j), radius=0.5, facecolor=colors[i, j]))
ax.axis('image')
# Set the limits of the axes
ax.set_xlim(-0.5, m-0.5)
ax.set_ylim(-0.5, n-0.5)
ticks = range(0,n)
#print(ticks)
ax.set_xticks(ticks)
ax.set_xticklabels(ticks)
ax.set_yticks(ticks)
ax.set_yticklabels(ticks)
# Set the title and labels
ax.set_xlabel("X")
ax.set_ylabel("Y")


# Show the plot
plt.show()

<Figure size 640x480 with 0 Axes>

epochs = np.arange(0, 50)
lr_decay = [0.01, 0.1, 0.5, 0.99]
fig,ax = plt.subplots(nrows=1, ncols=4, figsize=(15,4))
plt_ind = np.arange(4) + 141
for decay, ind in zip(lr_decay, plt_ind):
    plt.subplot(ind)
    learn_rate = np.exp(-epochs * decay)
    plt.plot(epochs, learn_rate, c='black')
    plt.title('decay rate: ' + str(decay))
    plt.xlabel('epochs $t$')
    plt.ylabel('$\eta^(t)$')
fig.subplots_adjust(hspace=0.5, wspace=0.3)
plt.show()

distance = np.arange(0, 30)
sigma_sq = [100,10,1,0.1]
fig,ax = plt.subplots(nrows=1, ncols=4, figsize=(15,4))
plt_ind = np.arange(4) + 141
for s, ind in zip(sigma_sq, plt_ind):
    plt.subplot(ind)
    f = np.exp(-distance ** 2 / 2 / s)
    plt.plot(distance, f, c='black')
    plt.title('$\sigma^2$ = ' + str(s))
    plt.xlabel('Distance')
    plt.ylabel('Neighborhood function $f$')
fig.subplots_adjust(hspace=0.5, wspace=0.3)
plt.show()