{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### How to ensemble clustering algorithms?\n\n> - Since this competition is about clustering \n> - And since it is necessary to ensemble in order to achieve maximum Kaggling. \n\nOn this notebook we show how to ensemble multiple clustering algorithms.\n\n![](https://i.ibb.co/Sy5xgnQ/Scheme-of-ensemble-clustering-approach.png)\n\n[[source](https://www.researchgate.net/figure/Scheme-of-ensemble-clustering-approach_fig1_277589770)]\n\n\nThe approach is simple: We create a sparse matrix of all pairs of samples and assign 1 to any pair of samples that are on the same cluster.\n\nWe then simply take the mean/median of all sparse matrices, apply a threshold, and convert it back to cluster ids.\n\nThis way we can use multiple clustering algorithms at the same time. \n\n\n_____\n\n##### Updated version\n\nMultiple upgrades and fixes had been proposed by [ehekatlact](https://www.kaggle.com/ehekatlact/) on the [this](https://www.kaggle.com/code/ehekatlact/tps2207-ultra-fast-ensemble/data) excellent notebook. \nThis notebook had been upgraded with the proposed fixes, it is now extremely fast to run and it also submit to the LB. \n\nAt the bottom of this notebook you can find a single function that performs the full computation end-to-end. \nCopy it wherever you want to perform easy clustering ensemble! \n\nHave fun! ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-04T15:41:53.439899Z","iopub.execute_input":"2022-07-04T15:41:53.440308Z","iopub.status.idle":"2022-07-04T15:41:53.45155Z","shell.execute_reply.started":"2022-07-04T15:41:53.440276Z","shell.execute_reply":"2022-07-04T15:41:53.450577Z"}}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom tqdm import trange\nimport matplotlib.pyplot as plt\nfrom sklearn.cluster import KMeans\nfrom collections import defaultdict\nfrom scipy.sparse import csr_matrix\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.metrics import adjusted_rand_score\n\ndf = pd.read_csv('../input/tabular-playground-series-jul-2022/data.csv')\ndf = df.drop(columns = 'id')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:36.514710Z","iopub.execute_input":"2022-07-08T07:54:36.515721Z","iopub.status.idle":"2022-07-08T07:54:50.844752Z","shell.execute_reply.started":"2022-07-08T07:54:36.515684Z","shell.execute_reply":"2022-07-08T07:54:50.843566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Now let's say we got two algorithms we want to ensemble\n\n- KMeans(n_clusters = 7)\n- Gaussian Mixture(n_components = 7)","metadata":{}},{"cell_type":"code","source":"clusters_k_means = KMeans(n_clusters = 7).fit_predict(df)\nclusters_mixture = GaussianMixture(n_components = 7).fit_predict(df)\n\nprint(clusters_k_means.shape)\nprint(clusters_mixture.shape)\n\nclusters_list = [clusters_k_means, clusters_mixture]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:36.514710Z","iopub.execute_input":"2022-07-08T07:54:36.515721Z","iopub.status.idle":"2022-07-08T07:54:50.844752Z","shell.execute_reply.started":"2022-07-08T07:54:36.515684Z","shell.execute_reply":"2022-07-08T07:54:50.843566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cls_tup_list = []\nfor cls_tup in zip(*clusters_list):\n    cls_tup_list.append(cls_tup)\nzipper = {x: i for i, x in enumerate(sorted(set(cls_tup_list)))}\nzipped_list = [zipper[x] for x in cls_tup_list]\nunzipper = defaultdict(set)\nfor idx, cls_tup in enumerate(cls_tup_list):\n    zipped = zipper[cls_tup]\n    unzipper[zipped].add(idx)\ncomp_clusters_list = [[-1]*len(zipper) for _ in range(len(clusters_list))]\nfor clusters, comp_clusters in zip(clusters_list, comp_clusters_list):\n    for i, cluster_i in enumerate(clusters):\n            value = zipped_list[i]\n            comp_clusters[value] = cluster_i","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:54:50.847090Z","iopub.execute_input":"2022-07-08T07:54:50.847782Z","iopub.status.idle":"2022-07-08T07:54:51.085411Z","shell.execute_reply.started":"2022-07-08T07:54:50.847736Z","shell.execute_reply":"2022-07-08T07:54:51.084213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### We create a sparse matrix of all pairs of samples and assign 1 to any pair of samples that are on the same cluster.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"def create_sparse_matrix(clusters):\n    n = len(clusters)\n    data = []\n    row = []\n    col = []\n    # O(n**2)\n    for i in trange(n):\n        for j in range(i+1, n):\n            if clusters[i] == clusters[j]:\n                data.append(1)\n                row.append(i)\n                col.append(j)\n    return csr_matrix((data, (row, col)), shape=(n, n))\n\nsparse_matrix_list = [create_sparse_matrix(comp_clusters) for comp_clusters in comp_clusters_list]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:51.088235Z","iopub.execute_input":"2022-07-08T07:54:51.088989Z","iopub.status.idle":"2022-07-08T07:54:51.109350Z","shell.execute_reply.started":"2022-07-08T07:54:51.088932Z","shell.execute_reply":"2022-07-08T07:54:51.108052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### We then simply take the mean/median of all sparse matrices, apply a threshold, and convert it back to cluster ids.\n\n![](https://i.ibb.co/84nnT31/1-s2-0-S153204641300021-X-gr1.jpg)\n[[source](https://www.sciencedirect.com/science/article/pii/S153204641300021X)]\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"sparse_matrix_mean = (sparse_matrix_list[0] * 0.5 + sparse_matrix_list[1] * 0.5) / 2","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:51.112484Z","iopub.execute_input":"2022-07-08T07:54:51.113443Z","iopub.status.idle":"2022-07-08T07:54:51.120311Z","shell.execute_reply.started":"2022-07-08T07:54:51.113399Z","shell.execute_reply":"2022-07-08T07:54:51.119050Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Now we apply the threshold\n\n> Special thanks to [ehekatlact](https://www.kaggle.com/ehekatlact) for coming up with this [fix](https://www.kaggle.com/code/ehekatlact/tps2207-ultra-fast-ensemble)! \n\nUse DSU to control the upper limit of each cluster size to prevent clusters from sticking together.\nIf you are interested in this phenomenon, try setting SIZ_MAX to float('inf') and experiment with it!","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# The lower the threshold the lower number of clusters: Since we are counting by the connected component.\n# It seems desirable to set it a little higher.\nthreshold = 0.5\nsparse_matrix_mean[sparse_matrix_mean < threshold] = 0","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:51.122070Z","iopub.execute_input":"2022-07-08T07:54:51.122733Z","iopub.status.idle":"2022-07-08T07:54:51.149595Z","shell.execute_reply.started":"2022-07-08T07:54:51.122691Z","shell.execute_reply":"2022-07-08T07:54:51.148511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Now we reconstruct our clusters and return an integer column","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"sparse_matrix_mean = sparse_matrix_mean.toarray()\nsparse_matrix_mean","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:51.151361Z","iopub.execute_input":"2022-07-08T07:54:51.152110Z","iopub.status.idle":"2022-07-08T07:54:51.164441Z","shell.execute_reply.started":"2022-07-08T07:54:51.152059Z","shell.execute_reply":"2022-07-08T07:54:51.163217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Intuition:** If sparse_matrix_mean[node1][node2]=1, then node1 and node2 belong to the same cluster\nDisjoint Set Union(DSU) is useful, with complexity of O(α(N)).","metadata":{}},{"cell_type":"code","source":"clusters_final = np.zeros(len(df))\nclusters_final_next_id = 0\n\nnode_end = len(comp_clusters_list[0])\nedge_list = []  # [(w, fr, to), ...]\nfor fr in range(node_end):\n    for to in range(fr, node_end):\n        w = sparse_matrix_mean[fr][to]\n        if w == 0:\n            continue\n        edge_list.append((w, fr, to))\nedge_list.sort(reverse=True)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:54:51.166159Z","iopub.execute_input":"2022-07-08T07:54:51.166988Z","iopub.status.idle":"2022-07-08T07:54:51.179297Z","shell.execute_reply.started":"2022-07-08T07:54:51.166941Z","shell.execute_reply":"2022-07-08T07:54:51.178141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SIZ_MAX = 18000","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:54:51.180405Z","iopub.execute_input":"2022-07-08T07:54:51.181260Z","iopub.status.idle":"2022-07-08T07:54:51.187932Z","shell.execute_reply.started":"2022-07-08T07:54:51.181224Z","shell.execute_reply":"2022-07-08T07:54:51.187139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Disjoint Set Union(DSU)","metadata":{}},{"cell_type":"code","source":"par = [i for i in range(node_end)]\nsiz = [len(unzipper[i]) for i in range(node_end)]\n\ndef find(x):\n    if par[x] == x: return x\n    par[x] = find(par[x])\n    return par[x]\n\ndef union(x, y):\n    x = find(x)\n    y = find(y)\n    if x == y:\n        return\n    if siz[x] > siz[y]: x, y = y, x\n    par[x] = y\n    siz[y] += siz[x]\n\ndef get_siz(x):\n    x = find(x)\n    return siz[x]\n\nfor w, fr, to in edge_list:\n    if (get_siz(fr)+get_siz(to)) > SIZ_MAX: continue\n    union(fr, to)\n\n    clusters_final = [0]*len(clusters_list[0])\nfor node in range(node_end):\n    cluster_id = find(node)\n    idx_list = unzipper[node]\n    for idx in idx_list:\n        clusters_final[idx] = cluster_id","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:56:04.436978Z","iopub.execute_input":"2022-07-08T07:56:04.437380Z","iopub.status.idle":"2022-07-08T07:56:04.466420Z","shell.execute_reply.started":"2022-07-08T07:56:04.437348Z","shell.execute_reply":"2022-07-08T07:56:04.465389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"compress cluster_final to replace it with an integer starting from 0.","metadata":{}},{"cell_type":"code","source":"zipper = {x: i for i, x in enumerate(sorted(set(clusters_final)))}\nclusters_final = [zipper[x] for x in clusters_final]","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:56:21.648753Z","iopub.execute_input":"2022-07-08T07:56:21.649174Z","iopub.status.idle":"2022-07-08T07:56:21.667910Z","shell.execute_reply.started":"2022-07-08T07:56:21.649141Z","shell.execute_reply":"2022-07-08T07:56:21.666275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis\ndef compare_clusterings(y1, y2, title=''):\n    \"\"\"Show the adjusted rand score and plot the two clusterings in color\"\"\"\n    ars = adjusted_rand_score(y1, y2)\n    n1 = y1.max() + 1\n    n2 = y2.max() + 1\n    argsort = np.argsort(y1*100 + y2) if n1 >= n2 else np.argsort(y2*100 + y1)\n    plt.figure(figsize=(16, 0.5))\n    for i in range(6, 11):\n        plt.scatter(np.arange(len(y1)), np.full_like(y1, i), c=y1[argsort], s=1, cmap='tab10')\n    for i in range(5):\n        plt.scatter(np.arange(len(y2)), np.full_like(y2, i), c=y2[argsort], s=1, cmap='tab10')\n    plt.gca().axis('off')\n    plt.title(f'{title}\\nAdjusted Rand score: {ars:.5f}')\n    plt.savefig(title + '.png', bbox_inches='tight')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:56:25.087521Z","iopub.execute_input":"2022-07-08T07:56:25.088046Z","iopub.status.idle":"2022-07-08T07:57:29.513093Z","shell.execute_reply.started":"2022-07-08T07:56:25.087976Z","shell.execute_reply":"2022-07-08T07:57:29.511848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for clusters in clusters_list: compare_clusterings(np.array(clusters), np.array(clusters_final))","metadata":{"execution":{"iopub.status.busy":"2022-07-08T07:56:25.087521Z","iopub.execute_input":"2022-07-08T07:56:25.088046Z","iopub.status.idle":"2022-07-08T07:57:29.513093Z","shell.execute_reply.started":"2022-07-08T07:56:25.087976Z","shell.execute_reply":"2022-07-08T07:57:29.511848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.read_csv(\"../input/tabular-playground-series-jul-2022/sample_submission.csv\")\nsubmission[\"Predicted\"] = clusters_final\nsubmission.to_csv(\"submission.csv\", index=False)\nsubmission","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-08T07:57:29.515521Z","iopub.execute_input":"2022-07-08T07:57:29.516152Z","iopub.status.idle":"2022-07-08T07:57:29.955933Z","shell.execute_reply.started":"2022-07-08T07:57:29.516105Z","shell.execute_reply":"2022-07-08T07:57:29.954889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"_____","metadata":{}},{"cell_type":"markdown","source":"### Single Function\n##### Let's make a single functions out of all of this - for ease of use","metadata":{}},{"cell_type":"code","source":"# Copy this cell wherever you want without thinking 👌\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import trange\nimport matplotlib.pyplot as plt\nfrom sklearn.cluster import KMeans\nfrom collections import defaultdict\nfrom scipy.sparse import csr_matrix\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.metrics import adjusted_rand_score\n\ndef clustering_ensemble(clusters_list, weights = None, threshold = 0.5, size_max = 18000):\n    \"\"\"\n    Parameters: \n                clusters_list:  List of numpy arrays, representing the cluster_id of each data point\n                weights:        List representing the voting weight, for example = [0.5, 0.5]. default None (uses uniform weights)\n                threshold:      float(0.0, 1.0), the threshold for determining an edge in the sparse matrix, default 0.5\n                size_max:       The maximum size of a cluster, default 18000\n    Returns: \n                clusters_final: An ensemble of the clustering algorithms predictions.\n    \"\"\"\n    cls_tup_list = []\n    for cls_tup in zip(*clusters_list):\n        cls_tup_list.append(cls_tup)\n    zipper = {x: i for i, x in enumerate(sorted(set(cls_tup_list)))}\n    zipped_list = [zipper[x] for x in cls_tup_list]\n    unzipper = defaultdict(set)\n    for idx, cls_tup in enumerate(cls_tup_list):\n        zipped = zipper[cls_tup]\n        unzipper[zipped].add(idx)\n    comp_clusters_list = [[-1]*len(zipper) for _ in range(len(clusters_list))]\n    for clusters, comp_clusters in zip(clusters_list, comp_clusters_list):\n        for i, cluster_i in enumerate(clusters):\n                value = zipped_list[i]\n                comp_clusters[value] = cluster_i                \n    def create_sparse_matrix(clusters):\n        n = len(clusters)\n        row = []\n        col = []\n        data = []\n        for i in trange(n):            \n            for j in range(i + 1, n):\n                if clusters[i] == clusters[j]:\n                    data.append(1)\n                    row.append(i)\n                    col.append(j)\n        return csr_matrix((data, (row, col)), shape=(n, n))\n    sparse_matrix_list = [create_sparse_matrix(comp_clusters) for comp_clusters in comp_clusters_list]\n    if weights is None: weights = [1 / len(sparse_matrix_list) for i in range(len(sparse_matrix_list))]\n    weights = np.asarray(weights) / np.sum(weights)\n    sparse_matrix_mean = (sparse_matrix_list[0] * weights[0] + sparse_matrix_list[1] * weights[1]) / len(sparse_matrix_list)\n    sparse_matrix_mean[sparse_matrix_mean < threshold] = 0\n    sparse_matrix_mean = sparse_matrix_mean.toarray()\n    clusters_final = np.zeros(len(clusters_list[0]))\n    clusters_final_next_id = 0\n    node_end = len(comp_clusters_list[0])\n    edge_list = []  # [(w, fr, to), ...]\n    for fr in range(node_end):\n        for to in range(fr, node_end):\n            w = sparse_matrix_mean[fr][to]\n            if w == 0:\n                continue\n            edge_list.append((w, fr, to))\n    edge_list.sort(reverse=True)\n    par = [i for i in range(node_end)]\n    siz = [len(unzipper[i]) for i in range(node_end)]\n    def find(x):\n        if par[x] == x: return x\n        par[x] = find(par[x])\n        return par[x]\n    def union(x, y):\n        x = find(x)\n        y = find(y)\n        if x == y:\n            return\n        if siz[x] > siz[y]: x, y = y, x\n        par[x] = y\n        siz[y] += siz[x]\n    def get_siz(x):\n        x = find(x)\n        return siz[x]\n    for w, fr, to in edge_list:\n        if (get_siz(fr)+get_siz(to)) > size_max: continue\n        union(fr, to)\n        clusters_final = [0]*len(clusters_list[0])\n    for node in range(node_end):\n        cluster_id = find(node)\n        idx_list = unzipper[node]\n        for idx in idx_list:\n            clusters_final[idx] = cluster_id\n    zipper = {x: i for i, x in enumerate(sorted(set(clusters_final)))}\n    clusters_final = [zipper[x] for x in clusters_final]            \n    return clusters_final","metadata":{"execution":{"iopub.status.busy":"2022-07-08T08:36:41.157420Z","iopub.execute_input":"2022-07-08T08:36:41.157813Z","iopub.status.idle":"2022-07-08T08:36:41.184516Z","shell.execute_reply.started":"2022-07-08T08:36:41.157779Z","shell.execute_reply":"2022-07-08T08:36:41.183680Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Simple usage example","metadata":{}},{"cell_type":"code","source":"clusters_1 = pd.read_csv(\"../input/tps-high-scoring-subs-2/final_df.csv\")[\"Predicted\"]\nclusters_2 = pd.read_csv(\"../input/tps-high-scoring-subs-2/submission_simple_soft_voting.csv\")[\"Predicted\"]\n\nclusters_list = [clusters_1, clusters_2]\n\nensemble_preds = clustering_ensemble(clusters_list, weights = None)\n\nfor clusters in clusters_list: compare_clusterings(np.array(clusters), np.array(ensemble_preds))","metadata":{"execution":{"iopub.status.busy":"2022-07-08T08:36:43.739143Z","iopub.execute_input":"2022-07-08T08:36:43.739545Z","iopub.status.idle":"2022-07-08T08:38:52.029672Z","shell.execute_reply.started":"2022-07-08T08:36:43.739511Z","shell.execute_reply":"2022-07-08T08:38:52.028411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv(\"../input/tabular-playground-series-jul-2022/sample_submission.csv\")\nsubmission[\"Predicted\"] = ensemble_preds\nsubmission.to_csv(\"submission.csv\", index=False)\nsubmission","metadata":{"execution":{"iopub.status.busy":"2022-07-08T08:34:31.319386Z","iopub.execute_input":"2022-07-08T08:34:31.319786Z","iopub.status.idle":"2022-07-08T08:34:31.562484Z","shell.execute_reply.started":"2022-07-08T08:34:31.319737Z","shell.execute_reply":"2022-07-08T08:34:31.561337Z"},"trusted":true},"execution_count":null,"outputs":[]}]}