{"cells":[{"metadata":{"_uuid":"522e1d3bb63641df96471b3620cf576220b84017"},"cell_type":"markdown","source":"# Unsupervised Learning - Clustering\nVisualization of the data clustering and the dimension reduction, from the dataset containing pixels whose intensity vary over the 48 000 images (saved in train_changing_pixels_df.csv)"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"54722f277b2160de863a8443a513913a3017a729"},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.cluster.hierarchy import linkage, dendrogram\nfrom scipy.sparse import csr_matrix","execution_count":1,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"c8097b36b9878f3ad448aa22ce82362b48eaa6af"},"cell_type":"code","source":"from sklearn.manifold import TSNE\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.cluster import KMeans\nfrom sklearn.svm import SVC\nfrom sklearn.preprocessing import StandardScaler, Normalizer, MaxAbsScaler\nfrom sklearn.decomposition import PCA, NMF, TruncatedSVD\nfrom sklearn.model_selection import train_test_split, RandomizedSearchCV, GridSearchCV","execution_count":2,"outputs":[]},{"metadata":{"collapsed":true,"scrolled":true,"trusted":false,"_uuid":"4774bce345dbf640155746f100734d29ccce7077"},"cell_type":"code","source":"#Loading datasets\ntrain_changing_pixels_df = pd.read_csv('../input/changing-pixels-csv/train_changing_pixels_DB.csv', index_col=0)\ntrain = pd.read_csv('../input/digit-recognizer/train.csv')\nholdout = pd.read_csv('../input/digit-recognizer/test.csv')\n#print(train_changing_pixels_df.head())\n#print(train_simplified_pixels_df.head())\n#print(train.head())\n\n#Pixels with constant intensities that have been dropped in the train_changing_pixels_df\nDROPPED_PIX = ['pixel0', 'pixel1', 'pixel2', 'pixel3', 'pixel4', 'pixel5', 'pixel6', 'pixel7', 'pixel8', 'pixel9', \n               'pixel10', 'pixel11', 'pixel16', 'pixel17', 'pixel18', 'pixel19', 'pixel20', 'pixel21', 'pixel22', \n               'pixel23', 'pixel24', 'pixel25', 'pixel26', 'pixel27', 'pixel28', 'pixel29', 'pixel30', 'pixel31', \n               'pixel52', 'pixel53', 'pixel54', 'pixel55', 'pixel56', 'pixel57', 'pixel82', 'pixel83', 'pixel84', \n               'pixel85', 'pixel111', 'pixel112', 'pixel139', 'pixel140', 'pixel141', 'pixel168', 'pixel196', \n               'pixel392', 'pixel420', 'pixel421', 'pixel448', 'pixel476', 'pixel532', 'pixel560', 'pixel644', \n               'pixel645', 'pixel671', 'pixel672', 'pixel673', 'pixel699', 'pixel700', 'pixel701', 'pixel727', \n               'pixel728', 'pixel729', 'pixel730', 'pixel731', 'pixel754', 'pixel755', 'pixel756', 'pixel757', \n               'pixel758', 'pixel759', 'pixel760', 'pixel780', 'pixel781', 'pixel782', 'pixel783']","execution_count":3,"outputs":[]},{"metadata":{"_uuid":"fa46395d228a8c4c8c07e64a68582b7479e69364"},"cell_type":"markdown","source":"## Clustering visualizations"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"a7e9dcea7938ca4fa28f0fe1a4d04c7241ef6c48"},"cell_type":"code","source":"#Preparing samples and labels arrays\ns = np.random.randint(low=0, high=42001, size=8200).tolist()\nsamples = train_changing_pixels_df.iloc[s, :].values\ndigits = train.iloc[s, :]['label'].tolist()\nprint(samples.shape)\n\n#Dendrogram of the clustering of data\nmergings = linkage(samples, method='complete')\n_ = plt.figure(figsize=(30,20))\n_ = dendrogram(mergings, labels=digits)\n_ = plt.title('Dendrogram of the sampled changing pixels (not scaled/normalized)')\n_ = plt.legend(digits)\n_ = plt.savefig('Dendrogram.png')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"57b6d85d3159efc94c1288c7b5a90d09a54ac165"},"cell_type":"code","source":"#Preparing samples and labels arrays\ns = np.random.randint(low=0, high=42001, size=8400).tolist()\nsamples = train_changing_pixels_df.iloc[s, :].values\ndigits = train.iloc[s, :]['label'].tolist()\nprint(samples.shape)\n\n#t-SNE visualization of the data\nmodel = TSNE(learning_rate=100)  #between 50-200\ntransformed = model.fit_transform(samples)\nxs = transformed[:, 0]\nys = transformed[:, 1]\n_ = plt.figure(figsize=(30,20))\n_ = plt.scatter(xs, ys, c=digits, alpha=0.6, cmap = \"nipy_spectral\")\n_ = plt.title('t-SNE of the sampled changing pixels (not scaled/normalized)')\n_ = plt.legend(digits)\n_ = plt.savefig('t-SNE.png')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"93fabe7a68db85bef6cda87cfcf0f0ad87dc09d0"},"cell_type":"markdown","source":"*For a sample of 8200 images (20% of the train set), the t-SNE displays 9 clusters, for the 10 digits of the dataset, that are clear. However, there are 2x2 digits that look harder to dissociate from each other.*\n\n*For a sample of half the data (21 000 images) the cluster are clearer with still 2 of them that have a common area.*"},{"metadata":{"_uuid":"ea28448023e210423e2f50b13a077cc6d15fa380"},"cell_type":"markdown","source":"## Dimension Reduction\n### Principal Components Analysis (PCA)"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"d43a3fb0194f86b90be9d3ed0afc16b98990ce8a"},"cell_type":"code","source":"#Preparing samples and labels arrays\n#s = np.random.randint(low=0, high=42001, size=1050).tolist()\nsamples = train_changing_pixels_df.values  #.iloc[s, :]\ndigits = train['label'].tolist()  #.iloc[s, :]\nprint(samples.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"9c89035d919a35a07cc7201dd6f8704f49c9effc"},"cell_type":"code","source":"#PCA model\npca = PCA()\npca.fit(samples)\n\n#PCA features variance visualization\npca_features = range(pca.n_components_)\n_ = plt.figure(figsize=(30,20))\n_ = plt.bar(pca_features, pca.explained_variance_)\n_ = plt.xticks(pca_features)\n_ = plt.title('Principal Components Analysis for Dimension Reduction')\n_ = plt.xlabel('PCA features')\n_ = plt.ylabel('Variance of the PCA feature')\n_ = plt.savefig('PCA features variance.png')\nplt.show()\n\n#PCA features variance visualization - ZOOM in\nl= 100\nx = range(l)\n_ = plt.figure(figsize=(30,20))\n_ = plt.bar(x, pca.explained_variance_[:l])\n_ = plt.xticks(x)\n_ = plt.title('Principal Components Analysis for Dimension Reduction - Zoom In {} first features'.format(l))\n_ = plt.xlabel('PCA features')\n_ = plt.ylabel('Variance of the PCA feature')\n_ = plt.savefig('PCA features variance_zoom.png')\nplt.show()","execution_count":1,"outputs":[]},{"metadata":{"_uuid":"d58947558f021b073a321e8a8bb47810e8a7218d"},"cell_type":"markdown","source":"*The best number of PCA features to keep (reduced dimension) is <80. The 1st PCA feature represents a lot of the variance, as the 5 following ones.*"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"a57c169a049ed525dd28cbcf9c297150ad71580e"},"cell_type":"code","source":"#Visualization of the variance of the data carried by the number of PCA features\nn_components = np.array([1,2,3,4,5,6, 10, 30, 60, 80, 100, 200, 400, 700])\ncumul_variance = np.empty(len(n_components))\nfor i, n in enumerate(n_components):\n    pca = PCA(n_components=n)\n    pca.fit(samples)\n    cumul_variance[i] = np.sum(pca.explained_variance_ratio_)\n\nprint(cumul_variance)\n\n_ = plt.figure(figsize=(30,20))\n_ = plt.grid(which='both')\n_ = plt.plot(n_components, cumul_variance, color='red')\n_ = plt.xscale('log')\n_ = plt.xlabel('Number of PCA features', size=20)\n_ = plt.ylabel('Cumulated variance of data (%)', size=20)\n_ = plt.title('Data variance cumulated vs number of PCA features', size=20)\nplt.savefig('Cumulated variance_pca features.png')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c94ac50ef270025373ca89a21e49aed21609d505"},"cell_type":"markdown","source":"*50% of the data variance is carried by 10 PCA features, and 90% of the data variance is carried by 80 PCA features (vs 708 features initially). The reduction is huge !*\n\n*However, the PCE features have no interpretable meaning.*\n\n*--> The use of PCA() in the pipleine along with StandardScaler() and SVC() increased the accuracy up to 0.97 trained on half of the samples.*"},{"metadata":{"_uuid":"a054f54494d9c79d5681c5b687c358e205c9df64"},"cell_type":"markdown","source":"### Non-negative Matrix Factorization (NMF)"},{"metadata":{"trusted":false,"_uuid":"66d3b35e695844728a4d155bbd16569112784f5c"},"cell_type":"code","source":"#Preparing samples and labels arrays\ns = np.random.randint(low=0, high=42001, size=8200).tolist()\nsamples = train.drop(columns='label').values  #.iloc[s, :]\ndigitsa = train['label'].tolist()  #.iloc[s, :]\nprint(samples.shape)","execution_count":57,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"a66f1f715aba20d20c41ac8a4e3dd45ac1249511"},"cell_type":"code","source":"#Creating the NMF model and the features & components\nnmf = NMF(n_components=16)\nnmf_features = nmf.fit_transform(samples)\nnmf_components = nmf.components_\nprint(\"Shape of NMF features: {}, shape of NMF components: {}\".format(nmf_features.shape, nmf_components.shape))\n\n#Visualization of the features\nfor i, component in enumerate(nmf_components):\n    N = nmf_components.shape[0]\n    ax = plt.subplot((N//3)+1, 3, i+1)\n    bitmap = component.reshape((28,28))\n    plt.imshow(bitmap, cmap='gray')\n    plt.xticks([])\n    plt.yticks([])\nplt.title('NMF Components from the original images')\nplt.show()","execution_count":58,"outputs":[]},{"metadata":{"_uuid":"161c19d200cebf9a9046b1a87e2d2aa6e02195fe"},"cell_type":"markdown","source":"*The samples have been split between NMF features and components, in a way that Features x Components = Samples (with matrices product)*\n\n*The images have been decomposed into combinations of the NMF components displaed above.*"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"1ec91c20c3227ce6082bfb9e51972afa3d24e96e"},"cell_type":"code","source":"%%time\n\n#Sample randomly the dataset using discrete_uniform pick-up to reduce the amount of data\n#sample = np.random.randint(low=0, high=42001, size=8400).tolist()\nX = train_changing_pixels_df.values  #.iloc[sample, :]\nX = X / 255.0\ny = train['label'].values  #.iloc[sample, :]\nprint(\"Shape of X and Y arrays: {}\".format((X.shape, y.shape)))\n\n#Split the training set\nX_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2, random_state=4, stratify=y)\n\n#Yielding the scores according to number of NMF components\ncomponents = np.arange(1, 100)\nscores = np.empty(len(components))\nfor n in components:\n    pipeline = make_pipeline(NMF(n_components=n), SVC(kernel='rbf', cache_size=1000))\n    pipeline.fit(X_train, y_train)\n    scores[n-1] = pipeline.score(X_test, y_test)\n\n#Plotting of the scores evlution\n_ = plt.figure(figsize=(30,20))\n_ = plt.grid(which='both')\n_ = plt.plot(components, scores)\n_ = plt.xlabel('Number of NMF components', size=20)\n_ = plt.ylabel('Score obtained', size=20)\n_ = plt.title('Evolution of SVC classification score (samples={})'.format(len(y)), \n              size=30)\nplt.savefig('Score vs components NMF.png')\nplt.show()\n\nprint(\"Best score {} obtained for {} components\".format(scores.max(), scores.argmax()+1))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f7a7e09b254c656b55586f4ccf7ad6bee4da8289"},"cell_type":"markdown","source":"*Highest score of 0.83 obtained for 16 components*"},{"metadata":{"_uuid":"f861a09ba35bb587399c2fd505f6adac6b9ceae5"},"cell_type":"markdown","source":"### Sparce matrix\nAs the pixel inteisties present a lot of 0 values, we can turn the train 2D array into a sparse matrix to reduce the memory size\n\nhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html\n- TruncatedSVD for the equivalent of PCA()\n- NMF() directly with the csr_matrix"},{"metadata":{"scrolled":false,"trusted":false,"_uuid":"c0b40dd1d24d318e55dffad5ed22e728bb822fdd"},"cell_type":"code","source":"#Sparsity visuaization & sparse matrix creation\nsamples = train_changing_pixels_df.values\n_ = plt.figure(figsize=(10,100))\n_ = plt.spy(samples)\nplt.show()\n\nsparse_samples = csr_matrix(samples)\n\n#Memory Size comparison\ndense_size = samples.nbytes/1e6\nsparse_size = (sparse_samples.data.nbytes + \n               sparse_samples.indptr.nbytes + sparse_samples.indices.nbytes)/1e6\nprint(\"From {} to {} Mo in memory usage with the sparse matrix\".format(dense_size, sparse_size))\n\n#Dimension reduction using PCA equivalent for sparse matrix\nmodel = TruncatedSVD(n_components=10)\nmodel.fit(sparse_samples)\nreduced_sparse_samples = model.transform(sparse_samples)\nprint(reduced_sparse_samples.shape)","execution_count":4,"outputs":[]},{"metadata":{"_uuid":"f1cee7b339a5255cc0ee3bcc1750d1d3b7b4cd8a"},"cell_type":"markdown","source":"*The sparsity visualization confirms the relevance of the use of a csr matrix for the samples*"},{"metadata":{"_uuid":"4c3b554b20428bb4a0a74445c255805859e574fe"},"cell_type":"markdown","source":"---"},{"metadata":{"_uuid":"56b3d878debe49bb85770135c770bb858c068676"},"cell_type":"markdown","source":"## Clustering models\n### with PCA dimension reduction"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"816e5fa4f0a990bf5dc041ec82394240a3b716d8"},"cell_type":"code","source":"#Preparing samples and labels arrays\n#s = np.random.randint(low=0, high=42001, size=42000).tolist()\nsamples = train_changing_pixels_df.values  #.iloc[s, :]\ndigits = train['label'].tolist()  #.iloc[s, :]\nprint(\"Shape of samples: {}\".format(samples.shape))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"scrolled":true,"trusted":false,"_uuid":"1b01c745551902b0ba76f71d5474241e27f7515d"},"cell_type":"code","source":"%%time\n\n#K-Means clustering\npipeline = make_pipeline((StandardScaler()),\n                         PCA(n_components=80),    \n                         KMeans(n_clusters=10))  #10 digits to identify\npipeline.fit(samples)\nlabels = pipeline.predict(samples)\n\n#Predictions on holdout set\nX_holdout = holdout.drop(columns=DROPPED_PIX).values\npredictions = pipeline.predict(X_holdout)  #/!!\\ predicted labels don't correspond to digits !\n\n#Cross-tabulation\ndf = pd.DataFrame({'clusters_labels': labels, 'digits': digits })\ncross_tab= pd.crosstab(df['clusters_labels'], df['digits'])\nprint(cross_tab)\n\n#cluster-digits matching dict\ncluster_digit_dict = {}\nfor i, row in cross_tab.iterrows():\n    cluster_digit_dict[i] = row.idxmax()\nprint(cluster_digit_dict)\n\n#Submission form\nsubmission_df = pd.DataFrame({'ImageId': range(1,28001), 'cluster_label': predictions})\nsubmission_df['Label'] = submission_df['cluster_label'].map(cluster_digit_dict)\nsubmission_df.drop(columns=['cluster_label'], inplace=True)\n#submission_df.to_csv('submission_pca_kmeans_DB.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7ce57b51931363c43896ffc685671dfc79d06772"},"cell_type":"markdown","source":"*Digits 1, 4, 5 and 7 are harder to classify as they are represented in several clusters closely*"},{"metadata":{"_uuid":"378bae37e64b02e1e52b67d13a14e1de6e2d5a14"},"cell_type":"markdown","source":"### with NMF dimension reduction\n/!\\ NMF doesn't accept negative values -> StandardScaler() not possible here "},{"metadata":{"trusted":false,"_uuid":"c039a02d58bbb6e1bd8d738f0889025eb5c5a0cd"},"cell_type":"code","source":"#Preparing samples and labels arrays\n#s = np.random.randint(low=0, high=42001, size=42000).tolist()\nsamples = train_changing_pixels_df.values  #.iloc[s, :]\ndigits = train['label'].tolist()  #.iloc[s, :]\nprint(\"Shape of samples: {}\".format(samples.shape))","execution_count":15,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"bdae76b7786eb88a9ee6551e3ba213ed6bc1722b"},"cell_type":"code","source":"%%time\n\n#K-Means clustering\npipeline = make_pipeline(MaxAbsScaler(), \n                         NMF(n_components=18),\n                         KMeans(n_clusters=10))  #10 digits to identify\n\npipeline.fit(samples)\nlabels = pipeline.predict(samples)\n\n#Predictions on holdout set\nX_holdout = holdout.drop(columns=DROPPED_PIX).values\npredictions = pipeline.predict(X_holdout)  #/!!\\ predicted labels correspond to clusters !\n\n#Cross-tabulation\ndf = pd.DataFrame({'clusters_labels': labels, 'digits': digits })\ncross_tab= pd.crosstab(df['clusters_labels'], df['digits'])\nprint(cross_tab)\n\n#cluster-digits matching dict\ncluster_digit_dict = {}\nfor i, row in cross_tab.iterrows():\n    cluster_digit_dict[i] = row.idxmax()\nprint(cluster_digit_dict)\n\n#Submission form\nsubmission_df = pd.DataFrame({'ImageId': range(1,28001), 'cluster_label': predictions})\nsubmission_df['Label'] = submission_df['cluster_label'].map(cluster_digit_dict)\nsubmission_df.drop(columns=['cluster_label'], inplace=True)\n#submission_df.to_csv('submission_pca_kmeans_DB.csv', index=False)","execution_count":17,"outputs":[]},{"metadata":{"collapsed":true,"_uuid":"42a009db8464202efebc4718b8e075b26a403a99"},"cell_type":"markdown","source":"## Cosine Similarity\nhttps://towardsdatascience.com/building-improving-a-k-nearest-neighbors-algorithm-in-python-3b6b5320d2f8"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"f48d31842393fcb87399194589598affb72aaa7b"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.4"}},"nbformat":4,"nbformat_minor":1}