# ---
# 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
# ---

# + [markdown] id="oo4hlNqu43zp"
# #Homework 6: Clustering (K-means, DP-means)
# by <NAME> and <NAME>
#
# **Due date**: March 24, Wednesday by 11:59pm
#
# **Late** due date: March 27, Saturday by 11:59pm
#
# This homework has a total of 100 points, with a 10-point bonus question at the end. **Online code along with library functions for K-means or DP-means are prohibited.**
#
# For all the functions you need to implement, feel free to change their signatures (inputs and returns) but do NOT change their names.

# + [markdown] id="jIr7WQmm43zt"
# ##**Question 1:**  K-means (35 points)
#
# Here you will apply K-means on simple Gaussian data. You are given some code to assist you in implementing this clustering method.

# + [markdown] id="8MVORa2A43zu"
# ###**a.** Create a dataset and plot it (10 points)
# Generate $3$ $2$D Gaussian clusters of data with the following means and covariances:
#
# $\boldsymbol{\mu}_0 = [2,2]^T, \boldsymbol{\mu}_1 = [-2,2]^T, \boldsymbol{\mu}_2 = [0,-2]^T$,
#
# $\Sigma_0 = [[0.02,0];[0,0.02]]$, $\Sigma_1 = [[0.2,0];[0,0.2]]$, $\Sigma_2 = [[0.05,0];[0,0.05]]$ 
#
# Let each cluster of data have $50$ points. Complete the scatter_by_category function and use it create a colored scatter plot of the generated Gaussian data.
#
# **Important:** Assume **row** vectors for feature vectors and **column** vectors for mean vectors for the **entire** homework 6. This is a convention that makes matrix multiplications and dot products easier.
#
# Hints: 
# 1. Use rng.multivariate_normal() to generate the dataset.
# 2. Use the imported seaborn library for plotting. Relevant documentation: https://seaborn.pydata.org/generated/seaborn.scatterplot.html

# + id="NGS7AyWBY0R2"
import numpy as np
import pandas as pd
# %matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import numexpr as ne
from scipy import sparse
from numba import njit
import time

# + colab={"base_uri": "https://localhost:8080/", "height": 337} id="iLqXr-qb43zu" outputId="3c5a9db8-9a0d-4669-cf9b-6e2a5956abf4"
# Do NOT change the seed in any problem, unless you are explicitly allowed to
# Use rng instances instead of using np.random directly, as recommended 
# by numpy manuals
rng = np.random.default_rng(0)

# Generate 2D Gaussian Data with 3 clusters of 3 different means
# WRITE CODE HERE:

# m is an optional argument for cluster means; when supplied, the function
# should plot them in red
# As mentioned, the means are column vectors, so m has shape (d, k), 
# where d is the number of features
def scatter_by_category(X, labels, m=None):
    plt.figure(figsize=(5,5))
    plt.xlim([-4, 3])
    plt.ylim([-3, 4])
    # WRITE CODE HERE:

scatter_by_category(DATA, DATA_labels)


# + [markdown] id="uTwm23dUdELr"
# ###**b.** K-means objective (5 points)
# Implement the K-means objective (within-cluster sum of squares) function: $\sum\limits_{i=1}^k Var(C_i) = \sum\limits_{i=1}^k \sum\limits_{\mathbf x \in C_i} ||\mathbf x-\mathbf{\mu_i}||_2^2$
#
# Note: we recommend that you either use njit or not use loops here. More note on njit can be found in question 3.b where a loop over the data is inevitable.

# + id="WJAufnXpdE0B"
@njit(fastmath=True,cache=True)
def WCSS(X, labels, m):
    # WRITE CODE HERE:


# + [markdown] id="SyA7SoGv43zw"
# ###**c.** Implement K-means (20 points)
# Implement K-means and apply it to the Gaussian data above. The code below contains comments of what needs to be coded up at a specific part of the code.  Please write your code below these comments where necessary. 
#
# Plot the created clusters per iteration in different colors. For convergence of the algorithm (stopping criterion), check whether the objective value becomes stationary (change less than a suitable threshold). Finally, print the final objective value.
#
# Notes: 
# 1. As always, try to use vectorized routines instead of loops whenever possible. This will be important in a later question.

# + colab={"base_uri": "https://localhost:8080/", "height": 1000} id="YEeEok_YFGBP" outputId="9673776c-1dee-4cbd-a9b3-4f0e9d9416d9"
## K-means implementation:
# X is the points as row vectors, m is the initial centers as column vectors
# maxT helps it terminate in case something goes wrong
# Feel free to change the function signatures
def K_means(X, m, maxT=10, threshold=.001, plots=False):
    (d, k), n = m.shape, X.shape[0]
    # WRITE CODE HERE:
    for t in range(maxT):
        # 1. Assignment step - Assign each sample to the nearest mean vector:

        # 2. Conditionally plot clustering results:
        if plots: scatter_by_category(X, labels, m)
        # 3. Mean updating - Update the cluster means
        # Note: be careful with empty clusters; you don't want division by 0
        # If there are empty clusters, simply don't update their centers
        # 4. Check for convergence
        if 0 <= threshold: break # CHANGE THIS
    return labels, obj, t+1
    
# Do NOT change the seed in any problem, unless you are explicitly allowed to
rng = np.random.default_rng(5)
num_clusters = 3
m = np.asfortranarray(rng.random((2, num_clusters)))
ne.evaluate('6*m - 3', out=m) # Magic formula ;)
labels, obj, iters = k_means(DATA, m, plots=True)
print('The algorithm took', iters, 'iterations to terminate.')
print('WCSS:', obj)
print(m)

# + [markdown] id="L9zUlW4Sm_ib"
# ##**Question 2:** The MNIST Dataset (25 points)
# Here comes your first decently-sized dataset this semester, the MNIST Dataset. The dataset contains 70000 28x28 preprocessed images of handwritten didits. Before proceeding further, you should familiarize yourself with the dataset, and we have provided a few cells to help you. You can modify them as long as you don't mess up the dataset.
#
# Due to the large size of MNIST, the following cells might take a while to run, so be patient.

# + colab={"base_uri": "https://localhost:8080/"} id="pJhSt35qmOLK" outputId="1f4aafd1-9057-492e-83e0-25e3923250f4"
# Fetching the MNIST dataset takes a while
# You only need to run this cell once per session
start = time.process_time()
from sklearn.datasets import fetch_openml
X, y = fetch_openml('mnist_784', return_X_y=True, as_frame=False)
print(time.process_time()-start, 'seconds')

# + colab={"base_uri": "https://localhost:8080/", "height": 362} id="Ky7NAdWTQHOU" outputId="1502a37d-1099-4588-c07c-cb1823840d06"
print(min(X[0]), max(X[0]))
plt.imshow(X[0].reshape(28, 28), 'gray')
plt.show()

print('Fraction of X that is non-zero:', np.count_nonzero(X)/X.size)
print(X.shape)
print(X.nbytes, 'bytes')
print('Are row elements contiguous?', X.data.c_contiguous)
print(X.dtype)

# + [markdown] id="hvV_x4EJQegC"
# ###**a.** Preprocessing for speed (5 points)
# Preprocessing is critical in ML for many reasons. One of them is that it can make the data more friendly to both your code and your hardware, potentially saving you tons of time. Based on the observations above, what things can you do to the MNIST data to speedup your algorithms without compromising the data? Explain why they are beneficial for speed for full credit.
#
# Notes: 
# 1. Although recommended, you do not have to apply all of your answers here in code.
# 2. Assume your code will be run on typical modern CPU.

# + colab={"base_uri": "https://localhost:8080/"} id="HLioIfSEHKdv" outputId="5534a5f8-ba04-44bb-8381-146de8dde86f"
print(max(X[0])) # A hint to one of the things we did to X

# + [markdown] id="Pof7JDXSG60w"
# ###**b.** K-means on MNIST (10 points)
# Now peform K-means on MNIST. We have provided the initializations for you, in the variable m2. This dataset will stress test your K-means implementation, so it's important that you use optimized routines. You can lower maxT if it takes too long, but no lower than 15. If your algorithm can't terminate in 60 seconds in our testing, there will be some penalty. If you really have trouble getting it to terminate in a reasonable amount of time, take a subset of X by performing something like Xs = X[:10000], but there will be some penalty.
#
# Hints:
# 1. Avoid loops that iterate over the dataset. If you really hav issue vectorizing, try to use the nopython JIT compiler (njit) from Numba to speed up your loops. See Question 3 code comments for more details on njit. This is not recommeneded however, as writing a jitted function usually requires more work than simply finding vectorized equivalents. You are also unlikely to match the speed of Numpy's vectorized routines this way.
# 2. $||\mathbf a - \mathbf b||_2^2 = \mathbf{a \cdot a} + \mathbf{b \cdot b} - 2\mathbf{a \cdot b}$
# 3. $\underset{\mathbf b \in \mathbf B}{\arg\min} ||\mathbf a - \mathbf b||_2^2 = \underset{\mathbf b \in \mathbf B}{\arg\min}~\mathbf{b \cdot b} - 2\mathbf{a \cdot b}$

# + colab={"base_uri": "https://localhost:8080/"} id="kJV4Im7qYit_" outputId="804dafdf-1b79-4971-95d7-ea7b537ac022"
# Definitely do NOT try naive printing/plotting this dataset
# Do NOT set plots=True when calling k_means
rng = np.random.default_rng(70)
# USE THESE parameters; note the names are different from those
# in question 1 to prevent accidental plotting
m2 = np.asfortranarray(rng.random((784, 10), np.float32))
m2/= 4
maxT = 30
threshold = 16
start = time.process_time()
# Run K-means on MNIST HERE:
print(time.process_time()-start, 'seconds')
print('The algorithm took', iters, 'iterations to terminate.')
print('WCSS:', obj)

# + [markdown] id="u1uVLyzzPHC2"
# ###**c.** Show the final centers as images (5 points)
# Hint: use subplots
#

# + colab={"base_uri": "https://localhost:8080/", "height": 313} id="eFhUNdHCVzaQ" outputId="5f1f8a60-2469-4f94-b8f9-71d08ef6685b"
# WRITE CODE HERE:
figsize=(12, 5)

# + [markdown] id="GlvkcnA7Pe6B"
# ###**d.** Avoiding bad initializations (5 points)
# A simple way to mitigate the risk of bad initializations is to use multiple random ones to obain multiple solutions, compare them, and keep the best one according to an objective function. Run K-means on MNIST 3 more times, each with a different initialization, and keep only the current best solution. Finally, show the best solution's centers as images and report WCSS.
#
# Note:
# 1. there is no time limit for this part, but you need to run your code to generate your own images. The run time shown in the preview is just for reference.
# 2. Use initializations methods analogous to the one in part b if you want to get the same results as shown.

# + colab={"base_uri": "https://localhost:8080/"} id="8RbkeoVzC1iM" outputId="24782fc3-e4f7-41c6-b5a3-0ee746ed016a"
start = time.process_time()
for i in range(3):
    m_tmp = np.asfortranarray(rng.random((784, 10), np.float32))
    m_tmp/= 4
    labels_tmp, obj_tmp, iters = k_means(X, m_tmp, 30, 16, Xs=Xs, X2=X2)
    if obj_tmp < obj: obj, m2 = obj_tmp, m_tmp

print(time.process_time()-start, 'seconds')
print('WCSS:', obj) # The lowest WCSS from the four solutions

# + colab={"base_uri": "https://localhost:8080/", "height": 313} id="Ik0cr5wAYkNo" outputId="1571cc7f-9221-4c82-d197-dc4d24f20835"
# WRITE CODE HERE:
figsize=(12, 5)


# + [markdown] id="e0N9X6iA43zx"
# ##**Question 3:** DP-means (40 points)
#
# Here, you will implement DP-means and test it out on the dataset from question 1.
#
#

# + [markdown] id="TrlNsiJB43zx"
# ###**a.** The role of $\lambda$ (5 points)
# Describe what the $\lambda$ parameter does in the DP-means algorithm discussed in class. 

# + [markdown] id="bkh8Rd3y43zx"
#

# + [markdown] id="uENwG6Zs43zz"
# ###**b.** Implement DP-means (25 points)
# Implement DP-means and apply it to the 2D Gaussian data from question 1. The code below contains comments of what needs to be coded up at a specific part of the code.  Please write your code below these comments where necessary. 
#
# Plot the created clusters per iteration in different colors. Similar to K-means, check how much the objecive value changed to determine convergence. The DP-means objective is given by $\lambda k+ \sum\limits_{i=1}^k Var(C_i)$
#
# Use $\lambda$ = 4 and plots=True so your function produce plots over the iterations.

# + id="pi7vvo1iVh5E"
def DP_means_obj(X, m, Lambda, labels, X2, m2):
    # WRITE CODE HERE:

## DP-means sub-routine:
# Remove @njit if you have issues with it and don't want to do the bonus 
# question. Numba can speedup for loops and Numpy sigificantly, but it often 
# has issues with python lists, scipy, pyplot, etc. It should be able to handle
# everything you need to do in onepass. Numba really prefers row vectors, 
# so passing in m.T instead of m is highly recommended. As always, you are
# welcome to change te function signature as long as you make it work.
@njit(fastmath=True, cache=True)
def onepass(X, mT, Lambda, labels, X2, m2): # welcome to change the signature
    (k, d), n = mT.shape, X.shape[0]
    # WRITE CODE HERE:
    # Iterating over X is required by the DP-means algorithm
    for i in range(n):
        # 1. Calculate distance from current point to all centers
        
        # 2. Check how far the clostest center is; add the currect point as
        # a center if necessary (add it to mT)
        if 0 > Lambda: # CHANGE THIS
            pass
        else:
            # assign point to the closest center
            pass
    # m should be returned because the updates performed may not be in place
    return mT.T # reversing the transpose


# + id="O15FHyzF7quq"
## DP-means main function:
def DP_means(X, m, Lambda=2, maxT=10, threshold=.001, plots=False):
    (d, k), n = m.shape, X.shape[0]
    # WRITE CODE HERE:
    for t in range(maxT):
        # 1. Implement the onepass function and use it here

        # 2. Conditionally plot clustering results:
        if plots: scatter_by_category(X, labels, m)
        # 3. Mean updating - Update the cluster means
        # Similar to K-means, be careful about empty clusters

        # 4. Check for convergence based on how much k and m changed
        if False: # CHANGE THIS
            break
    return m, t+1


# + colab={"base_uri": "https://localhost:8080/", "height": 707} id="3Gvjfg2Q7hdP" outputId="fd825ab6-bee8-4f09-c6cc-f26cdf2004d5"
rng = np.random.default_rng(0)
m = np.asfortranarray(rng.random((2, 1)))
# WRITE CODE HERE:
print(m)
print('Value of the objective:', obj)

# + [markdown] id="XpWnlOxH7unT"
# ###**c.** Reflection (4 points)
# If your implementation runs correctly, you should see the solution having one more center than the "expected" 3. **Explain why this happens and suggest a simple, generalizable fix for this.** (generlizable means you can apply this fix easily on other datasets). Also demonstrate this fix in code for full credit.

# + [markdown] id="rLph8NyQbdIa"
# This is because the initial center is provided at random and is not close enough to any point. Choosing a random point from the dataset to serve as the initial center would solve this issue.

# + colab={"base_uri": "https://localhost:8080/", "height": 723} id="ob8mNUcXodmM" outputId="99082d78-7e7e-496f-b8e8-b9a146522d94"
# WRITE CODE HERE:
print(m)
print('The algorithm took', iters, 'iterations to terminate.')
print('Value of the objective:', obj)

# + [markdown] id="Ppy9jyH843z0"
# ###**d.** Try different $\lambda$ values (6 points)
# Run DP-means with the following $\lambda$ values:
#
# $\underline{\lambda} = [1,9,16]$
#
# For each $\lambda$ value, plot the resulting clusters and report the final objective value.
#
# Note: The solution shown is obtained without the expected fix in part c.

# + colab={"base_uri": "https://localhost:8080/", "height": 1000} id="Buw3140OcmvU" outputId="da6870b3-fde6-4bb7-c14a-c152ca5b2f75"
rng = np.random.default_rng(50)
print('Lambda = 1')
# WRITE CODE HERE:
print('The algorithm took', iters, 'iterations to terminate.')
print('Value of the objective:', obj)

# + colab={"base_uri": "https://localhost:8080/", "height": 1000} id="qedjA4sWgm_b" outputId="26bc03ad-e74e-447f-af3e-56b556572fa1"
print('Lambda = 9')
# WRITE CODE HERE:
print(m)
print('The algorithm took', iters, 'iterations to terminate.')
print('Value of the objective:', obj)

# + colab={"base_uri": "https://localhost:8080/", "height": 740} id="JQlLRgztguZS" outputId="6174ce89-e110-4e74-b4ef-5467fe142b33"
print('Lambda = 16')
# WRITE CODE HERE:
print(m)
print('The algorithm took', iters, 'iterations to terminate.')
print('Value of the objective:', obj)

# + [markdown] id="fEchEfeeyQo-"
# ##**Question 4:** DP-means on MNIST (Bonus 10 points)
# This extra credit is fairly open-ended, just get **your** DP-means implementation to run reasonably well on MNIST.
#
# You may run DP-means on a subset of MNIST, but the subset should contain at least 10000 samples.
#
# Aim to achieve the following (ranked by importance):
# 1. Show the final centers as images. Your images need to look like those shown in the preview.
# 2. Find a suitable $\lambda$ and initial m that produce 10 to 30 clusters.
# 3. maxT must be at least 10. (Our solution uses 30)
# 4. The run time should be within 60 seconds.
#
# Notes:
# 1. $\lambda$ needs to be relatively large. It's better to start high and then decrease to find an appropriate value.
# 2. All solution previews in this homework are generated from running on the entire dataset.

# + colab={"base_uri": "https://localhost:8080/"} id="tVkK_EITCBnI" outputId="7b556c5a-d3dd-4682-acc3-a6bf700438ff"
# Initialize HERE:

start = time.process_time()
# Run DP-means on MNIST HERE:
print(time.process_time() - start, 'seconds')
print('Shape of the means matrix:', m2.shape)
print('The algorithm took', iters, 'iterations to terminate.')
print('Value of the objective:', obj)

# + colab={"base_uri": "https://localhost:8080/", "height": 485} id="im5l4MWmer8e" outputId="8d668af7-b54f-423d-8a97-b455247a28f2"
# WRITE CODE HERE:
figsize=(10, 8)

StarcoderdataJupyterScriptsDedupFiltered