<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 K-Means Clustering
The UCI Machine Learning Repository is a "go-to" place for fascinating data sets often used in machine learning education and research.  
Let's check out the Stone Flakes data set.  Find and download the StoneFlakes.dat file by clicking here or by running wget like below.
Step1:  Let's look at the first few lines.
Step2:  Read about the column names and the meaning of the ID values at the data set's web site.  
Notice that values are separated by commas, except for the first column.  Also notice that there are question marks where data is missing.  How can we read this?  Well, the usual answer is to "google" for the answer.  Try seaching for "read data set numpy"
numpy includes functions for reading csv files.  However, the pandas package offers more powerful functions for parsing data.  Let's try its read_csv function.
Step3:  Let's just replace commas with spaces, using unix.  Read about the tr unix command at Linux TR Command Examples
Step4:  To see this data, let's try plotting each column as a separate curve on the same axes.
Step5:  Each sample has 8 attribues, so each sample is a point in 8-dimensional space.  I wonder how well the samples "clump" in those 8 dimensions.  Let's try clustering them with the k-means algorithm.
First, let's try to find two clusters, so $k=2$.  We must initialize the two means of the two clusters.  Let's just pick two samples at random.
Step6:  Now we must find all samples that are closest to the first center, and those that are closest to the second sample.
Check out the wonders of numpy broadcasting.
Step7:  But what if we want to subtract every element of a with every element of b?
Step8:  However, we can ask numpy to do this duplication for us if we reshape a to be a column vector and leave b as a row vector.
$$ \begin{pmatrix}
1\
2\
3
\end{pmatrix}
-
\begin{pmatrix}
10 & 20 & 30
\end{pmatrix}
\;\; = \;\;
\begin{pmatrix}
1 & 1 & 1\
2 & 2 & 2\
3 & 3 & 3
\end{pmatrix}
-
\begin{pmatrix}
10 & 20 & 30\
10 & 20 & 30\
10 & 20 & 30
\end{pmatrix}
$$
Step9:  The single row vector a is duplicated for as many rows as there are in b!  We can use this to calculate the squared distance between a center and every sample.
Step10:  Let's define $J$ from the book, which is a performance measure being minimized by k-means. It is defined on page 424 as
$$
J = \sum_{n=1}^N \sum_{k=1}^K r_{nk} ||\mathbf{x}_n - \mathbf{\mu}_k||^2
$$
where $N$ is the number of samples, $K$ is the number of cluster centers, $\mathbf{x}_n$ is the $n^{th}$ sample and $\mathbf{\mu}_k$ is the $k^{th}$ center, each being an element of $\mathbf{R}^p$ where $p$ is the dimensionality of the data.
The sums can be computed using python for loops, but for loops are much slower than matrix operations in python, as the following cells show.
Step11:  So, the matrix form is 10 times faster!
Now, back to our problem.  How can we use matrix operations to calculate the squared distance between two centers and, say, five data samples?  Let's say both are two-dimensional.
Step12:  This will be a little weird, and hard to understand, but by adding an empty dimension to the centers array, numpy broadcasting does all the work for us.
Step13:  Now we have a 2 x 5 array with the first row containing the squared distance from the first center to each of the five data samples, and the second row containing the squared distances from the second center to each of the five data samples.
Now we just have to find the smallest distance in each column and sum them up.
Step14:  Let's define a function named calcJ to do this calculation.
<ASSISTANT_TASK:>
Python Code:
!wget http://archive.ics.uci.edu/ml/machine-learning-databases/00299/StoneFlakes.dat
Explanation: K-Means Clustering
The UCI Machine Learning Repository is a "go-to" place for fascinating data sets often used in machine learning education and research.  
Let's check out the Stone Flakes data set.  Find and download the StoneFlakes.dat file by clicking here or by running wget like below.
End of explanation
!head StoneFlakes.dat
Explanation: Let's look at the first few lines.
End of explanation
import pandas
d = pandas.read_csv(open('StoneFlakes.dat'))
d[:5]
d = pandas.read_csv(open('StoneFlakes.dat'),sep=',')
d[:5]
Explanation: Read about the column names and the meaning of the ID values at the data set's web site.  
Notice that values are separated by commas, except for the first column.  Also notice that there are question marks where data is missing.  How can we read this?  Well, the usual answer is to "google" for the answer.  Try seaching for "read data set numpy"
numpy includes functions for reading csv files.  However, the pandas package offers more powerful functions for parsing data.  Let's try its read_csv function.
End of explanation
! tr -s ' ' ',' < StoneFlakes.dat > StoneFlakes2.dat
! head StoneFlakes2.dat
d = pandas.read_csv(open('StoneFlakes2.dat'))
d[:5]
d = pandas.read_csv(open('StoneFlakes2.dat'),na_values='?')
d[:5]
d = pandas.read_csv(open('StoneFlakes2.dat'),na_values='?',error_bad_lines=False)
d[:5]
d[:5].isnull()
d[:5].isnull().any(axis=1)
d[:5].isnull().any(axis=1) == False
print(d.shape)
d = d[d.isnull().any(axis=1)==False]
print(d.shape)
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
data = d.iloc[:,1:].values
data.shape
data[:5,:]
Explanation: Let's just replace commas with spaces, using unix.  Read about the tr unix command at Linux TR Command Examples
End of explanation
plt.plot(data);
Explanation: To see this data, let's try plotting each column as a separate curve on the same axes.
End of explanation
np.random.choice(range(data.shape[0]),2, replace=False) # data.shape[0] is number of rows, or samples
np.random.choice(range(data.shape[0]),2, replace=False)
centersIndex = np.random.choice(range(data.shape[0]),2, replace=False)
centersIndex
centers = data[centersIndex,:]
centers
Explanation: Each sample has 8 attribues, so each sample is a point in 8-dimensional space.  I wonder how well the samples "clump" in those 8 dimensions.  Let's try clustering them with the k-means algorithm.
First, let's try to find two clusters, so $k=2$.  We must initialize the two means of the two clusters.  Let's just pick two samples at random.
End of explanation
a = np.array([1,2,3])
b = np.array([10,20,30])
a, b
a-b
Explanation: Now we must find all samples that are closest to the first center, and those that are closest to the second sample.
Check out the wonders of numpy broadcasting.
End of explanation
np.resize(a,(3,3))
np.resize(b, (3,3))
np.resize(a,(3,3)).T
np.resize(a,(3,3)).T - np.resize(b,(3,3))
Explanation: But what if we want to subtract every element of a with every element of b?
End of explanation
a[:,np.newaxis]
a.reshape((-1,1))
a.shape, a[:,np.newaxis].shape
a[:,np.newaxis] - b
a = np.array([1,2,3])
b = np.array([[10,20,30],[40,50,60]])
print(a)
print(b)
b-a
Explanation: However, we can ask numpy to do this duplication for us if we reshape a to be a column vector and leave b as a row vector.
$$ \begin{pmatrix}
1\
2\
3
\end{pmatrix}
-
\begin{pmatrix}
10 & 20 & 30
\end{pmatrix}
\;\; = \;\;
\begin{pmatrix}
1 & 1 & 1\
2 & 2 & 2\
3 & 3 & 3
\end{pmatrix}
-
\begin{pmatrix}
10 & 20 & 30\
10 & 20 & 30\
10 & 20 & 30
\end{pmatrix}
$$
End of explanation
centers[0,:]
np.sum((centers[0,:] - data)**2, axis=1)
np.sum((centers[1,:] - data)**2, axis=1) > np.sum((centers[0,:] - data)**2, axis=1)
centers
centers[:,np.newaxis,:].shape, data.shape
(centers[:,np.newaxis,:] - data).shape
np.sum((centers[:,np.newaxis,:] - data)**2, axis=2).shape
np.argmin(np.sum((centers[:,np.newaxis,:] - data)**2, axis=2), axis=0)
cluster = np.argmin(np.sum((centers[:,np.newaxis,:] - data)**2, axis=2), axis=0)
cluster
data[cluster==0,:].mean(axis=0)
data[cluster==1,:].mean(axis=0)
k = 2
for i in range(k):
    centers[i,:] = data[cluster==i,:].mean(axis=0)
centers
def kmeans(data, k = 2, n = 5):
    # Initial centers
    centers = data[np.random.choice(range(data.shape[0]),k, replace=False), :]
    # Repeat n times
    for iteration in range(n):
        # Which center is each sample closest to?
        closest = np.argmin(np.sum((centers[:,np.newaxis,:] - data)**2, axis=2), axis=0)
        # Update cluster centers
        for i in range(k):
            centers[i,:] = data[closest==i,:].mean(axis=0)
    return centers
kmeans(data,2)
kmeans(data,2)
Explanation: The single row vector a is duplicated for as many rows as there are in b!  We can use this to calculate the squared distance between a center and every sample.
End of explanation
a = np.linspace(0,10,30).reshape(3,10)
a
b = np.arange(30).reshape(3,10)
b
result = np.zeros((3,10))
for i in range(3):
    for j in range(10):
        result[i,j] = a[i,j] + b[i,j]
result
%%timeit
result = np.zeros((3,10))
for i in range(3):
    for j in range(10):
        result[i,j] = a[i,j] + b[i,j]
result = a + b
result
%%timeit
result = a + b
Explanation: Let's define $J$ from the book, which is a performance measure being minimized by k-means. It is defined on page 424 as
$$
J = \sum_{n=1}^N \sum_{k=1}^K r_{nk} ||\mathbf{x}_n - \mathbf{\mu}_k||^2
$$
where $N$ is the number of samples, $K$ is the number of cluster centers, $\mathbf{x}_n$ is the $n^{th}$ sample and $\mathbf{\mu}_k$ is the $k^{th}$ center, each being an element of $\mathbf{R}^p$ where $p$ is the dimensionality of the data.
The sums can be computed using python for loops, but for loops are much slower than matrix operations in python, as the following cells show.
End of explanation
centers = np.array([[1,2],[5,4]])
centers
data = np.array([[3,2],[4,6],[7,3],[4,6],[1,8]])
data
Explanation: So, the matrix form is 10 times faster!
Now, back to our problem.  How can we use matrix operations to calculate the squared distance between two centers and, say, five data samples?  Let's say both are two-dimensional.
End of explanation
centers[:,np.newaxis,:]
centers[:,np.newaxis,:].shape
data.shape
diffsq = (centers[:,np.newaxis,:] - data)**2
diffsq
diffsq.shape
np.sum(diffsq,axis=2)
Explanation: This will be a little weird, and hard to understand, but by adding an empty dimension to the centers array, numpy broadcasting does all the work for us.
End of explanation
np.min(np.sum(diffsq,axis=2), axis=0)
np.sum(np.min(np.sum(diffsq,axis=2), axis=0))
Explanation: Now we have a 2 x 5 array with the first row containing the squared distance from the first center to each of the five data samples, and the second row containing the squared distances from the second center to each of the five data samples.
Now we just have to find the smallest distance in each column and sum them up.
End of explanation
def calcJ(data,centers):
    diffsq = (centers[:,np.newaxis,:] - data)**2
    return np.sum(np.min(np.sum(diffsq,axis=2), axis=0))
calcJ(data,centers)
def kmeans(data, k = 2, n = 5):
    # Initialize centers and list J to track performance metric
    centers = data[np.random.choice(range(data.shape[0]),k,replace=False), :]
    J = []
    
    # Repeat n times
    for iteration in range(n):
        
        # Which center is each sample closest to?
        sqdistances = np.sum((centers[:,np.newaxis,:] - data)**2, axis=2)
        closest = np.argmin(sqdistances, axis=0)
        
        # Calculate J and append to list J
        J.append(calcJ(data,centers))
        
        # Update cluster centers
        for i in range(k):
            centers[i,:] = data[closest==i,:].mean(axis=0)
            
    # Calculate J one final time and return results
    J.append(calcJ(data,centers))
    return centers,J,closest
centers,J,closest = kmeans(data,2)
J
plt.plot(J);
centers,J,closest = kmeans(data,2,10)
plt.plot(J);
centers,J,closest = kmeans(data,3,10)
plt.plot(J);
small = np.array([[8,7],[7,6.6],[9.2,8.3],[6.8,9.2], [1.2,3.2],[4.8,2.3],[3.4,3.2],[3.2,5.6],[1,4],[2,2.2]])
plt.scatter(small[:,0],small[:,1]);
c,J,closest = kmeans(small,2,n=2)
c
closest
plt.scatter(small[:,0], small[:,1], s=80, c=closest, alpha=0.5);
plt.scatter(c[:,0],c[:,1],s=80,c="green",alpha=0.5);
c,J,closest = kmeans(small,2,n=2)
plt.scatter(small[:,0], small[:,1], s=80, c=closest, alpha=0.5);
plt.scatter(c[:,0],c[:,1],s=80,c="green",alpha=0.5);
J
import gzip
import pickle
with gzip.open('mnist.pkl.gz', 'rb') as f:
    train_set, valid_set, test_set = pickle.load(f, encoding='latin1')
    # zero = train_set[0][1,:].reshape((28,28,1))
    # one = train_set[0][3,:].reshape((28,28,1))
    # two = train_set[0][5,:].reshape((28,28,1))
    # four = train_set[0][20,:].reshape((28,28,1))
X = train_set[0]
T = train_set[1].reshape((-1,1))
Xtest = test_set[0]
Ttest = test_set[1].reshape((-1,1))
X.shape, T.shape, Xtest.shape, Ttest.shape
c,J,closest = kmeans(X, k=10, n=20)
plt.plot(J)
c.shape
for i in range(10):
    plt.subplot(2,5,i+1)
    plt.imshow(-c[i,:].reshape((28,28)), interpolation='nearest', cmap='gray')
    plt.axis('off')
c,J,closest = kmeans(X, k=10, n=20)
plt.plot(J)
plt.figure()
for i in range(10):
    plt.subplot(2,5,i+1)
    plt.imshow(-c[i,:].reshape((28,28)), interpolation='nearest', cmap='gray')
    plt.axis('off')
c,J,closest = kmeans(X, k=20, n=20)
plt.plot(J)
plt.figure()
for i in range(20):
    plt.subplot(4,5,i+1)
    plt.imshow(-c[i,:].reshape((28,28)), interpolation='nearest', cmap='gray')
    plt.axis('off')
Explanation: Let's define a function named calcJ to do this calculation.
End of explanation
<END_TASK>