{"cells":[{"metadata":{"_uuid":"48f90084ded00a2e6d2487ea88fd0a764538b3b0"},"cell_type":"markdown","source":"# Simple Clustering Analysis in the Frequency Domain\n\n**In this kernel, we are going to present some simple observations based on clusters created based on the frequency spectrum of the signals.  \nThis competition provides a good challenge (and a lot of fun!) regarding the stability of the training and the differences beetwen results in the local CV and LB score. While the discussions on adversarial validation can help with the issues, maybe the information in the clusters provide some complementary insights.**"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom joblib import Parallel, delayed\nimport pyarrow.parquet as pq\nfrom sklearn.metrics import silhouette_score\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.cluster import KMeans, DBSCAN\nfrom sklearn.decomposition import PCA","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"478b5b408b041554c657db86dfe2d97c1f0428c3"},"cell_type":"code","source":"N_THREADS = 4\nFEATURES_DIM = 1000\nCLUSTER_DIM = 100\nN_SAMPLES = 3\nRANDOM_SEED = 2019\n\nnp.random.seed(RANDOM_SEED)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8958fe55c0fbbe1aae56eb09dce35fa5463fcde5"},"cell_type":"markdown","source":"**Let's start by defining the features that will be used to cluster the signals.  \nWe are going to use a simple median of the slices of the frequency spectrum from a discrete Fourrier transform. It's a very simple feature space to start the analysis, but let's Keep It Sweet & Simple!  \nWe can expect a very skewed distribution of the data in the spectrum slices, but we are not going to transform the data this time.  \nLet's stop to write and start to coding. Here is the function to extract the features of a list of `signal_id`'s**\n\n    def cluster_features(signals, dataset='train'):\n        if dataset == 'train':\n            data = pq.read_pandas('../input/train.parquet', columns=[str(s) for s in signals]).to_pandas().values\n        else:\n            data = pq.read_pandas('../input/test.parquet', columns=[str(s) for s in signals]).to_pandas().values\n        features = np.zeros((data.shape[1], FEATURES_DIM))\n        for i, signal in enumerate(data.T):\n            fft = np.fft.rfft(signal)\n            fft = np.abs(fft)\n            fft = np.array_split(fft, FEATURES_DIM)\n            features[i] = [np.median(d) for d in fft]\n    return features\n    \n**In this case, the spectrum is been splitted in ranges of 20 KHz. Later, we will use PCA to reduce the dimension of the feature space.**"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"def load_signal(signal_id):\n    if signal_id <= 8711:\n        signal = pd.read_parquet('../input/train.parquet', columns=[str(signal_id)])\n    else:   \n        signal = pd.read_parquet('../input/test.parquet', columns=[str(signal_id)])\n    return np.squeeze(signal.values)\n\ndef plot_grid(id_measurements=None, signal_id=None, frame_size=(5,3), preprocessing=None):\n    meta_df = pd.read_csv('../input/metadata_train.csv')\n    meta_df2 = pd.read_csv('../input/metadata_test.csv')\n    meta_df = pd.concat((meta_df, meta_df2), axis=0, ignore_index=True, sort=False)\n    if signal_id is not None:\n        id_measurements = meta_df.loc[meta_df['signal_id'].isin(signal_id), 'id_measurement'].unique()\n    n_imgs = len(id_measurements)\n    fig = plt.figure(figsize=(frame_size[0]*3, frame_size[1]*n_imgs))\n    img_cursor = 1\n    for i in id_measurements:\n        signal_ids = meta_df.loc[meta_df['id_measurement']==i, 'signal_id'].values\n        targets = meta_df.loc[meta_df['id_measurement']==i, 'target'].values\n        for s,t in zip(signal_ids, targets):\n            color = '#007fff'\n            if signal_id is not None and s not in signal_id:\n                color = '#c0c0c0'\n            ax = fig.add_subplot(n_imgs, 3, img_cursor)\n            signal = load_signal(s)\n            if preprocessing is not None: \n                ax.plot(signal, color)\n                ax.plot(preprocessing(signal))\n            else:\n                ax.plot(signal, color) \n            ax.ticklabel_format(style='sci',scilimits=(-3,4),axis='x')\n            ax.grid(True)\n            ax.set_title('signal_id = {} , target = {}'.format(str(s), t))\n            img_cursor += 1\n    plt.tight_layout()\n    \ndef cluster_features(signals, dataset='train'):\n    if dataset == 'train':\n        data = pq.read_pandas('../input/train.parquet', columns=[str(s) for s in signals]).to_pandas().values\n    else:\n        data = pq.read_pandas('../input/test.parquet', columns=[str(s) for s in signals]).to_pandas().values\n    features = np.zeros((data.shape[1], FEATURES_DIM))\n    for i, signal in enumerate(data.T):\n        fft = np.fft.rfft(signal)\n        fft = np.abs(fft)\n        fft = np.array_split(fft, FEATURES_DIM)\n        features[i] = [np.median(d) for d in fft]\n    return features","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4d96efce2a80d08cc0c6f651f20d7a2357a3e41b"},"cell_type":"markdown","source":"**Let's load the features for the training and test datasets (s2 joblib!)**"},{"metadata":{"trusted":true,"_uuid":"a9fa139cc6dc03732733da246afff82f89cff6da"},"cell_type":"code","source":"train_meta = pd.read_csv('../input/metadata_train.csv')\ntrain_meta = train_meta.loc[train_meta['phase'] == 0]\ntrain_sig_ids = train_meta['signal_id'].values\ntrain_feat = Parallel(n_jobs=N_THREADS, verbose=2)(delayed(cluster_features)(s, 'train') for s in np.array_split(train_sig_ids, 3*N_THREADS))\ntrain_feat = np.concatenate(train_feat, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2588184a8d666982d38d5c10365ba2025a1f77ab"},"cell_type":"code","source":"test_meta = pd.read_csv('../input/metadata_test.csv')\ntest_meta = test_meta.loc[test_meta['phase'] == 0]\ntest_sig_ids = np.asarray([str(s) for s in test_meta['signal_id'].values])\ntest_feat = Parallel(n_jobs=N_THREADS, verbose=2)(delayed(cluster_features)(s, 'test') for s in np.array_split(test_sig_ids, 3*N_THREADS))\ntest_feat = np.concatenate(test_feat, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7331ede9099363ed7c88b51f3137bb856d49dbcd"},"cell_type":"markdown","source":"**The number of components in the PCA was enough to explain almost 99% of the variance of the training data and, as we were expecting, the features present a lot of outliers.**"},{"metadata":{"trusted":true,"_uuid":"00c0e70f7e30c9ad618ecd9530fde04d745de413","_kg_hide-input":true},"cell_type":"code","source":"features_idx = np.random.choice(np.arange(CLUSTER_DIM), 10)\nnorm_feat = MinMaxScaler().fit_transform(train_feat)\npca = PCA(n_components=CLUSTER_DIM).fit(norm_feat)\nnorm_feat = pca.transform(norm_feat)\nprint('PCA Explained Variance: {:.2f}%'.format(pca.explained_variance_ratio_.sum()*100))\nfig = plt.figure(figsize=(12,6))\nax = fig.add_subplot('111')\nax.set_title ('Random Subset of Features')\n_ = sns.boxplot(data=norm_feat[:, features_idx], ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"713ce59dcff5a300e7eb65ab05c2e2aeca581bea"},"cell_type":"markdown","source":"**Now, we are going to select the number of clusters based on the [Silhouette Score](https://en.wikipedia.org/wiki/Silhouette_%28clustering%29). For now, just the training data is been used.**"},{"metadata":{"trusted":true,"_uuid":"3601e429e6a68aa5deeb0bb5e7bc5083eefd49c0","_kg_hide-input":true},"cell_type":"code","source":"n_clusters = np.arange(2, 10)\n\ntrain_norm_feat = MinMaxScaler().fit_transform(train_feat)\npca = PCA(n_components=CLUSTER_DIM).fit(train_norm_feat)\ntrain_norm_feat = pca.transform(train_norm_feat)\nprint('PCA Explained Variance: {:.2f}%'.format(pca.explained_variance_ratio_.sum()*100))\ndef get_silhouette(n):\n    clt = KMeans(n_clusters=n, random_state=RANDOM_SEED).fit(train_norm_feat)\n    #clt = DBSCAN(0.5*n).fit(train_norm_feat)\n    try:\n        result = silhouette_score(train_norm_feat, clt.labels_)\n    except:\n        result = -1\n    return result\nsilhouette = Parallel(n_jobs=N_THREADS)(delayed(get_silhouette)(n) for n in n_clusters)\nplt.plot(np.arange(2, len(silhouette)+2), silhouette)\nplt.xlabel('N_Clusters')\nplt.ylabel('Silhoutte Score')\nplt.grid(True)\n_ = plt.title('Clustering Score - Training Set')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"395627b955f318eef7e1b89aa41c926c65bc172a"},"cell_type":"markdown","source":"**It's clear that five cluster work well in our training data. Now we can see how the `target` are spread across these clusters: **"},{"metadata":{"trusted":true,"_uuid":"661574b25ae70d29dd319c8257537d683ee9c5bc","_kg_hide-input":true},"cell_type":"code","source":"best_n_clusters = 5\nclt = KMeans(n_clusters=best_n_clusters, random_state=RANDOM_SEED).fit(train_norm_feat)\ntrain_clusters = pd.DataFrame(clt.labels_, columns=['cluster'], index=pd.Index([int(x) for x in train_sig_ids], name='signal_id'))\ntrain_clusters['target'] = train_meta.loc[train_meta['phase']==0, 'target'].values\nstats_df = train_clusters.groupby('cluster')['target'].agg(['count','sum'])\nstats_df.columns = ['count', 'pos_count']\nstats_df['pos_rate'] = stats_df['pos_count']/stats_df['count']\nstats_df['pos_rate'] = stats_df['pos_rate'].apply(lambda x: '{:.1f} %'.format(x*100))\ndisplay(stats_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a8291f33a57bddc8346ea372a707c0e7949974dd"},"cell_type":"markdown","source":"**Interesting to note the clusters 4 and 2. The former virtually has no positive samples while tha latter has almost 50% of positive samples.**  \n**Let's plot a small sample of the signals in each cluster to get a feel of how they look like.**"},{"metadata":{"_uuid":"50918d09ffadb76bc18c58e63917e011f620a173"},"cell_type":"markdown","source":"### **Cluster 0**"},{"metadata":{"trusted":true,"_uuid":"c14ad950e478a0a1e1a6b271c34af33c9b26293f","_kg_hide-input":true},"cell_type":"code","source":"plot_ids = train_clusters.loc[train_clusters['cluster']==0].sample(N_SAMPLES).index\nplot_grid(signal_id=plot_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af11af29719dc75528c8596051cc18a91692e82c"},"cell_type":"markdown","source":"### **Cluster 1**"},{"metadata":{"trusted":true,"_uuid":"61b738e8e1a15cb0613221b7de93a48f56a0d7c0","_kg_hide-input":true},"cell_type":"code","source":"plot_ids = train_clusters.loc[train_clusters['cluster']==1].sample(N_SAMPLES).index\nplot_grid(signal_id=plot_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d49b80881816b6079b61d2f033ab4b2f14279451"},"cell_type":"markdown","source":"### **Cluster 2**"},{"metadata":{"trusted":true,"_uuid":"f30dc7b4951f03c6b62e1d933870fa97d876f4e9","_kg_hide-input":true},"cell_type":"code","source":"plot_ids = train_clusters.loc[train_clusters['cluster']==2].sample(N_SAMPLES).index\nplot_grid(signal_id=plot_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c244c3a4e9051a60a2ffa8b88269103160737cbf"},"cell_type":"markdown","source":"### **Cluster 3**"},{"metadata":{"trusted":true,"_uuid":"f844ba1883176c54275342936ccebaec60d3c29a","_kg_hide-input":true},"cell_type":"code","source":"plot_ids = train_clusters.loc[train_clusters['cluster']==3].sample(N_SAMPLES).index\nplot_grid(signal_id=plot_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b5e619659d0db07d2a75518bb36cebbd31047ae4"},"cell_type":"markdown","source":"### **Cluster 4**"},{"metadata":{"trusted":true,"_uuid":"b8cd2ce256613ea29228f7779403caa5b1b9874c","_kg_hide-input":true},"cell_type":"code","source":"plot_ids = train_clusters.loc[train_clusters['cluster']==4].sample(N_SAMPLES).index\nplot_grid(signal_id=plot_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"942407a83475ceda230072af1d13197426a1a6b7"},"cell_type":"markdown","source":"**Now, the clusters will be created using just the test set.**"},{"metadata":{"trusted":true,"_uuid":"acedc76aa86c84444b47e154574f1e6597af9bf8","_kg_hide-input":true},"cell_type":"code","source":"test_norm_feat = MinMaxScaler().fit_transform(test_feat)\npca = PCA(n_components=CLUSTER_DIM).fit(test_norm_feat)\ntest_norm_feat = pca.transform(test_norm_feat)\nprint('PCA Explained Variance: {:.2f}%'.format(pca.explained_variance_ratio_.sum()*100))\ndef get_silhouette(n):\n    clt = KMeans(n_clusters=n, random_state=RANDOM_SEED).fit(test_norm_feat)\n    return silhouette_score(test_norm_feat, clt.labels_)\nsilhouette = Parallel(n_jobs=N_THREADS)(delayed(get_silhouette)(n) for n in n_clusters)\nplt.plot(np.arange(2, len(silhouette)+2), silhouette)\nplt.xlabel('N_Clusters')\nplt.ylabel('Silhoutte Score')\nplt.grid(True)\n_ = plt.title('Clustering Score - Test Set')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"61955891b6147aa2cd78d17bb00e84ac5ea1200e"},"cell_type":"markdown","source":"**Good! We've got five cluster again.**  \n**The number of samples in each cluster is the following: **"},{"metadata":{"trusted":true,"_uuid":"766878371b2ce476b883269b6566503efb122b56","_kg_hide-input":true},"cell_type":"code","source":"best_n_clusters = 5\nnorm_feat = MinMaxScaler().fit_transform(test_feat)\nclt = KMeans(n_clusters=best_n_clusters, random_state=RANDOM_SEED).fit(norm_feat)\ntest_clusters = pd.DataFrame(clt.labels_, columns=['cluster'], index=test_meta['signal_id'])\nstats_df = test_clusters.reset_index().groupby('cluster').count()\nstats_df.columns = ['count']\ndisplay(stats_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9366389eb999c12f26742ae5bf924586db648c68"},"cell_type":"markdown","source":"**Let's straight to the clusters using the whole dataset (training + test). We hope to find the same clusters in both dataset.**"},{"metadata":{"trusted":true,"_uuid":"9c6d16784c3658f540c40d04225dcd248c00191d","_kg_hide-input":true},"cell_type":"code","source":"norm_feat = MinMaxScaler().fit_transform(np.concatenate((train_feat, test_feat), axis=0))\npca = PCA(n_components=CLUSTER_DIM).fit(norm_feat)\nnorm_feat = pca.transform(norm_feat)\nprint('PCA Explained Variance: {:.2f} %'.format(pca.explained_variance_ratio_.sum()*100))\n\ndef get_silhouette(n):\n    clt = KMeans(n_clusters=n, random_state=RANDOM_SEED).fit(norm_feat)\n    return silhouette_score(norm_feat, clt.labels_)\nsilhouette = Parallel(n_jobs=N_THREADS)(delayed(get_silhouette)(n) for n in n_clusters)\nplt.plot(np.arange(2, len(silhouette)+2), silhouette)\nplt.xlabel('N_Clusters')\nplt.ylabel('Silhoutte Score')\nplt.grid(True)\n_ = plt.title('Clustering Score - Test Set')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"53391ab72ab9dfb972e6324ed9df238a1fc691fe"},"cell_type":"markdown","source":"**Now, it seems the data have 6 clusters! We have some cluster in a dataset that don't matches a cluster in the other one**\n**Getting some statistics in each cluster: **"},{"metadata":{"trusted":true,"_uuid":"78dff001e09d56f6bdb7491976676b62c3413e26","_kg_hide-input":true},"cell_type":"code","source":"best_n_clusters = 6\nnorm_feat = MinMaxScaler().fit_transform(np.concatenate((train_feat, test_feat), axis=0))\npca = PCA(n_components=CLUSTER_DIM).fit(norm_feat)\nnorm_feat = pca.transform(norm_feat)\nclt = KMeans(n_clusters=best_n_clusters, random_state=RANDOM_SEED).fit(norm_feat)\nclusters = pd.DataFrame(clt.labels_, columns=['cluster'], index=pd.Index([int(x) for x in np.concatenate((train_sig_ids, test_sig_ids), 0)], name='signal_id'))\nclusters['dataset'] = 0\nclusters.iloc[train_sig_ids.shape[0]:, 1] = 1\nclusters['target'] = 0\nclusters.iloc[:train_sig_ids.shape[0], -1] = train_meta['target'].values\nstats_df = pd.DataFrame()\nstats_df['count'] = clusters['cluster'].value_counts()\nstats_df['train'] = clusters.groupby('cluster')['dataset'].apply(lambda x: (x==0).sum())\nstats_df['test'] = clusters.groupby('cluster')['dataset'].apply(lambda x: (x==1).sum())\nstats_df['pos_count'] = clusters.groupby('cluster')['target'].sum()\nstats_df['pos_rate'] = clusters.groupby('cluster')['target'].sum()/stats_df['train']\nstats_df['pos_rate'] = stats_df['pos_rate'].apply(lambda x: '{:.2f} %'.format(x*100))\ndisplay(stats_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3bdf19624ba4cc6e95b0a169d3fe15df3b242d96"},"cell_type":"markdown","source":"**We have a couple of interesting obsevations in the above dataframe:**\n1. The cluster 5 matches the cluster 2 in the training analysis and it has no proper equivalent in the test set\n2. The cluster 2 matches the cluster 1 in the test analysis and it has no proper equivalent in the training set\n3. The proportion of training/test samples in cluster 0 matches the proportion of training/test samples in the dataset => Can we expect the same positive rate in the test set - cluster 0 ?  \n\n**Let's take a look in the clusters: **"},{"metadata":{"_uuid":"a8534b6adcd715c05bb46218bd2c79e976e88399"},"cell_type":"markdown","source":"### **Cluster 0**"},{"metadata":{"trusted":true,"_uuid":"37b93602bdb6d89a3d3f128a787c2917dc717870","_kg_hide-input":true},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==0].sample(N_SAMPLES).index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8e9fdf292160ae33f184691e0e9d84523d710d0f"},"cell_type":"markdown","source":"### **Cluster 1**"},{"metadata":{"trusted":true,"_uuid":"931df7a746ba04e0cca114979d3a1fc22cc34602","_kg_hide-input":true},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==1].sample(N_SAMPLES).index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"90754083ad5971c885bc1376d898032bf8a1aa3a"},"cell_type":"markdown","source":"### **Cluster 2 (Problematic)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"43df8d49b454bb3e0bb8c3d38c601989b19b6911"},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==2].index[:6])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"aa4c25e44e2b23ee7e1b29e59d3dea98f3149aa5"},"cell_type":"markdown","source":"### **Cluster 3**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"a922f3674cfe6abcddd78fe79c925f244f5b98f9"},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==3].sample(N_SAMPLES).index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a614e01e45c83cd4a814814392ae5bc257356eaa"},"cell_type":"markdown","source":"### **Cluster 4**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"2367a357bd654feb7a48730974021a5b8efa481f"},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==4].sample(N_SAMPLES).index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2392ca5ba587d120ade37ce60fc0ad0fd43001ea"},"cell_type":"markdown","source":"### **Cluster 5 (Problematic)**"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"2546e68179c6978bb473e4fc05d4850b372ac675"},"cell_type":"code","source":"plot_grid(signal_id=clusters.loc[clusters['cluster']==5].index[-6:])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"665ca89fa3707708fff0f25169eaa10a0453db75"},"cell_type":"code","source":"train_clusters.to_csv('train_clusters.csv')\ntest_clusters.to_csv('test_clusters.csv')\nclusters.to_csv('clusters.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1c7acbf910a039f41e5694915d9a4df2a7202c2f"},"cell_type":"markdown","source":"# Final Notes\n\n**This analysis still have a wide margin to improvements. Note the clusters are imperfect, e.g the test examples in the cluster 5 seems more similar to the cluster 4 and the clusters 2 and 3 are similar except the number of peaks.  \nYet, it can maybe help in a post processing phase. As an example, all positive predictions in my current solution are placed in the cluster 0. It can also be used to stratify the folds in the training **\n\n### **That's all Folks!**\n### **Please, leave your suggestions, comments or feedback**"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}