{"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":"# ___Graph-Based Ensemble Clustering___\n\nIn this notebook, I will build upon some of the other notebooks demonstrating the technique of __ensemble clustering__. I would highly encourage folks to take a look at these notebooks (https://www.kaggle.com/code/ehekatlact/tps2207-ultra-fast-ensemble, https://www.kaggle.com/code/nagsdata/ensemble-clustering-algorithms-tsne-visualization, https://www.kaggle.com/code/thedevastator/how-to-ensemble-clustering-algorithms-updated) for some other versions of enemble clustering and leave them an upvote. \n\nFor those interested in the ensemble clustering literature, I would strongly recommend the Strehl and Ghosh paper that outlines the idea of cluster ensembling (https://dl.acm.org/doi/10.1162/153244303321897735) as well Zhang's survey paper on techniques (https://arxiv.org/pdf/1910.02433.pdf).\n\nFrom the original Strehl and Ghosh paper, cluster ensembling is\n\n```\nthe problem of combining multiple partitionings of a set of objects into a single consolidated clustering without accessing the features or algorithms that determined these partitionings.\n```\n\nSo, ideally, an ensemble clustering should produce a clustering of the data that is 1) Robust to stochastic fluctations in clusterings, as many clustering algorithms are stochastic in nature, 2) produce a clustering solution that is better than any given clustering solution, and 3) can take any clustering methodology into the ensemble. To solve this problem, we will treat the idea of combining the multiple partitions as a __graph partitioning__ problem.","metadata":{}},{"cell_type":"code","source":"!pip install scikit-network","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-07-22T14:53:52.581431Z","iopub.execute_input":"2022-07-22T14:53:52.582160Z","iopub.status.idle":"2022-07-22T14:54:07.832402Z","shell.execute_reply.started":"2022-07-22T14:53:52.582062Z","shell.execute_reply":"2022-07-22T14:54:07.831114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os, random, time, umap, tqdm\nimport numpy as np\nimport pandas as pd\nfrom pandas.api.types import CategoricalDtype\n\nfrom sklearn.cluster import KMeans, OPTICS\nfrom sklearn.mixture import GaussianMixture, BayesianGaussianMixture\nfrom sklearn.metrics import silhouette_score\nfrom sklearn.preprocessing import OrdinalEncoder, MinMaxScaler, StandardScaler, OneHotEncoder, PowerTransformer, QuantileTransformer\nfrom sklearn.model_selection import train_test_split, KFold\nfrom sklearn.experimental import enable_iterative_imputer\nfrom sklearn.impute import IterativeImputer, SimpleImputer, KNNImputer\nfrom sklearn.compose import ColumnTransformer\nfrom category_encoders import MEstimateEncoder\nfrom sknetwork.clustering import Louvain\nfrom sklearn.base import ClusterMixin\n\nimport seaborn as sns\nfrom matplotlib import pyplot as plt\n\nplt.style.use(\"seaborn-whitegrid\")\nplt.rc(\"figure\", autolayout=True)\nplt.rc(\n    \"axes\",\n    labelweight=\"bold\",\n    labelsize=\"large\",\n    titleweight=\"bold\",\n    titlesize=14,\n    titlepad=10,\n)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:07.834518Z","iopub.execute_input":"2022-07-22T14:54:07.834937Z","iopub.status.idle":"2022-07-22T14:54:36.901581Z","shell.execute_reply.started":"2022-07-22T14:54:07.834898Z","shell.execute_reply":"2022-07-22T14:54:36.899683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Helper function to do one-hot encoding\n\ndef one_hot_encode(df):\n    X = df.copy()\n    for colname in X.select_dtypes([\"category\", \"object\", \"int\"]):\n        X = X.join(pd.get_dummies(X[colname], prefix=colname))\n        X = X.drop(colname, axis=1)\n    return X","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:36.904130Z","iopub.execute_input":"2022-07-22T14:54:36.905374Z","iopub.status.idle":"2022-07-22T14:54:36.915018Z","shell.execute_reply.started":"2022-07-22T14:54:36.905316Z","shell.execute_reply":"2022-07-22T14:54:36.913606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Exploratory Data Analysis\n\nFor the exploratory analysis, I will only do an import and some simple visualizations. I would highly reocmmend looking at [this](https://www.kaggle.com/code/javigallego/outliers-eda-clustering-tutorial) and [this](https://www.kaggle.com/code/ashaykatrojwar/eda-pca-bayesian-gaussian-mixture) for more thorough EDAs, that consider things like variable distributions and covariance.","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv(\"../input/tabular-playground-series-jul-2022/data.csv\", index_col='id')","metadata":{"papermill":{"duration":5.607504,"end_time":"2022-04-24T16:44:27.769593","exception":false,"start_time":"2022-04-24T16:44:22.162089","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-22T14:54:36.921587Z","iopub.execute_input":"2022-07-22T14:54:36.923852Z","iopub.status.idle":"2022-07-22T14:54:38.584583Z","shell.execute_reply.started":"2022-07-22T14:54:36.923728Z","shell.execute_reply":"2022-07-22T14:54:38.583233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head()","metadata":{"papermill":{"duration":0.050641,"end_time":"2022-04-24T16:44:27.842607","exception":false,"start_time":"2022-04-24T16:44:27.791966","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2022-07-22T14:54:38.586690Z","iopub.execute_input":"2022-07-22T14:54:38.587488Z","iopub.status.idle":"2022-07-22T14:54:38.629007Z","shell.execute_reply.started":"2022-07-22T14:54:38.587445Z","shell.execute_reply":"2022-07-22T14:54:38.627803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.info()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:38.631435Z","iopub.execute_input":"2022-07-22T14:54:38.631949Z","iopub.status.idle":"2022-07-22T14:54:38.670070Z","shell.execute_reply.started":"2022-07-22T14:54:38.631909Z","shell.execute_reply":"2022-07-22T14:54:38.668670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# look at the continuous variables\ndf.select_dtypes(include=[\"float\"]).describe().T.style.background_gradient(cmap='Blues')","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:38.672834Z","iopub.execute_input":"2022-07-22T14:54:38.673876Z","iopub.status.idle":"2022-07-22T14:54:39.084242Z","shell.execute_reply.started":"2022-07-22T14:54:38.673819Z","shell.execute_reply":"2022-07-22T14:54:39.078685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_distributions(df):\n    '''\n    code from https://www.kaggle.com/code/javigallego/outliers-eda-clustering-tutorial\n    '''\n    figure = plt.figure(figsize = (16,8))\n    for i in range(29):\n        feature_name = 'f_0{}'.format(i) if i < 10 else 'f_{}'.format(i) \n        plt.subplot(5, 6, i+1)\n        if df[feature_name].dtype == 'int': \n            sns.kdeplot(df[feature_name], fill = True, color = 'red')\n        else: \n            sns.kdeplot(df[feature_name], fill = True)        \n    figure.tight_layout(h_pad=1.0, w_pad=0.5)\n    plt.suptitle('Distribution Plots', y=1.02)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:39.087052Z","iopub.execute_input":"2022-07-22T14:54:39.087703Z","iopub.status.idle":"2022-07-22T14:54:39.103894Z","shell.execute_reply.started":"2022-07-22T14:54:39.087636Z","shell.execute_reply":"2022-07-22T14:54:39.102110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_distributions(df)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:39.106185Z","iopub.execute_input":"2022-07-22T14:54:39.107089Z","iopub.status.idle":"2022-07-22T14:54:56.752869Z","shell.execute_reply.started":"2022-07-22T14:54:39.107030Z","shell.execute_reply":"2022-07-22T14:54:56.751914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" ### Some observations\n - data consists of two group of continuous variable columns and a group of integer variable columns; integer variable columns may be categories\n - no missing data or other unusual data\n - correlations exist between 07 through 13, 22 through 28, 07-13 with 22-28, and a weak one between 19 and 07","metadata":{}},{"cell_type":"markdown","source":"# 2. Create Ensemble Members\n\nAn important first step to create a cluster ensembling is to create ensemble members. To do this, I first build upon previous work showing the utility of Gaussian Mixture Models for this challenge, [here](https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis), [here](https://www.kaggle.com/code/ashaykatrojwar/eda-pca-bayesian-gaussian-mixture), and [here](https://www.kaggle.com/code/azminetoushikwasi/different-clustering-techniques-algorithms). These analysis showed that 1) GMM's were good clustering methods for this data (to include just subsetting the data for [certain columns](https://www.kaggle.com/code/akioonodera/tps-jul2022-bgmm), 2) there are multiple ways to transform the data, and 3) that there are around 7 clusters in the data. So, to build the ensemble, I will use all of these aspects.","metadata":{}},{"cell_type":"code","source":"# set up possible cluster numbers around 7, along with several at 7 for the stochatsic clustering algorithms\n\nn_cluster_possibilities = [6,7,7,7,7,7,7,7,8,9]","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:54:56.755716Z","iopub.execute_input":"2022-07-22T14:54:56.756374Z","iopub.status.idle":"2022-07-22T14:54:56.760712Z","shell.execute_reply.started":"2022-07-22T14:54:56.756334Z","shell.execute_reply":"2022-07-22T14:54:56.759734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# set up data transformations, with different scalings of the data and subsets of the data\n\nscaled_all_features_df = pd.DataFrame(StandardScaler().fit_transform(df), columns=df.columns)\n\ntwo_block_columns = ['f_07', 'f_08', 'f_09', 'f_10', 'f_11', 'f_12', 'f_13']\none_hot_two_columns = one_hot_encode(df[two_block_columns])\npower_scaled_two_columns = pd.DataFrame(PowerTransformer().fit_transform(df[two_block_columns]), columns=two_block_columns)\nnon_two_scaled_features_df = scaled_all_features_df[scaled_all_features_df.columns[~scaled_all_features_df.columns.isin(['f_07', 'f_08', 'f_09', 'f_10', 'f_11', 'f_12', 'f_13'])]]\n\nscaled_and_one_hot_two_block_df = non_two_scaled_features_df.join(one_hot_two_columns)\nscaled_and_power_two_block_df = non_two_scaled_features_df.join(power_scaled_two_columns)\nmost_useful_scaled_and_two_block_power_df = scaled_and_power_two_block_df[['f_07', 'f_08', 'f_09', 'f_10', 'f_11', 'f_12', 'f_13', 'f_22', 'f_23', 'f_24', 'f_25', 'f_26', 'f_27', 'f_28']]\n\n\ndata_transformations = {\"standard\":scaled_all_features_df, \n                        \"one_hot\":scaled_and_one_hot_two_block_df, \n                        \"power\":scaled_and_power_two_block_df, \n                        \"most_useful\": most_useful_scaled_and_two_block_power_df \n                       }\n","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:59:26.000584Z","iopub.execute_input":"2022-07-22T14:59:26.001030Z","iopub.status.idle":"2022-07-22T14:59:27.202792Z","shell.execute_reply.started":"2022-07-22T14:59:26.000986Z","shell.execute_reply":"2022-07-22T14:59:27.201715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Try some different clustering models, across all the possible cluster numbers and data transformations\n\nclusterings = pd.DataFrame(index=df.index)\n\nfor data_label, data in data_transformations.items():\n    for n_clusters in tqdm.tqdm(n_cluster_possibilities):\n        model = GaussianMixture(n_components=n_clusters, n_init=3, max_iter=400)\n        clusterings[\"gmm_{}_{}\".format(n_clusters, data_label)] = model.fit_predict(data)\n        model = BayesianGaussianMixture(n_components=n_clusters, n_init=3, max_iter=400)\n        clusterings[\"bgmm_{}_{}\".format(n_clusters, data_label)] = model.fit_predict(data)","metadata":{"execution":{"iopub.status.busy":"2022-07-22T14:57:47.275204Z","iopub.execute_input":"2022-07-22T14:57:47.275624Z","iopub.status.idle":"2022-07-22T14:58:51.488625Z","shell.execute_reply.started":"2022-07-22T14:57:47.275591Z","shell.execute_reply":"2022-07-22T14:58:51.487051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusterings.to_csv(\"ensemble_of_clusterings.csv\", index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Do Ensemble Clustering and Submit Solutions\n\nTo do a Graph-based ensembling, I will demonstrate two technique. The first is  __B__ipartite __G__raph __P__artitioning __A__lgorithm, BGPA, which is from the following [paper](https://dl.acm.org/doi/10.1145/1015330.1015414). This technique creates a bipartite graph between objects and the cluster labels, across all of the base clusterings, and then cluster that graph. Note that in constrast to the original paper, I use BiLouvain from [sci-kit network](https://scikit-network.readthedocs.io/en/latest/tutorials/clustering/louvain.html) to cluster the object-by-cluster bipartite graph. The second technique __L__ocally __W__eighted __B__ipartite __G__raph Partitioning Algorithm, LWBG, works much the same way except that it applirs a weighting term to edges in the bipartite graph to form a better graph to cring out the clusters. Its paper can be found [here](https://arxiv.org/pdf/1605.05011.pdf)","metadata":{}},{"cell_type":"code","source":"clusterings = pd.read_csv(\"ensemble_of_clusterings.csv\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class BGPA(ClusterMixin):\n    \"\"\"Bipartite Graph Partitioning Algorithm\n    \n    BGPA+ : treat the objects-by-bse clsuterings\n    as a bipartite graph partitioning problem\n    \n    Parameters\n    ----------\n    base_clusters : list of array_like or pandas dataframe\n        List of cluster labels where each set of labels must be in the same order\n        and have shape (n_objects,). Note: standard output from sci-kit learn \n        clustering algorithms will return this format. If pandas dataframe, \n        should be shape (n_objects, n_clusters)\n    \n    Returns\n    -------\n    z : ndarray\n        cluster labels of shape (n,).\n    \"\"\"\n    \n    \n    def __init__(self, metaclustering_alg = 'louvain', n_clusters = (None,None),\n                 refine_clusters = False, ensemble_iterations=10):\n        \n        self.metaclustering_alg = metaclustering_alg\n        self.n_clusters = n_clusters\n        self.refine_clusters = refine_clusters\n        self.ensemble_iterations = ensemble_iterations\n    \n    def fit_predict(self, base_clusters):\n        if self.refine_clusters:\n            #note: using a refinement process in BGPA is still in development\n            #for now, its best to not use it.\n            converged = False\n            ba_matrix = self._create_ba_matrix(base_clusters)\n            while not converged:\n                final_clusterer = meta_alg()\n                clusters = [final_clusterer.cluster(ba_matrix, self.n_clusters, self.metaclustering_alg)\n                            for _ in range(self.ensemble_iterations)]\n                ba_matrix = self._create_ba_matrix(clusters)\n                C = ba_matrix @ ba_matrix.T\n                if np.all(C[np.nonzero(C)] >= len(base_clusters)):\n                    converged = True\n            return clusters\n            \n        else:\n            ba_matrix = self._create_ba_matrix(base_clusters)\n            final_clusterer = Louvain()\n            final_clusterer.fit(ba_matrix)\n            self.obj_clusters, self.cluster_clusters = final_clusterer.labels_row_, final_clusterer.labels_col_\n            return self.obj_clusters\n        \n    def _create_ba_matrix(self, clusters):\n        '''\n        Internal function for converting list of clustering labels into a\n        binary-association matrix that is (object,clusters)\n        '''\n        \n        if isinstance(clusters, pd.DataFrame):\n            clusters = [clusters.loc[:,i].to_numpy(dtype=int) for i in clusters.columns]\n        else:\n            clusters = [clustering.astype(np.int64) for clustering in clusters]\n        \n        ba_matrices =[]\n        for base_cluster in clusters:\n            ba_matrix = np.zeros((base_cluster.size, base_cluster.max()+1))\n            ba_matrix[np.arange(base_cluster.size), base_cluster] = 1\n            ba_matrices.append(ba_matrix)\n            \n        return np.concatenate(ba_matrices, axis=1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class LWBG(ClusterMixin):\n    \"\"\"Locally Weighted Bipartite Graph Partitioning Algorithm \n    \n    LWBG : treat the objects-by-base clsuters as a bipartite graph \n    partitioning problem. In this case the bipartite graph is weighted by \n    information-thoeretic measures\n    \n    Parameters\n    ----------\n    base_clusters : list of array_like or pandas dataframe\n        List of cluster labels where each set of labels must be in the same order\n        and have shape (n_objects,). Note: standard output from sci-kit learn \n        clustering algorithms will return this format. If pandas dataframe, \n        should be shape (n_objects, n_clusters)\n    theta : float\n        Controls the impact of the local weighting. Default is 0.5\n    \n    Returns\n    -------\n    z : ndarray\n        cluster labels of shape (n,).\n    \"\"\"\n    \n    def __init__(self, theta = 0.5):\n        self.theta = theta\n    \n    def fit_predict(self, base_clusters):\n        ba_matrices =[]\n        idxs = []\n        idxs_iter = 0\n        \n        if isinstance(base_clusters, pd.DataFrame):\n            base_clusters = [base_clusters.loc[:,i].to_numpy(dtype=int) for i in base_clusters.columns]\n        else:\n            base_clusters = [clustering.astype(np.int64) for clustering in base_clusters]\n        \n        for base_cluster in base_clusters:\n            ba_matrix = np.zeros((base_cluster.size, base_cluster.max()+1))\n            ba_matrix[np.arange(base_cluster.size), base_cluster] = 1\n            ba_matrices.append(ba_matrix)\n            idxs.append([idxs_iter, idxs_iter+ba_matrix.shape[1]])\n            idxs_iter += ba_matrix.shape[1]\n            \n        ba_matrix = np.concatenate(ba_matrices, axis=1)\n        \n        H_matrix = np.zeros((ba_matrix.shape[1], len(ba_matrices)))\n        \n        for col_idx in range(ba_matrix.shape[1]):\n            for mode in range(len(idxs)):\n                idx = idxs[mode]\n                if idx[0] <= col_idx <=idx[1]:\n                    pass\n                else:\n                    H_m = []\n                    for alt_idx in range(idx[0], idx[1]):\n                        C_i = ba_matrix[:,col_idx]\n                        C_j = ba_matrix[:,alt_idx]\n                        P_ij = np.dot(C_i, C_j)/ np.count_nonzero(C_i)\n                        if P_ij > 0:\n                            H_m.append(P_ij * np.log2(P_ij))\n                        else:\n                            H_m.append(0)\n                    H_matrix[col_idx, mode] = -1*np.sum(H_m)\n        \n        ECI = np.exp(-1* np.sum(H_matrix, axis=1)/(self.theta*len(ba_matrices)))\n        weighted_ba_matrix = ba_matrix * ECI\n        final_clusterer = Louvain()\n        final_clusterer.fit(weighted_ba_matrix)\n        self.obj_clusters, self.cluster_clusters = final_clusterer.labels_row_, final_clusterer.labels_col_\n        return self.obj_clusters","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ensemble_clustering_method = LWBG()\n\nensemble_clusters = ensemble_clustering_method.fit_predict(clusterings)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Take a look at how many clusters we end up with, and how many objects per cluster\n\nnp.unique(ensemble_clusters, return_counts=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(clusterings.index.values, columns=[\"Id\"])\nsubmission['Predicted'] = ensemble_clusters","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"submission.csv\", index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As always, I hope you have enjoyed this notebook and feel free to use any of the code or techniques here in your own work. I welcome an comments and would appreciate an upvote!","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}