# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Large-scale Distributed $k$-means Algorithm
#
# In this Notebook, we'll focus on the development of a simple distributed algorithm. As for the Notebook on SGD, we focus on iterative algorithms, which eventually converge to a desired solution.
#
# In what follows, we'll proceed with the following steps:
#
# * We first introduce formally the $k$-means algorithm
# * Then we focus on a serial implementation. To do this, we'll first generate some data using scikit. In passing, we'll also use the $k$-means implementation in scikit to have a baseline to compare against.
# * Subsequently, we will focus on some important considerations and improvements to the serial implementation of $k$-means.
# * At this point, we'll design our distributed version of the $k$-means algorithm using pyspark, and re-implement the enhancements we designed for the serial version 
#
# #### References:
# * https://en.wikipedia.org/wiki/K-means_clustering
# * http://nbviewer.jupyter.org/github/ethen8181/machine-learning/blob/master/clustering/kmeans.ipynb
# * https://apache.googlesource.com/spark/+/master/examples/src/main/python/kmeans.py
# * https://github.com/castanan/w2v/blob/master/ml-scripts/w2vAndKmeans.py
#

# #### Preamble code

# %matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.cluster import KMeans
from sklearn.metrics import pairwise_distances
from sklearn.datasets.samples_generator import make_blobs
from math import sqrt
import random as random
import time
import datetime as datetime
import warnings
from sklearn import cluster, datasets, mixture
from sklearn.neighbors import kneighbors_graph
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice

# # Preliminaries: the $k$-means algorithm
#
# $k$-means clustering aims to partition $n$ $d-$dimensional observations into $k$ clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. This results in a partitioning of the data space into Voronoi cells.
# The problem is computationally difficult (NP-hard); however, there are efficient heuristic algorithms that are commonly employed and converge quickly to a local optimum. In this Notebook we'll focus on one of them in particular: the Lloyd algorithm.
#
# The $k$-means problem can be formalized as follows. Given a set of observations $(x_1, x_2, \cdots, x_n)$, where each observation is a $d$-dimensional real vector, $k$-means clustering aims to partition the $n$ observations into $k \leq n$ sets $S = \{S_1, S_2, \cdots, S_k\}$ so as to minimize the within-cluster sum of squares (WCSS) (i.e. variance). 
# The objective is to find:
#
# $$
# \arg \min_S \sum_{i=1}^{k} \sum_{x \in S_i} || \boldsymbol{x} - \boldsymbol{\mu_i} ||^2
# $$

# The most common algorithm uses an iterative refinement technique. Given an initial set of $k$ centroids $\mu_1^{(1)}, \cdots, \mu_k^{(1)}$ , the algorithm proceeds by alternating between two steps: in the **assignment step**, observations are associated to the closest **centroid**, in terms of squared Euclidean distance; in the **update step** new centroids are computed based on the new points associated to each centroid. Note: $\mu_i^{(t)}$ stands for the $i$-th centroid as of the $t$-th iteration. So $\mu_1^{(1)}$ is the centroid 1 at iteration 1.

# ## Algorithm pseudo-code
# It is important to work on a principled approach to the design of large-scale algorithms, and this starts with using good data structures and scientific libraries, such as ```numpy``` and ```scipy```. In particular, we will focus on the use of ```numpy``` arrays, which come with efficient methods for array operations. A pseudo-code for the $k$-means algorithm is specified below:
#
# ```python
# def kmeans(X, k, maxiter, seed=None):
#     """
#     specify the number of clusters k and
#     the maximum iteration to run the algorithm
#     """
#
#     # randomly choose k data points as initial centroids
#     centroids = X[rand_indices]
#     
#     for itr in range(maxiter):
#         # ---------------
#         # ASSIGNMENT STEP
#         # ---------------
#         # compute the distance matrix between each data point and the set of centroids
#         distance_matrix = # row Index = data point Index; col Index = centroid Index; value=distance
#         # assign each data point to the closest centroid
#         cluster_assignment = # array Index = data point Index; array value = closest centroid Index
#         
#         # UPDATE STEP
#         # select all data points that belong to cluster i and compute
#         # the mean of these data points (each feature individually)
#         # this will be our new cluster centroids
#         new_centroids = ...
#         
#         # STOP CONDITION
#         # if centroids == new_centroids => stop
#  
#     
# ```

# # Synthetic data generation: working in bi-dimensional spaces
#
# Next, we use sklearn to generate some synthetic data to test our algorithm.

X, y = make_blobs(n_samples = 1000, centers = 4,
                  random_state = 0, cluster_std = 0.4)

# +
# change default figure and font size
plt.rcParams['figure.figsize'] = 8, 6 
plt.rcParams['font.size'] = 12

# scatter plot
plt.scatter(X[:, 0], X[:, 1], s = 50)
plt.ylim(-2, 10)

plt.xlim(-6, 6)
plt.gca().set_aspect('equal')
plt.show()


# -

# <div class="alert alert-block alert-info">
# <strong>Question 1.</strong> Implement your own version of k-means, as a serial algorithm.
#
# Follow the guidelines below:
#
# <ul>
#     <li> Define a function to perform k-means clustering. The function should accept as inputs: the training data x, the number of clusters k, and the iteration budget you allocate to the algorithm. Additional arguments might include the use of a random seed to initialize centroids.</li>
#     <li>The function should output the centroids, and the cluster assignment, that is, to which centroid each data point is assigned to</li>
#     <li> Optionally, keep track of the position of the centroids, for each iteration.</li>
# </ul>
# </div>

def kmeans(data, k, maxiter, seed=None):
    """
    specify the number of clusters k and
    the maximum iteration to run the algorithm
    """
    
    # randomly choose k data points as initial centroids
    if(seed is not None):
        np.random.seed(seed)
        
    rand_indices = np.random.randint(len(data), size=k)
    centroids = np.array(data[rand_indices])
    
    centroids_history = []

    for itr in range(maxiter):
        centroids_history.append(centroids)
        # ---------------
        # ASSIGNMENT STEP
        # ---------------
        # compute the distance matrix between each data point and the centroids
        #the metric used to compute the distance between a centroic and a data point is sqeuclidean: squared Euclidean distance
        distance_matrix = pairwise_distances(data, centroids, metric='sqeuclidean')
        
        # assign each data point to the closest centroid
        cluster_assignment = np.argmin(distance_matrix, axis = 1)
        #print(cluster_assignment)
        
        # UPDATE STEP
        # select all data points that belong to cluster i and compute
        # the mean of these data points (each feature individually)
        # this will be our new cluster centroids
        new_centroids = np.empty([k, 2])
        for i in range(k):
            new_centroids[i] = (np.mean(data[cluster_assignment == i], axis=0))

        # STOP CONDITION
        # stop if centroids are equals for two iteration
        if(np.array_equal(centroids, new_centroids)):
            break
        #Set the centroids to the new updated ones
        centroids = new_centroids
    #print(iter,itr)    
    return (centroids, cluster_assignment, np.array(centroids_history))


# <div class="alert alert-block alert-info">
# Once the ```kmeans``` function is defined, you can generate input data according to the cell above, that uses scikitlearn.
# The output of your cell should contain the following information:
# <ul>
#     <li> Print the number of data points that belong to each cluster</li>
#     <li> Plot the clustered data points:</li>
#     <ul>
#         <li>Using different colors for each cluster</li>
#         <li>Plot the centroid positions for each cluster</li>
#     </ul>
# </ul>
#  </div>   

def plot_clusters(dataset, centroids, cluster_assignment, centroids_history=None):
    import matplotlib.cm as cm
    
    fig = plt.figure(figsize=(15,8))
    #Initial data sets
    ax1 = fig.add_subplot(131)
    ax1.scatter(dataset[:, 0], dataset[:, 1], s = 40)
    if(centroids_history is not None):
        ax1.scatter(centroids_history[0,:,0], centroids_history[0,:,1], s = 80, c='black')
    plt.title("Initial Data generation")
    #Clusters 
    ax2 = fig.add_subplot(132)
    ax2.scatter(dataset[:, 0], dataset[:, 1], s=40, c = cluster_assignment, cmap=cm.rainbow, edgecolor='')
    ax2.scatter(centroids[:, 0], centroids[:, 1], s = 80, c='black')
    plt.title("K-means Clusters")
    #Plot the evolution of the centroids 
    if(centroids_history is not None):
        ax3 = fig.add_subplot(133)
        ax3.scatter(dataset[:, 0], dataset[:, 1], s=50, c = cluster_assignment, cmap=cm.rainbow, edgecolor='')
        for i in range(len(centroids)):
            ax3.plot(centroids_history[:,i,0], centroids_history[:,i,1], '->', color='black', linewidth=2.5)
            plt.title("Centroids Evolution")
    plt.gca().set_aspect('equal')
    plt.show()


# +
k = 4
kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history = kmeans(X, k, 300)

### OUTPUT ###
print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)
# -

# <div class="alert alert-block alert-info">
# <strong>Question 2.</strong> Use the built-in k-means implementation in sklearn and determine centroids and clusters.
#
# Follow the guidelines below:
#
# <ul>
#     <li> Use the KMeans algorithm from sklearn</li>
#     <li> Use the fit_predict method to cluster data</li>
#     <li> Use the cluster_centers_ method to retrieve centroids</li>
# </ul>
#
# The output of your cell should contain the following information:
# <ul>
#     <li> Plot the clustered data points, using the same code your have produced for Question.1</li>
#     <ul>
#         <li>Using different colors for each cluster</li>
#         <li>Plot the centroid positions for each cluster</li>
#     </ul>
# </ul>
# </div>

# +
k = 4
#Built-in k-means implementation in sklearn
kmeans_sklearn = KMeans(k, 
                init='random', 
                n_init=10, 
                max_iter=500, 
                tol=0.0001, 
                precompute_distances='auto', 
                verbose=0, 
                random_state=None, 
                copy_x=True,
                n_jobs=1)
#Cluster data with fit predict
kmeans_pred = kmeans_sklearn.fit_predict(X)
#Retrieve the centroids
kmeans_centroids = kmeans_sklearn.cluster_centers_

print("Number of points in each Cluster:", np.bincount(kmeans_pred))

plot_clusters(X, kmeans_centroids, kmeans_pred)

# -

# <div class="alert alert-block alert-info">
# <strong>Question 3.</strong> Use the sklearn dataset API to generate alternative synthetic data to test your k-means algorithm implementation.
#
# Follow the guidelines from this document: http://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html#sphx-glr-auto-examples-cluster-plot-cluster-comparison-py
#
#
# The output of your cell should contain the following information:
# <ul>
#     <li> Plot the new synthetic dataset you generated
#     <li> Plot the clustered data points, using the same code your have produced for Question.1</li>
#     <ul>
#         <li>Using different colors for each cluster</li>
#         <li>Plot the centroid positions for each cluster</li>
#     </ul>
# </ul>
# </div>

# +
n_samples = 1500

X, y = make_blobs(n_samples = 1000, centers = 4,
                  random_state = 0, cluster_std = 0.4)

#noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
k = 4
kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history = kmeans(X, k, 300)

### OUTPUT ###
print("Blobs")
print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)


X, y = datasets.make_circles(n_samples = 1000, factor=.5, noise=.05)

#noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
k = 4
kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history = kmeans(X, k, 300)

### OUTPUT ###
print("Circles")
print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)

X, y = datasets.make_moons(n_samples = 1000, noise=.05)

#noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
k = 4
kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history = kmeans(X, k, 300)

### OUTPUT ###
print("Moons")
print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)


# -

# ***

# # A simplified analysis of algorithm convergence
#
# One well-known weakness of $k$-means is that the algorithm's performance is closely tied with the randomly generated initial centroids' quality. If the algorithm starts with a set of bad inital centers, it will get stuck in a local minimum.
#
# Instead of taking a formal approach to study the convergence of $k$-means, let's study it with an experimental approach. One thing we can do is to build a measure of clustering quality: intuitively, a good clustering result should produce clusters in which data points should be very close to their centroids, and very far from other centroids. In this Notebook, we'll look at a metric called the **total within Sum of Squares**, which is sometimes referred ot as heterogeneity. Mathematically, we define heterogeneity as:
#
# $$
# \mathcal{H} = \sum_{j=1}^{k} \sum_{i: z_i=j} || \boldsymbol{x_i} - \boldsymbol{\mu_j}||_{2}^{2}
# $$
#
# Where $k$ denotes the total number of clusters, $x_i$ is the $i$-th data point, $\mu_j$ is the $j$-th centroid, and $|| \cdot ||_{2}^{2}$ denotes the squared L2 norm (Euclidean distance) between the two vectors.

# <div class="alert alert-block alert-info">
# <strong>Question 4.</strong> Modify your own version of k-means, to compute heterogeneity as defined above.
#
# Follow the guidelines below:
#
# <ul>
#     <li>Use the same method template you used in Question 1</li>
#     <li>Add the code required to compute heterogeneity</li>
#     <li>The function should return, in addition to the same return values as for the baseline version, the computed heterogeneity
# </ul>
# </div>

def heter_kmeans(data, k, maxiter, seed):
        d ={}
        kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history = kmeans(X, k, maxiter,seed)
        heterogeneity = 0
        for i in range(k):
            # note that pairwise_distance only accepts 2d-array
            cluster_data = X[kmeans_cluster_assignment == i]
            distances = pairwise_distances(cluster_data, [kmeans_centroids[i]], metric = 'euclidean')
            heterogeneity += np.sum(distances ** 2)
        d["seed"] = seed
        d["heterogeneity"] = heterogeneity
        return (kmeans_centroids, kmeans_cluster_assignment, d, kmeans_centroids_history)


# <div class="alert alert-block alert-info">
# <strong>Question 5.</strong> Using the modified k-means method you designed, study algorithm convergence as a function of heterogeneity.
#
# Follow the guidelines below:
#
# <ul>
#     <li>Run the modified k-means for at least 5 different initial seed values</li>
#     <li>Prepare a dictionary data structure containing: key = random seed, value = heterogeneity</li>
#     <li>Print seed, heterogeneity values</li>
# </ul>
#
# Add your personal comment about the convergence properties of the $k$-means algorithm.
#
# </div>

X, y = make_blobs(n_samples = 1000, centers = 4,
                  random_state = 0, cluster_std = 0.4)
k = 4
for s in [10,100,500,1000,5000,10000,100000]:
    kmeans_centroids, kmeans_cluster_assignment, d, kmeans_centroids_history = heter_kmeans(X, k, 300, s)

    ### OUTPUT ###
    print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
    print("Output:", d)
    print("---------------------------------")
    #plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)


# <div class="alert alert-success">
# We can see clearly that for some seed values, the k-means algorithms gets stuck at a  at a local minimum where the heterogenity (squared distance sums) are around 304,14 or 1263 regardless of what the initial seed is. This is happening because the dataset is already well clustered to begin with. As the dataset size is not big enough, it does not have enough time to converge to a better stage before clustering ends(because no new update can be done). Note that, the seed affects a lot the distribution of random selection of points and centroids, so it in fact does affect the clustering. To remove this issue, we are looking for solutions that will converge to the best stage regardless of initial centroid selections or seed being different on each time. We would like to remove the randomness out of the equation while still benefit from it to prove the dataset it being clustered properly on different initial conditions.
# </div>

# ***

# # A technique for a smart centroid initialization: $k$-means++
#
# One effective way to produce good initial centroids to feed to $k$-means is to proceed as follows: instead of randomly generating initial centroids, we will try to spread them out in our $d$-dimensional space, such that they are not "too close" to eachother. If you are interested in the details of this technique, you should refer to the link to the original research paper below: in summary, the $k$-means++ technique allows to improve the quality of the local optima in which $k$-means settles, and to reduce the average runtime of the algorithm.
#
# k-means++: the advantages of careful seeding, by <NAME> and <NAME> (Stanford University) https://dl.acm.org/citation.cfm?id=1283494
#
# A simplified workflow of the $k$-means++ approach is as follows:
#
#
# * Choose a data point at random from the dataset, this serves as the first centroid
# * Compute the squared euclidean distance of all other data points to the randomly chosen first centroid
# * To generate the next centroid, each data point is chosen with the probability (weight) of its squared distance to the chosen center in the current round, divided by the the total squared distance (this is just a normalization to make sure the probability adds up to 1). In other words, a new centroid should be as far as possible from the other centroids
# * Next, recompute the probability (weight) of each data point as the minimum of the distance between it and all the centers that are already generated (e.g. for the second iteration, compare the data point's distance between the first and second center and choose the smaller one)
# * Repeat step 3 and 4 until we have $k$ initial centroids to feed to the $k$-means algorithm
#

# <div class="alert alert-block alert-info">
# <strong>Question 6.</strong> Modify your own version of k-means, to introduce the smart initialization technique described above. Don't forget to keep track of heterogeneity as well! The whole point is to measure if k-means++ really improves on this metric.
#
# <br>
# <br>
#
# Follow the guidelines below:
#
# <ul>
#     <li>Use the same method template you used in Question 1</li>
#     <li>Add the code required to compute the initial clusters according to k-means++</li>
#     <li>Add the code required to compute heterogeneity</li>
#     <li>The function should return, in addition to the same return values as for the baseline version, the computed heterogeneity
# </ul>
# </div>

def kmeans_pp(data, k, maxiter, seed=None):
    """
    specify the number of clusters k and
    the maximum iteration to run the algorithm
    """
    
    # randomly choose k data points as initial centroids
    if(seed is not None):
        np.random.seed(seed)
        
    rand_indices = np.random.randint(len(data), size=k)
    centroids = np.array(data[rand_indices])
    
    centroids_history = []
    
    # compute distances from the first centroid chosen to all the other data points
    distances = pairwise_distances(X, [centroids[0]], metric = 'sqeuclidean').flatten()
    
    for i in range(1, k):
            # choose the next centroid, the probability for each data point to be chosen
            # is directly proportional to its squared distance from the nearest centroid
            prob = distances ** 2
            rand_index = np.random.choice(len(data), size = 1, p = prob / np.sum(prob))
            centroids[i] = X[rand_index]
            if i == k - 1:
                break
            
            # if we still need another cluster,
            # compute distances from the centroids to all data points
            # and update the squared distance as the minimum distance to all centroid
            distances_new = pairwise_distances(X, [centroids[i]], metric = 'euclidean').flatten()
            distances = np.min(np.vstack((distances, distances_new)), axis = 0)
    

    for itr in range(maxiter):
        centroids_history.append(centroids)
        # ---------------
        # ASSIGNMENT STEP
        # ---------------
        # compute the distance matrix between each data point and the centroids
        #the metric used to compute the distance between a centroic and a data point is sqeuclidean: squared Euclidean distance
        distance_matrix = pairwise_distances(data, centroids, metric='sqeuclidean')
        
        # assign each data point to the closest centroid
        cluster_assignment = np.argmin(distance_matrix, axis = 1)
        
        
        # UPDATE STEP
        # select all data points that belong to cluster i and compute
        # the mean of these data points (each feature individually)
        # this will be our new cluster centroids
        new_centroids = np.empty([k, 2])
        for i in range(k):
            new_centroids[i] = (np.mean(data[cluster_assignment == i], axis=0))

        # STOP CONDITION
        # stop if centroids are equals for two iteration
        if(np.array_equal(centroids, new_centroids)):
            break
        #Set the centroids to the new updated ones
        centroids = new_centroids
        
        d ={}
        heterogeneity = 0
        for i in range(k):
            # note that pairwise_distance only accepts 2d-array
            cluster_data = X[cluster_assignment == i]
            distances = pairwise_distances(cluster_data, [centroids[i]], metric = 'euclidean')
            heterogeneity += np.sum(distances ** 2)
        d["seed"] = seed
        d["heterogeneity"] = heterogeneity
        
    return (centroids, cluster_assignment,d, np.array(centroids_history))

# <div class="alert alert-block alert-info">
# <strong>Question 7.</strong> Similarly to question 5, using the modified k-means method you designed, study algorithm convergence as a function of heterogeneity.
# <br>
# <br>
#
#
# Follow the guidelines below:
#
# <ul>
#     <li>Run the modified k-means for at least 5 different initial seed values</li>
#     <li>Prepare a dictionary data structure containing: key = random seed, value = heterogeneity</li>
#     <li>Print seed, heterogeneity values</li>
# </ul>
#
# <br>
#
# One additional question to answer is the following: print the average heterogeneity for the baseline k-means algorithm, and the average heterogeneity when using the k-means++ initialization. Compare and comment with your own words.
# </div>

# +
print("The old heter_kmeans method:")
kmeans_centroids, kmeans_cluster_assignment,d, kmeans_centroids_history = heter_kmeans(X, 4, 300, 12345)
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)

print("The new kmeans_pp method:")
kmeans_centroids, kmeans_cluster_assignment, d, kmeans_centroids_history = kmeans_pp(X, 4, 300, 12345)
plot_clusters(X, kmeans_centroids, kmeans_cluster_assignment, kmeans_centroids_history)
# -

# We can see that with kmeans_pp the initial centroids as much better chosen, so we don't need many iterations to reach the best clusters.

k = 4
het_kpp = []
het_k = []
for s in [10,100,500,1000,5000,10000,100000]:
    _, kmeans_cluster_assignment, d_kpp, _ = kmeans_pp(X, k, 300, s)
    het_kpp.append(d_kpp["heterogeneity"])
    ### OUTPUT ###
    print("Number of points in each Cluster:", np.bincount(kmeans_cluster_assignment))
    print("Output:", d_kpp)
    
    _, kmeans_cluster_assignment, d_k, _ = heter_kmeans(X, k, 300, s)
    het_k.append(d_k["heterogeneity"])
    print("---------------------------------")
print('kmeans average: ',np.mean( np.array(het_k )))
print('kmeanspp average: ',np.mean( np.array(het_kpp )))

# <div class="alert alert-success">
# Please notice the kmeans++ algorithms in average actually performed much better than normal kmeans we had on Q1. The reason is simple, we have introduced a logic that allows kmeans++ to make more beneficial, logical initial centroids to begin the algorithm rather than a blind randomness. It works same as the probability wheel logic we learnt in MALIS for genetic algorithm's parent selection, the points on more distance are more likely to be chosen to become the centroids. Also keep in mind that there might be a seed value that chooses terrible initial points because the randomness actually chooses points with the least possibility to be chosen. We can also eliminate this issue by using different logics, such as logarithmic probability selection, or choosing a set of confidence interval and removing the points under that limit so they cannot be chosen as centroids to begin with.
# </div>

# ***

# # Determining the value of k: a simple and visual approach, called the Elbow method
#
# Another problem of $k$-means is that we have to specify the number of clusters $k$ before running the algorithm, which we often don't know a priori. 
# There are many different heuristics for choosing a suitable value for $k$, the simplest one being the **Elbow method**. Essentially, the idea is to run the $k$-means algorithm using different values of $k$ and plot the corresponding heterogeneity. This measure will decrease as the number of clusters increases, because each cluster will be smaller and tighter. 
# By visual inspection of the plot heterogeneity vs. $k$, we will (hopefully!) see that the curve flattens out at some value of $k$: this is what we call an "elbow", and we'll select the value of $k$ corresponding to the "elbow" position.
#

# <div class="alert alert-block alert-info">
# <strong>Question 8.</strong> Using the k-means algorithm you implemented, including the smart initialization technique, collect into an array the value of heterogeneity as a function of the number of clusters k, where k is to be selected in the range [2,10].
# <br>
# <br>
#
#
# Your output cell should contain:
#
# <ul>
#     <li>The plot of heterogeneity vs. k</li>
#     <li>A discussion on your visual inspection of the curve, together with a justification for an appropriate choice of the value k</li>
# </ul>
# </div>

# +
het = []
def elbow_method(data,maxiter):
    for k in list(range(2,11)):
        start_time = time.clock()
        kmeans_centroids, kmeans_cluster_assignment, d, kmeans_centroids_history = kmeans_pp(data, k, maxiter, seed=None)
        het.append(d["heterogeneity"])
        end_time = time.clock()
        print("For k=",k," Computation Time: ", str(end_time-start_time))
    return (het)


heterogeneity_values = elbow_method(X, 300)
fig = plt.figure(figsize=(18,9))
plt.plot(list(range(2,11)), heterogeneity_values, marker = 'o')
plt.xlabel('Number of clusters')
plt.ylabel('heterogeneities')
plt.show()


# -

# <div class="alert alert-success">
# First thing first, there was not much difference on computation time when X dataset was only 1000 values. Therefore, we pushed it up to 100K in order to see the actual computational time difference.
#
# Now please notice the difference of computation time when we increase the number of clusters. Having more k points might lower the totoal distance of points to centroids but comes with a great computational complexity and possible over clustering issues as the data might have less than that many clusters initially. Similarly clustering on very few k will result in very very high total distance of points causing terrible results in the end.
#
# We have concluded that 3 or 4 clusters would be the best approach here. As the elbow being on k=4, making it the best solution as the heterogenety decrease a lot between 3 to 4, but not that much after k=4.
# </div>

# ***

# # Distributed $k$-means with PySpark
#
# We're now ready to work on a distributed implmentation of the $k$-means algorithm, using the PySpark API.
#
# By now, you should be rather familiar with the $k$-means algorithm, which means we can focus on its parallel and distributed design.

# ## Distributed algorithm pseudo code
#
# The basic idea of distributed $k$-means is as follows: data points to be clustered should be stored as a distributed dataset, namely a RDD. As in the Notebook on distributed SGD, we will take a shortcut and avoid using HDFS RDDs: rather, we'll use sklearn to generate the data points, similary to the serial version of the algorithms, then use the ```parallelize()``` method to create an RDD, and determine the number of partitions.
#
# We also need to manipulate the centroids array: indeed, all machines should hold a copy of the centroid vector, such that they can proceed independently and in parallel in the first phase of the $k$-means algorithm, that is the **assignment step**. Specifically, every worker has a set of data points, and it will use a copy of the centroid vector to compute cluster assignement: we compute the distance between each data point and each centroid, to assign data points to their closest centroid.
#
# Once the assignement step is done, we need to recompute new centroids based on the assignement, that is, we execute the **update step**. Clearly, we will need to **shuffle** data over the network such that we will have, for each current centroid, the list of all data points that have been assigned to it. If you think about it, this problem should be familiar!! This is very similar to what we do in the Word Count example. As such, you will need to make sure the output of the update step is cast to a ```<key, value>``` type, where the key corresponds to a centroid identifier, and the value contains the list of data points associated to that centroid. The framework will take care of the distributed group by operation, and organize data according to the semantic of our algorithm.
#
# **NOTE:** since we will (potentially) work on large dataset sizes, we don't want our algorithm to return the final assignement after convergence, for otherwise we would need to collect a large amount of data in the driver machine, which has a finite and somehow limited amount of RAM.
#
# The pseudo code of the algorithm you need to implement is as follows:
#
# ```python
# datapoints = # Use sklearn, as usual, and work on blobs
# centroids = # Random initialization
#
# for itr in range(maxiter): # This for loop is executed by the driver
#     bcCentroids = sc.broadcast(centroids) # Use broadcast variables
#     
#     closest = datapoints.mapPartition(assignement_step) # This should happen in parallel
#     
#     centroids = closest.reduceByKey(update_step_sum). \ # This should happen in parallel
#         map(update_step_mean). \ # This should happen in parallel
#         collect() # Here we collect new centroids in the driver
# ```
#
# As you can see from the pseudo code, you need to figure out how to implement the ```assignement_step``` function and the update_step function. For the latter, the pseudo code gives you a big hint! Remember what we discussed in class about computing the mean!

# <div class="alert alert-block alert-info">
# <strong>Question 9.</strong> Implement the distributed version of the k-means algorithm, following the guidelines in the pseudo code.
# <br>
# <br>
#
# Your output cell should contain:
#
# <ul>
#     <li>The value of the centroids once the algorithm converges</li>
#     <li>The total runtime of the distributed algorithm, in seconds</li>
#     <li>A visualization of the data points and the computed centroids</li>
# </ul>
# </div>

def kmeans_distributed(datapoints, k, NbrPartitions, maxiter, seed=None):

    def assignement_step(datapoints):

        for x_point in datapoints:
            # ---------------
            # ASSIGNMENT STEP
            # ---------------
            # compute the distance matrix between each data point and the centroids
            #the metric used to compute the distance between a centroic and a data point is sqeuclidean: squared Euclidean distance
            distance_matrix = pairwise_distances(x_point.reshape(1, -1), bcCentroids.value, metric='sqeuclidean')
            # assign each data point to the closest centroid
            point_assignment = np.argmin(distance_matrix)
            yield (point_assignment, [x_point])
    
    random.seed(seed)
    #partition data across the machines we have in the cluster. 
    rdd_data = sc.parallelize(datapoints, NbrPartitions)#.cache() # we have 2 workers in our Cluster
    #print("Number of partitions: {}".format(rdd_data.getNumPartitions()))
    # Random initialization
       
    rand_indices = np.random.randint(datapoints.shape[0], size=k)
    centroids = np.array(datapoints[rand_indices])
    #print(centroids)
    
    centroids_history = []

    for itr in range(maxiter): # This for loop is executed by the driver
        bcCentroids = sc.broadcast(centroids) # Use broadcast variables
        centroids_history.append(centroids)
        closest = rdd_data.mapPartitions(assignement_step) # This should happen in parallel
        # This should happen in parallel
        closest_reduced = closest.reduceByKey(lambda x,y: (x+y))
        new_centroids = closest_reduced.map(lambda x: np.mean(np.array(x[1]), axis=0)).collect() # Here we collect new centroids in the driver
        # Check if we are improving
        if np.all(np.array(new_centroids) == centroids):
            break
        centroids = new_centroids
        
    return (centroids, centroids_history, itr)


sc
# Use sklearn, as usual, and work on blobs
datapoints, y = make_blobs(n_samples = 1000, centers = 4,
                  random_state = 0, cluster_std = 0.4)

# +
start_time = time.clock()
centroids, centroids_history, iterations = kmeans_distributed(datapoints, 4, 6,100, seed=None)
end_time = time.clock()
print("For k=",4," Computation Time: ", str(end_time-start_time))
print("centroids:",centroids)
#final cluster asignement
distance_matrix = pairwise_distances(datapoints, centroids, metric='sqeuclidean')
# assign each data point to the closest centroid
cluster_assignment = np.argmin(distance_matrix, axis = 1)

plot_clusters(datapoints, np.array(centroids), cluster_assignment, np.array(centroids_history))
# -

# <div class="alert alert-block alert-info">
# <strong>Question 10.</strong> Answer the following questions:
# <br>
# <br>
# <ul>
#     <li>How many partitions did you use? Why?</li>
#     <li>What is the size of the dataset you generate? Did you cache the dataset? What's the RAM occupation?</li>
#     <li>What is the size of the shuffle data over the network? How does it compare to the dataset size?</li>
# </ul>
# </div>

#Impact of number of partitions on computation time 
#data_size = 1000
k=4
for NumberofPartitions in [2,4,6,8,10,12]:
        start_time = time.clock()
        centroids, centroids_history, itertations = kmeans_distributed(datapoints, k, NumberofPartitions,100, seed=None)
        end_time = time.clock()
        print("For  NumberofPartitions=",NumberofPartitions," Computation Time: ", str(end_time-start_time))

# <div class="alert alert-success">
# The cluster we have is composed of 2 workers each with 3 cores, each with 9 GB RAM to use. The best partitions that we can opt for is with having 6 partitions, to use the whole configuration that we have.
# It depends also on the data size, in this cas that data size size we have is relatively small, but with bigger datasets, we can opt for increasing the number of partitions, for better perfomance, but always with a number that is a multiple of 6.
#     
# As seen above, having a number of partions of 4 (2 for each worker), 6 ( 3 for each worker) gives the best computation time. But having a number of partition of 4, is not an optimal idea, because it won't use the whole configuration and it will add some overhead in the case of bigger datasets.
# </div>

# +
#Impact of data size WITHOUT CACHE on computation time 
k=4
NumberofPartitions = 6

for N_samples in [100,1000,10000,100000]:
        datapoints, y = make_blobs(n_samples = N_samples, centers = 4,random_state = 0, cluster_std = 0.4)
        start_time = time.clock()
        centroids, centroids_history, itertations = kmeans_distributed(datapoints, k, NumberofPartitions,100, seed=None)
        end_time = time.clock()
        print("For  N_sample=",N_samples," ,Data size= ",sys.getsizeof(datapoints)/1024,"KB, Computation Time: ", str(end_time-start_time))


# +
#Impact of data size WITH CACHE on computation time 
k=4
NumberofPartitions = 6

for N_samples in [100,1000,10000,100000]:
        datapoints, y = make_blobs(n_samples = N_samples, centers = 4,random_state = 0, cluster_std = 0.4)
        start_time = time.clock()
        centroids, centroids_history, itertations = kmeans_distributed(datapoints, k, NumberofPartitions,100, seed=None)
        end_time = time.clock()
        print("For  N_sample=",N_samples," ,Data size= ",sys.getsizeof(datapoints)/1024,"KB, Computation Time: ", str(end_time-start_time))

# -

# <div class="alert alert-success">
# the data size changes when we change the n-samples to generate data,we did not cache the dataset because,it is already into the memory (in fact we use the parallelize function to distribute the data). However, in order to avoid multiple reads from the disk, it would have been a good idea to cache the dataset in case it was already on the disk.
# The RAM occupation is 29 KB for each worker; it's possibile to see that on the Spark UI.
#
# With comparing results of distributed computation with and without cache, we can see clearly the effect of caching the reduce the computation time and the I/O from the disk, only when the datasets are big enough
# </div>

# <div class="alert alert-block alert-info">
# <strong>Question 11.</strong> Comparison between serial and distributed implementations.
# <br>
# <br>
# Given the dataset size you chose for your experiments, answer the following questions:
# </div>

# <li>Which is "faster", the serial or distributed implementation of k-means?</li>
#

# <div class="alert alert-success">
# For dataset with small number of examples, from 100 to 1000, we can see clearly from the test above that the serial algorithm is much faster than the distributed algorithm, and thus due the the smaal dataset, that can be read and manipulated quickly, with less computation, whereas for the distributed algorithm we have partitioned data for the workers, and we have I/O calls with computations on the workers.
# </div>
#

# <li>What is a dataset size for which the distributed implementation is clearly faster than the serial one?</li>
#

# <div class="alert alert-success">
# It is clear that for big dataset, starting from 100000, the distributed algorithm performs better that the serial one, in these cases we have dataset that are bigger than the a the available memory (RAM) and the dataset size does not fit anymore in a single machine memory. In this case, the serial algorithm will be facing memory management problems, whereas the distributed algorithm is made essentially to deal with these problems and to parallelise big amoun of data for computation.
#     </div>

# <li>What would be different in your code, should the input dataset reside on disk? Clearly, the input RDD would be reading from HDFS. Any other differences with respect to partitions?</li>

# <div class="alert alert-success">
# The main difference is that we would send our dataset to the workers through the parallelize function but we would have to read it through HDFS RDDs. In this case is it better to use caching in order to avoid to trigger multiple I/O operation from HDFS.
#
# The number of partition will be related to the number of blocks of the dataset. To guarantee a good load balancing, it may be necessary to do a repartitioning, by using the repartition() function.
#    </div> 

StarcoderdataJupyterScriptsDedupFiltered