{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Uncover hidden image secrets! \n\nI used this notebook to generate datasets for the image statistics for the [SIIM-ISIC Melanoma Classification competition](https://www.kaggle.com/c/siim-isic-melanoma-classification). It's related to the EDA that I have done to highlight the differences and challenges in train and test data in my notebook [Don't turn into a Smoothie after the Shake-Up](https://www.kaggle.com/allunia/don-t-turn-into-a-smoothie-after-the-shake-up). \n\nThe intend of this notebook is:\n\n* to provide the code for generating image statistics\n* to extend the EDA to gain more insights \n\nHave fun! :-)\n\n\n\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"You are probably using external datasets and like to redo or complete the computation of image statistics. I have been asked to share my notebook for data generation. So here it is. ;-) \n\nIt's very likely that I will work on this notebook in a few days to expand the EDA. Stay tuned! ;-)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Table of contents\n\n1. [Prepare to start](#prepare)\n2. [Generating image statistics data](#data_generation) \n3. [About the differences of dicom and jpeg](#dcmjpg)\n4. [Explore the insights](#EDA) \n * [The image shape matters!](#image_shape)\n * [Open your eyes with mean and std](#open_eyes) \n * [Can we group images by shape?](#shape_groups) \n5. [Clustering based on image statistics](#stats_groups) \n * [Starting easy with K-Means](#k_means)\n * [Why should we bother about \"how K-Means works\"?](#preprocessing)\n * [How many clusters should we choose?](#num_clusters)\n * [Running K-means](#run_kmeans)\n * [The cluster patchwork quilt](#patchwork_clusters)\n6. [Clustering with GMM](#gmm)\n * [Why are GMM and Kmeans similar but different?](#gmm_kmeans)\n * [Running GMM on image statistics](#run_gmm)\n * [Exploring clusters](#cluster_eda)\n * [Exploring anomalies](#anomalies)\n7. [Fitting catboost and submission](#catboost)\n * [Data preparation](#data_prep)\n * [Validation strategy](#validation)\n * [Fitting](#fitting)\n * [Feature importances](#shap_values)\n * [Submission](#submission)\n8. [Conclusion](#conclusion)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Prepare to start ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Loading packages...","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport seaborn as sns\nsns.set()\n\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objects as go\n\nfrom PIL import Image\nimport pydicom\nfrom skimage.io import imread\n\nfrom sklearn.preprocessing import MinMaxScaler, StandardScaler\nfrom sklearn.cluster import KMeans\nfrom sklearn.mixture import GaussianMixture\nfrom scipy.stats import boxcox\n\nfrom catboost import CatBoostClassifier, Pool, cv\n\nimport shap\n# load JS visualization code to notebook\nshap.initjs()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What kind of datasets are given?","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from os import listdir\nlistdir(\"../input/\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To compute the statistics we need the dicom folder:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"basepath = \"../input/siim-isic-melanoma-classification/\"\ntrain_image_path = basepath + \"/train/\"\nsome_files = listdir(train_image_path)[0:5]\nsome_files","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info = pd.read_csv(basepath + \"train.csv\")\ntrain_info.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_info = pd.read_csv(basepath + \"test.csv\")\ntest_info.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info[\"dcm_path\"] = basepath + \"/train/\" + train_info.image_name + \".dcm\"\ntest_info[\"dcm_path\"] = basepath + \"/test/\" + test_info.image_name + \".dcm\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info.head(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_info.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generating image statistics data \n\nIf you like to redo the computation of image statistics, here is the code. Otherwise you could also load the [dataset](https://www.kaggle.com/allunia/siimisic-melanoma-classification-image-stats) I have published. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"extract_train_0 = False\nextract_train_1 = False\nextract_train_2 = False\nextract_test = False","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"from tqdm.notebook import tqdm \nfrom scipy.stats import skew\n\n\ndef extract_shapes(df):\n all_paths = df.dcm_path.values\n image_eda = pd.DataFrame(index=np.arange(len(df)),\n columns=[\"path\", \"rows\", \"columns\",\n \"channels\", \"img_mean\", \"img_std\",\n \"img_skew\", \"red_mean\", \"green_mean\",\n \"blue_mean\"])\n for i in tqdm(range(0, len(df))):\n path = all_paths[i]\n\n dcm_file = pydicom.dcmread(path)\n image = dcm_file.pixel_array\n\n image_eda.iloc[i][\"path\"] = path\n image_eda.iloc[i][\"rows\"] = image.shape[0]\n image_eda.iloc[i][\"columns\"] = image.shape[1]\n image_eda.iloc[i][\"channels\"] = image.shape[2]\n\n # some image stats\n image_eda.iloc[i][\"img_mean\"] = np.mean(image.flatten())\n image_eda.iloc[i][\"img_std\"] = np.std(image.flatten())\n image_eda.iloc[i][\"img_skew\"] = skew(image.flatten())\n\n image_eda.iloc[i][\"red_mean\"] = np.mean(image[:,:,0].flatten())\n image_eda.iloc[i][\"green_mean\"] = np.mean(image[:,:,1].flatten())\n image_eda.iloc[i][\"blue_mean\"] = np.mean(image[:,:,2].flatten())\n\n return image_eda\n\n\nif extract_train_0:\n train_shapes = extract_shapes(train_info.iloc[0:10000])\n train_shapes.to_csv(\"train_image_stats_10000.csv\", index=False)\nelif extract_train_1:\n train_shapes = extract_shapes(train_info.iloc[10000:20000])\n train_shapes.to_csv(\"train_image_stats_20000.csv\", index=False)\nelif extract_train_2:\n train_shapes = extract_shapes(train_info.iloc[20000::])\n train_shapes.to_csv(\"train_image_stats_toend.csv\", index=False)\nelif extract_test:\n test_shapes = extract_shapes(test_info)\n test_shapes.to_csv(\"test_image_stats.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loading stats and meta features for train and test\n\nI have published a [dataset](https://www.kaggle.com/allunia/siimisic-melanoma-classification-image-stats) that contains the image statistics on jpeg-images and meta features for the original train and test data in one csv-file. Let's load these files instead of generating new one. You can also load the other files but they were computed using dicom!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"listdir(\"../input/siimisic-melanoma-classification-image-stats\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats = pd.read_csv(\"../input/siimisic-melanoma-classification-image-stats/train_stats_and_meta.csv\")\ntest_image_stats = pd.read_csv(\"../input/siimisic-melanoma-classification-image-stats/test_stats_and_meta.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_names = train_image_stats.image_name.values\ntest_image_names = test_image_stats.image_name.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats[\"img_area\"] = train_image_stats[\"rows\"] * train_image_stats[\"columns\"]\ntest_image_stats[\"img_area\"] = test_image_stats[\"rows\"] * test_image_stats[\"columns\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_image_stats.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# About the differences of dicom and jpeg \n\n","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"jpeg_path = basepath + \"/jpeg/train/\"\n\nfig, ax = plt.subplots(3,2,figsize=(20,15))\nfor n in range(3):\n dcm_file = pydicom.dcmread(train_info.dcm_path.values[n])\n pixel_array = dcm_file.pixel_array\n \n image_path = jpeg_path + train_info.image_name.values[n] + \".jpg\"\n image = imread(image_path)\n \n sns.distplot(pixel_array.flatten(), ax=ax[n,0], color=\"sienna\")\n ax[n,0].set_title(\"Distribution of values in dicom pixelarrays\")\n \n sns.distplot(image.flatten(), ax=ax[n,1], color=\"deepskyblue\")\n ax[n,1].set_title(\"Distribution of values in jpeg images\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this notebook I'm using statistics on jpeg files!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Explore the insights! \n\nLet's do some more EDA to gain more and more insights about the training and test images. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## The image shape matters! \n\nIn the first notebook I wrote one can already see that there are some image groups in the test data that differ from the training images and one of this group holds ~15 % in test! These groups can be formed and summarized by the image shape. Images of the same shape show similar statistics. \n\nTo uncover more groups, let's take a look at a scatter plot of row and column values in train and test:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nax[0].scatter(train_image_stats[\"rows\"].values, train_image_stats[\"columns\"].values, c=\"orangered\")\nax[1].scatter(test_image_stats[\"rows\"].values, test_image_stats[\"columns\"].values, c=\"lightseagreen\")\n\nax[0].set_title(\"Train images\")\nax[1].set_title(\"Test images\")\n\nfor n in range(2):\n ax[n].set_xlabel(\"row value\")\n ax[n].set_ylabel(\"column value\")\n ax[n].set_xlim([0,6500])\n ax[n].set_ylim([0,6500])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\nCan you see it? There are two very interesting insights here:\n\n1. There is one very clear linear straight line of image shapes in train and test and 2-3 more could be imagined as well. Is there some relationship between images with shapes on the line? Can we use simple rescaling with images on the line or do they still differ too much in their color information?\n2. There are several groups in train that are not present in test and vice versa there are also some groups in test that are not given in train!\n\nLet's overlap the train spots with the test spots to see groups in test that are not given in train:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.scatter(test_image_stats[\"rows\"].values, test_image_stats[\"columns\"].values, c=\"lightseagreen\")\nplt.scatter(train_image_stats[\"rows\"].values, train_image_stats[\"columns\"].values, c=\"orangered\");\nplt.title(\"Uncover groups in test that are not present in train\");\nplt.xlabel(\"row value\")\nplt.ylabel(\"column value\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are at least 4 groups that are not present in train and do not live on the straight line.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Open your eyes with mean and std ","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"from plotly.subplots import make_subplots\nfig = make_subplots(rows=1, cols=2, subplot_titles=(\"Train image stats\", \"Test image stats\"))\n\n\ntrace0 = go.Scatter(\n x = train_image_stats.img_std.values,\n y = train_image_stats.img_mean.values,\n mode='markers',\n text=train_image_stats[\"rows\"].values,\n marker=dict(\n color=train_image_stats[\"columns\"].values,\n colorscale='Jet',\n opacity=0.4,\n size=4\n )\n)\n\ntrace1 =go.Scatter(\n x = test_image_stats.img_std,\n y = test_image_stats.img_mean,\n mode='markers',\n text=test_image_stats[\"rows\"],\n marker=dict(\n color=test_image_stats[\"columns\"],\n colorscale='Jet',\n colorbar=dict(thickness=10, len=1.1, title=\"image columns\"),\n opacity=0.4,\n size=4\n )\n)\n\nfig.add_trace(trace0, row=1, col=1)\nfig.add_trace(trace1, row=1, col=2)\n\nfig.update_xaxes(title_text=\"Image std\", row=1, col=1)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=1)\nfig.update_xaxes(title_text=\"Image std\", row=1, col=2)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=2)\n\nfig.update_layout(height=425, width=850, showlegend=False)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Plotting the image mean and standard deviation for each image in train and test we can easily see that there is a big group in test (~140 mean, ~60 std) that is not present in train. This is the group that holds ~15% of the test data.\n* Furthermore we can see that the mean and standard deviation depends on the row and column value. There are clear clusters and groups that are related to these quantities. The columns are given by the color and the rows by the text you can see when browsing over a spot. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Can we group images by shape? \n\nWe have seen that there are several groups of images that show different image statistics like mean and standard deviation. As you can see by the colors and texts of the scatter plot above, these groups also depend on the shape of the image. A first good idea is therefore to group by shape. But some images have uncommon shapes and they need to be assigned to a similar group that shows same kind of image information. For this reason a clustering could help us more than just grouping by shape. ;-)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How often can we find images whose shape is unique and not part of a bigger, more common group?","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(train_image_stats[\"rows\"].value_counts().values,\n kde=False, ax=ax[0], color=\"magenta\")\nsns.distplot(train_image_stats[\"columns\"].value_counts().values,\n kde=False, ax=ax[1], color=\"mediumvioletred\")\nax[0].set_title(\"Distribution of row value counts\")\nax[1].set_title(\"Distribution of column value counts\")\nfor n in range(2):\n ax[0].set_xlabel(\"value counts\")\n ax[0].set_ylabel(\"frequency\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"By looking at the value counts of row and column values we can see that most values are only present once. In contrast groups with a few hundereds to thousand images having the same row or column value are much smaller! If you are not sure how to interpret these plots take a look at some examples instead:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats[\"rows\"].value_counts().sort_values().iloc[0:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is only one image with a row value of 1955. In contrast there are 14703 images with a row value of 4000:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats[\"rows\"].value_counts().sort_values(ascending=False).iloc[0:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If we would like to work with image groups we need to find a better method for grouping them then just using the shape as it is.","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(train_image_stats.img_area, ax=ax[0], kde=False, bins=30, color=\"orangered\")\nsns.distplot(test_image_stats.img_area, ax=ax[1], kde=False, bins=30, color=\"lightseagreen\")\n\nfor n in range(2):\n ax[n].set_xlabel(\"image area\")\n ax[n].set_ylabel(\"counts\")\n\nax[0].set_title(\"Image area distribution in train\")\nax[1].set_title(\"Image area distribution in test\");","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,3,figsize=(20,12))\n\nsns.distplot(train_image_stats.img_mean, ax=ax[0,0], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.img_mean, ax=ax[0,0], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.img_std, ax=ax[0,1], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.img_std, ax=ax[0,1], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.img_skew, ax=ax[0,2], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.img_skew, ax=ax[0,2], color=\"lightseagreen\", label=\"test\")\n\nsns.distplot(train_image_stats.red_mean, ax=ax[1,0], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.red_mean, ax=ax[1,0], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.green_mean, ax=ax[1,1], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.green_mean, ax=ax[1,1], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.blue_mean, ax=ax[1,2], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.blue_mean, ax=ax[1,2], color=\"lightseagreen\", label=\"test\")\n\n\nfor n in range(3):\n for m in range(2):\n ax[m,n].set_ylabel(\"density\");\n ax[m,n].legend()\n \nax[0,0].set_title(\"Image means\")\nax[0,1].set_title(\"Image stds\")\nax[0,2].set_title(\"Image skewnesses\");\n\nax[1,0].set_title(\"Red channel mean\")\nax[1,1].set_title(\"Green channel mean\")\nax[1,2].set_title(\"Blue channel mean\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* All distributions are skewed somehow and show outliers. We need to perform some preprocessing to prevent a bad influence on the cluster center computation.\n* We can clearly see that there are differences in train and test image statistics.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"With preprocessing we should try to transform our distributions such that they look more normal. One way to do this is the [boxcox-transformation](https://en.wikipedia.org/wiki/Power_transform). The two parameters constant and lambda have to be set carefully and one needs to play a bit with them to find a good choice:","execution_count":null},{"metadata":{"_kg_hide-input":false,"trusted":true,"_kg_hide-output":false},"cell_type":"code","source":"def preprocess_k_means(train_data, test_data, feature, constant, lam):\n min_max_scaler = MinMaxScaler()\n scaled_train_feature = min_max_scaler.fit_transform(train_data[feature].values.reshape(-1, 1))\n scaled_test_feature = min_max_scaler.transform(test_data[feature].values.reshape(-1,1))\n \n boxcox_train_feature = boxcox(scaled_train_feature[:,0] + constant, lam)\n boxcox_test_feature = boxcox(scaled_test_feature[:,0] + constant, lam)\n\n scaler = StandardScaler()\n preprocessed_train_feature = scaler.fit_transform(boxcox_train_feature.reshape(-1,1))\n preprocessed_test_feature = scaler.fit_transform(boxcox_test_feature.reshape(-1,1))\n \n train_data.loc[:, \"preprocessed_\" + feature] = preprocessed_train_feature\n test_data.loc[:, \"preprocessed_\" + feature] = preprocessed_test_feature\n return train_data, test_data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's preprocess the skewness as well as the red and green mean:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats, test_image_stats = preprocess_k_means(train_image_stats,\n test_image_stats, \n \"red_mean\",\n constant=1, lam=10)\ntrain_image_stats, test_image_stats = preprocess_k_means(train_image_stats,\n test_image_stats, \n \"green_mean\",\n constant=0.5,lam=2)\ntrain_image_stats, test_image_stats = preprocess_k_means(train_image_stats,\n test_image_stats, \n \"img_skew\",\n constant=0.05,lam=2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The preprocessed distributions are less skewed now and the effect of outliers is suppressed:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,5))\n\nsns.distplot(train_image_stats.preprocessed_red_mean, ax=ax[0], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.preprocessed_red_mean, ax=ax[0], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.preprocessed_green_mean, ax=ax[1], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.preprocessed_green_mean, ax=ax[1], color=\"lightseagreen\", label=\"test\")\nsns.distplot(train_image_stats.preprocessed_img_skew, ax=ax[2], color=\"crimson\", label=\"train\")\nsns.distplot(test_image_stats.preprocessed_img_skew, ax=ax[2], color=\"lightseagreen\", label=\"test\")\n\n\nfor n in range(3):\n ax[n].set_ylabel(\"density\");\n ax[n].legend()\n \nax[0].set_title(\"Image means\")\nax[1].set_title(\"Image stds\")\nax[2].set_title(\"Image skewnesses\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's take a look at all basic statistic features at once in a scatterplot:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig = make_subplots(rows=1, cols=2, subplot_titles=(\"Preprocessed train image stats\", \"Preprocessed test image stats\"))\n\n\ntrace0 = go.Scatter(\n x = train_image_stats.img_mean.values,\n y = train_image_stats.img_std.values,\n mode='markers',\n text=train_image_stats[\"columns\"].values,\n marker=dict(\n color=train_image_stats.preprocessed_img_skew.values,\n colorbar=dict(thickness=10, len=1.1, title=\"preprocessed skewness\"),\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\ntrace1 = go.Scatter(\n x = test_image_stats.img_mean.values,\n y = test_image_stats.img_std.values,\n mode='markers',\n text=test_image_stats[\"columns\"].values,\n marker=dict(\n color=test_image_stats.preprocessed_img_skew.values,\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\nfig.add_trace(trace0, row=1, col=1)\nfig.add_trace(trace1, row=1, col=2)\n\nfig.update_xaxes(title_text=\"Image std\", row=1, col=1)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=1)\nfig.update_xaxes(title_text=\"Image std\", row=1, col=2)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=2)\n\nfig.update_layout(height=425, width=850, showlegend=False)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can observe that the differences in the train and test distributions with mean values (70) and std (150) belongs to a hidden group of images!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## How many clusters should we choose? \n \nWe have already found that the image statistics depend on the shape of the images. For this reason we could set the number of clusters equal to the number of the most common groups find in train and test:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_shapes = train_image_stats.groupby(\n [\"rows\", \"columns\"]).size().sort_values(ascending=False) / train_image_stats.shape[0] * 100\ntest_shapes = test_image_stats.groupby(\n [\"rows\", \"columns\"]).size().sort_values(ascending=False) / test_image_stats.shape[0] * 100","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's take all groups into account that hold at least 0.2 % of all images. For example for the training data this would mean at least this number of images in one cluster:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.shape[0] * 0.2/100","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"common_train_shapes = set(list(train_shapes[train_shapes > 0.3].index.values))\ncommon_test_shapes = set(list(test_shapes[test_shapes > 0.3].index.values))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"common_shape_groups = common_train_shapes.union(common_test_shapes)\ncommon_shape_groups","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_clusters = len(common_shape_groups)\nnum_clusters","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Running K-Means ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"combined_stats = train_image_stats.append(test_image_stats)\ncombined_stats.head(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kmeans = KMeans(n_clusters=num_clusters, \n random_state=0)\n\nx = combined_stats.loc[:, [\"img_mean\", \"img_std\", \"preprocessed_img_skew\",\n \"preprocessed_red_mean\", \"preprocessed_green_mean\", \"blue_mean\"]].values #,\n #\"img_area\", \"rows\", \"columns\"]].values\ncluster_labels = kmeans.fit_predict(x)\ncombined_stats[\"cluster_label\"] = cluster_labels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats = combined_stats.iloc[0:train_image_stats.shape[0]]\ntest_image_stats = combined_stats.iloc[train_image_stats.shape[0]::]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig = make_subplots(rows=1, cols=2, subplot_titles=(\"Train image stats\", \"Test image stats\"))\n\n\ntrace0 = go.Scatter(\n x = train_image_stats.img_std.values,\n y = train_image_stats.img_mean.values,\n mode='markers',\n text=train_image_stats[\"cluster_label\"].values,\n marker=dict(\n color=train_image_stats.cluster_label.values,\n colorbar=dict(thickness=10, len=1.1, title=\"cluster label\"),\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\ntrace1 = go.Scatter(\n x = test_image_stats.img_std.values,\n y = test_image_stats.img_mean.values,\n mode='markers',\n text=test_image_stats[\"cluster_label\"].values,\n marker=dict(\n color=test_image_stats.cluster_label.values,\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\nfig.add_trace(trace0, row=1, col=1)\nfig.add_trace(trace1, row=1, col=2)\n\nfig.update_xaxes(title_text=\"Image std\", row=1, col=1)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=1)\nfig.update_xaxes(title_text=\"Image std\", row=1, col=2)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=2)\n\nfig.update_layout(height=425, width=850, showlegend=False)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* We have obtained a patchwork quilt of cluster points. This is a drawback of K-Means. It's not able to adapt well to the shape of elliptical clusters. This won't improve with further preprocessing and we need a better clustering algorithm to cover the shapes better! ;-)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"combined_stats = train_image_stats.append(test_image_stats)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## The cluster patchwork quilt ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at some example images within these 8 patchwork clusters. ","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"fig, ax = plt.subplots(num_clusters,8, figsize=(20, 2.5*num_clusters))\n\nfor cluster in range(num_clusters):\n selection = np.random.choice(combined_stats[combined_stats.cluster_label==cluster].image_path.values,\n size=8, replace=False)\n m=0\n for path in selection:\n image = imread(path)\n ax[cluster, m].imshow(image)\n ax[cluster, m].set_title(\"K-Means cluster {}\".format(cluster))\n ax[cluster, m].axis(\"off\")\n m+=1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n1. With K-Means we already found some cluster that contain similar images. \n3. Furthermore some images seem to fit not so well in their group.\n\nIn the next step I like to improve the clustering using a Gaussian Mixture Model. It has the advantage that it can adapt better to the elliptical cluster shapes and it's also powerful to detect anomalies in the data. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Clustering with GMM ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Why are GMM and Kmeans similar but different? \n\nBy clustering with gaussian mixture models we still follow an interative approach of updating cluster centers and assigning data spots to them but in a so called \"soft\" manner. This means now that each cluster is responsible for each data point but with different strenghts. The cluster that holds the highest responsibility for one data point can then be hard assigned as well. \n\nIn the case of **mixture models it's assumed that there exist some hidden latent variables that generate the distribution of data you observe**. In our case we hope that these latent variables walk along with different cameras that have taken the images at different conditions. Some cameras might always produce dark, small images whereas other ones might take large, but bright images. Do you get the idea? We don't know which camera was used and we don't know how the conditions looked like. This information is hidden. \n\n\nNow the gaussian mixture model descripes the following: **Each of our latent variable, each camera, causes one multivariate gaussian**, this is a normal distribution over more than one dimension. In our case we like to use the statistics per image like mean, standard deviation and skewness. Doing so our feature space is 3-dimensional and by looking at the 3d-scatterplot we can observe ellipsoids with different densities, locations and expansions. Then our distribution $p(x)$ of data we observe is said to be generated by these gaussian ellipsoids:\n\n$$ p(x) = \\sum_{k=1}^{K} \\pi_{k} \\cdot N(x|\\mu_{k}, \\Sigma_{k}) $$\n\n\nEach gaussian $k$ is somehow responsible that a single data spot is placed where it is. **One gaussian may be more responsible than the others but they are all working in a mixture to explain your data**. During learning all of these ellipsoid clusters are moving in space and varying their shape trying to match the distributed data $p(x)$ perfectly. **The learning procedure for models with latent variables is the expectation maximization**. This algorithm works in two consequtive steps that are repeated until all gaussians only have found their place and would only show slight changes in further steps. \n\n#### E-Step (expectation): \n\nAt the start of the model are gaussians are placed and shaped randomly (or pretrained by k-means locations). Then for each gaussian, that represents one cluster, the model calculates how responsible the cluster $z_{k}$ is for one data spot $x_{n}$.\n\n$$\\gamma_{nk} = \\frac {\\pi_{k} \\cdot N(x_{n}|\\mu_{k}, \\Sigma_{k})} {\\sum_{j=1}^{K} \\pi_{j} \\cdot N(x_{n}|\\mu_{j}, \\Sigma_{j})}$$\n\nThese responsibilities sum up to 1 over all clusters for one data spot. Here we can read out the winner at the end of the learning procedure that is assigned as the predicted cluster for that data spot $x_{n}$.\n\n#### M-Step (maximization):\n\nDuring the maximization step, all parameters of the gaussians, the locations given by the cluster center and the shapes covered by the covariance matrices are updated:\n\n$$\\mu_{k} = \\frac{1}{N_{k}} \\sum_{n=1}^{N} \\gamma_{nk} x_{n} $$\n\n$$\\Sigma_{k} = \\frac{1}{N_{k}} \\sum_{n=1}^{N} \\gamma_{nk} (x_{n} - \\mu_{k}) (x_{n} - \\mu_{k})^{T}$$ \n\n\nYou can see that we obtain each cluster location $\\mu_{k}$ by looking and summing over all data spots. This is a big difference between K-Means and GMM as the further only calculates the cluster centers using cluster members. But in our case, as already told, each gaussian is responsibile for all data points. Hence the new location is calculated by using them all, but weighting with the responsibility. Thisway points that could be well explained by this cluster contribute much more to the new location (mean) than the other data points. The same holds for the cluster shapes aka covariances. \n\n\nBut nontheless there is also one more similarity: We are taking a weighted mean and for this reason we also need to take about outliers. We need to preprocess our data by rescaling and trying to transform the features such that they look normally distributed to suppress their effect. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Running GMM on image statistics ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Let's now run the Gaussian Mixture Model to find clusters that suite better to the elliptical shape of the groups we can observe. Instead of using a number of clusters that roughly equals the most common shapes, let's use a number with one more cluster. I have often experienced that there will be at least once cluster that holds many outliers and for this reason it could be fruitful to have one more \"outlier-cluster\" as well:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"num_clusters = num_clusters + 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gmm = GaussianMixture(n_components=num_clusters, \n max_iter=100, \n n_init=10,\n random_state=0)\n\nfeatures = [\"img_mean\", \"img_std\", \"preprocessed_img_skew\",\n \"preprocessed_red_mean\", \"preprocessed_green_mean\", \"blue_mean\"]#,\n #\"img_area\", \"rows\", \"columns\"]\n\nx = combined_stats.loc[:, features].values\ncluster_labels = gmm.fit_predict(x)\ncombined_stats[\"gmm_cluster_label\"] = cluster_labels\ncombined_stats[\"gmm_logL\"] = gmm.score_samples(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gmm.converged_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Beside the cluster label we can obtain a measure on how anomalistic an image is by computing the per sample log likelihood. This is a nice feature to find outliers in the data!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats = combined_stats.iloc[0:train_image_stats.shape[0]]\ntest_image_stats = combined_stats.iloc[train_image_stats.shape[0]::]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, now let's take a look at the cluster assigments:","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig = make_subplots(rows=1, cols=2, subplot_titles=(\"Train image stats\", \"Test image stats\"))\n\n\ntrace0 = go.Scatter(\n x = train_image_stats.img_std,\n y = train_image_stats.img_mean,\n mode='markers',\n text=train_image_stats[\"gmm_cluster_label\"],\n marker=dict(\n color=train_image_stats.gmm_cluster_label,\n colorbar=dict(thickness=10, len=1.1, title=\"cluster label\"),\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\ntrace1 = go.Scatter(\n x = test_image_stats.img_std,\n y = test_image_stats.img_mean,\n mode='markers',\n text=test_image_stats[\"gmm_cluster_label\"],\n marker=dict(\n color=test_image_stats.gmm_cluster_label,\n colorscale='Jet',\n opacity=0.4,\n size=2\n )\n)\n\nfig.add_trace(trace0, row=1, col=1)\nfig.add_trace(trace1, row=1, col=2)\n\nfig.update_xaxes(title_text=\"Image std\", row=1, col=1)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=1)\nfig.update_xaxes(title_text=\"Image std\", row=1, col=2)\nfig.update_yaxes(title_text=\"Image mean\", row=1, col=2)\n\nfig.update_layout(height=425, width=850, showlegend=False)\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* The patchwork quilt of clusters has been removed. \n* The shapes of the clusters suite better to what we can observe with our own eyes.\n* Furthermore the test data seems to have one cluster more than the training data.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Exploring clusters ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Looking at examples","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"I'm very curious how these clusters look like with GMM:","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(num_clusters,8, figsize=(20, 2.5*8))\n\nfor cluster in range(num_clusters):\n selection = np.random.choice(combined_stats[combined_stats.gmm_cluster_label==cluster].image_path.values,\n size=8, replace=False)\n m=0\n for path in selection:\n image = imread(path)\n ax[cluster, m].imshow(image)\n ax[cluster, m].set_title(\"Cluster {}\".format(cluster))\n ax[cluster, m].axis(\"off\")\n m+=1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Even though we are only using some basic statistics like the mean, std and skewness of the overall pixel values of an image as well as color channel means, we can find pretty nice clusters.\n* The result looks way better and more ordered than with K-Means.\n* We can clearly see that one cluster, that holds the big hidden image group in test that is not really present in train, differs a lot in the way the images look like compared to all other clusters. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Cluster image shapes & statistics","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In my first notebook for this competition [\"Don't turn into a Smoothie after the Shape-Up!\"](https://www.kaggle.com/allunia/don-t-turn-into-a-smoothie-after-the-shake-up) one can already see by the 3d-scatter-plot of image statistics and colored by column, texted by row that the statistics depend on the image shape. For this reason we can expect that the clusters found by GMM based on image statistics should cover specific shapes as well: ","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,5))\nsns.boxplot(combined_stats.gmm_cluster_label, combined_stats[\"rows\"], ax=ax[0], palette=\"Greens\")\nax[0].set_title(\"Rows in train and test\")\nsns.boxplot(combined_stats.gmm_cluster_label, combined_stats[\"columns\"], ax=ax[1], palette=\"Blues\")\nax[1].set_title(\"Columns in train and test\")\nsns.boxplot(combined_stats.gmm_cluster_label, combined_stats[\"img_area\"], ax=ax[2], palette=\"Reds\")\nax[2].set_title(\"Areas in train and test\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Some clusters hold very specific image shapes and this is great to see. It confirms the assumption that our different image types show different shapes. \n* Some clusters hold a broad range of rows, columns and areas. Here we have found images that look similar but show very different shapes. One way to find out wheather these are already nice groups or should be improved, is to add more features to the clustering process or to add the column and row values themselves.","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,3,figsize=(20,10))\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.img_mean, ax=ax[0,0], palette=\"Purples\")\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.img_std, ax=ax[0,1], palette=\"Oranges\")\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.preprocessed_img_skew, ax=ax[0,2], palette=\"Greys\");\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.preprocessed_red_mean, ax=ax[1,0], palette=\"Reds\")\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.preprocessed_green_mean, ax=ax[1,1], palette=\"Greens\")\nsns.violinplot(combined_stats.gmm_cluster_label, combined_stats.blue_mean, ax=ax[1,2], palette=\"Blues\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* These are the inner cluster distributions of the preprocessed features we used for clustering: the mean, the standard deviation and the skewness as well as the red, green and blue color channel means.\n* Some clusters show clearly defined scopes of the features, but there are some with broad ranges as well. \n* Furthermore some clusters hold many outlier values.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Cluster interactions ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Another great feature of mixture models is the possibility to obtain the responsibilities of each cluster per data point (in this case per image). This way we get some insight, which cluster may be an alternative for the hard assigned one. Looking at THE alternative cluster for all data points we might get an impression which image group is close to the group of interest - for example close to our hidden test group. :-)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"responsibilities = gmm.predict_proba(x)\ncombined_stats[\"alternative_cluster\"] = np.argsort(responsibilities, axis=1)[:,-2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"competition = combined_stats.groupby(\"gmm_cluster_label\").alternative_cluster.value_counts().unstack()\ncompetition.fillna(0, inplace=True)\nplt.figure(figsize=(15,15))\nsns.heatmap(competition, cmap=\"Greens\", annot=True, fmt=\"g\", cbar=False);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Using only image statistics (without further shape information) we can see that the cluster 7 with the hidden test images seems to be close to cluster 0 and 5. If you compare these clusters above you can see that only cluster 0 is really similar as we find more microscopes or light loss at image boundaries. In addition you can see that only the microscopes have shifted the test group to darker image mean values and higher standard deviations. If more groups would have microscopes we would be better able to simulate this hidden test group in our training data.\n* Furthermore we can find some clusters that have only one close cluster whereas other cluster have multiple counterparts. This could be interesting information if we would like to perform cluster-dependent augmentation. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Exploring anomalies ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"One very nice feature of GMM is the possibility to detect samples that live in regions with low density. For this purpose we need to compute the log likelihood values per sample. If you like to know more about anomaly detection with Gaussian Mixture Models, I would like to guide you to [my data science trainee notebook on GMMS](https://www.kaggle.com/allunia/hidden-treasures-in-our-groceries) that explains a bit more on \"how it works\". ;-)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(combined_stats.gmm_logL, color=\"gold\");\nplt.title(\"Sample log-likelihoods in train and test\")\nplt.ylabel(\"Density\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At the moment it's sufficient to know that low logL numbers are related to points in low dense, \"outlier\"-regions. If we would like a hard assignment or decision what is meant by an outlier, we could use a quantile measurement like this one:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"np.quantile(combined_stats.gmm_logL.values, q=0.05)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"combined_stats[\"is_outlier\"] = np.where(combined_stats.gmm_logL <= np.quantile(\n combined_stats.gmm_logL.values, q=0.05), 1, 0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats = combined_stats.iloc[0:train_image_stats.shape[0]]\ntest_image_stats = combined_stats.iloc[train_image_stats.shape[0]::]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try to find out if there is a cluster full of outliers in our data:","execution_count":null},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"outliers_in_cluster = combined_stats[combined_stats.is_outlier==1].groupby(\"gmm_cluster_label\").size() * 100\noutliers_in_cluster /= combined_stats.groupby(\"gmm_cluster_label\").size() \n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.boxplot(x=combined_stats.gmm_cluster_label, y=combined_stats.gmm_logL, ax=ax[1])\nsns.barplot(x=outliers_in_cluster.index, y=outliers_in_cluster.values, ax=ax[0])\nax[0].set_ylabel(\"% outlier in cluster\");\nax[0].set_title(\"How many outliers can be found per cluster?\")\nax[1].set_title(\"How is logL distributed within the clusters?\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Given our quantile-definition there are only 2 clusters that are occupied with > 20% outliers - cluster 5 and 6. \n* If we compare their logL-distributions with all other clusters, we can see that they show the lowest values and even their medians are lower than of all other groups. \n* Furthermore cluster 0 and 4 hold a lot of outliers as well!\n\nLet's look at some example images:","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(num_clusters,8, figsize=(20, 2.5*num_clusters))\n\nfor cluster in range(num_clusters):\n outlier_selection = np.random.choice(combined_stats[\n (combined_stats.gmm_cluster_label==cluster) & (combined_stats.is_outlier==1)\n ].image_path.values, size=8, replace=True)\n \n m=0\n for path in outlier_selection:\n image = imread(path)\n ax[cluster, m].imshow(image)\n ax[cluster, m].set_title(\"Outlier in \\n cluster {}\".format(cluster))\n ax[cluster, m].axis(\"off\")\n m+=1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Insights\n\n* Some cluster have only a few or one outliers like cluster 9, 8 and 1. Their images are often repeated only to yield this plot. \n* In contrast cluster 0, 4, 5 and 6 yield an impression of how anomalies in the data look like and for which kind of images our model may start to struggle in making good predictions.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.loc[:, \"image_name\"] = train_image_names\ntest_image_stats.loc[:, \"image_name\"] = test_image_names","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_stats.to_csv(\"train_stats_meta_cluster.csv\", index=False)\ntest_image_stats.to_csv(\"test_stats_meta_cluster.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fitting catboost & submission ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Data preparation ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_image_stats.columns.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"to_drop = [\"image_path\", \"channels\", \"image_name\", \"patient_id\"]\ntest_image_stats = test_image_stats.drop(to_drop, axis=1)\ntrain_image_stats = train_image_stats.drop(to_drop, axis=1)\ntest_image_stats = test_image_stats.drop(\"target\", axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"use_features = test_image_stats.columns.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"targets = train_image_stats.target.values\ntrain_df = train_image_stats[use_features].copy()\ntest_df = test_image_stats[use_features].copy()\ntrain_df[\"target\"] = targets","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for df in [train_df, test_df]:\n df[\"rows\"] = df[\"rows\"].astype(np.str)\n df[\"columns\"] = df[\"columns\"].astype(np.str)\n df[\"img_area\"] = df[\"img_area\"].astype(np.str)\n df[\"gmm_cluster_label\"] = df[\"gmm_cluster_label\"].astype(np.str)\n df[\"cluster_label\"] = df[\"cluster_label\"].astype(np.str)\n df[\"age_approx\"] = df[\"age_approx\"].astype(np.str)\n df[\"is_outlier\"] = df[\"is_outlier\"].astype(np.str)\n\ncat_features = np.where(test_df.dtypes==\"object\")[0]\ncat_features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.columns.values[cat_features]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df.fillna(\"NaN\")\ntest_df = test_df.fillna(\"NaN\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Validation strategy ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split, StratifiedKFold\n\nN_SPLITS = 10\nskf = StratifiedKFold(n_splits=N_SPLITS, shuffle=True, random_state=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Fitting ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class_weights = [1, train_df[train_df.target==0].shape[0] / train_df[train_df.target==1].shape[0]]\nclass_weights","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {\n 'iterations': 4000,\n 'learning_rate': 0.01,\n 'eval_metric': 'AUC',\n 'random_seed': 42,\n 'logging_level': 'Silent',\n 'use_best_model': True,\n 'loss_function': 'Logloss',\n 'od_type': 'Iter',\n 'od_wait': 1500,\n 'l2_leaf_reg': 100,\n 'depth': 5,\n 'rsm': 0.6,\n 'random_strength': 2,\n 'bagging_temperature': 10,\n 'class_weights': class_weights,\n 'random_seed': 0\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_pred = np.zeros(len(test_df))\noof = pd.DataFrame(data=np.zeros(len(train_df)), index=train_df.index.values, columns=[\"pred\"])\nfeature_importance_df = pd.DataFrame(index=use_features)\n\nm = 0\nfor train_idx, dev_idx in skf.split(train_df.drop(\"target\", axis=1).index.values, train_df.target.values):\n \n x_train, x_dev = train_df.loc[train_idx].drop(\"target\", axis=1), train_df.loc[dev_idx].drop(\"target\", axis=1)\n y_train, y_dev = train_df.loc[train_idx].target, train_df.loc[dev_idx].target\n \n train_pool = Pool(x_train, y_train, cat_features=cat_features)\n dev_pool = Pool(x_dev, y_dev, cat_features=cat_features)\n\n model = CatBoostClassifier(**params)\n model.fit(train_pool, eval_set=dev_pool, plot=True)\n \n oof.loc[dev_idx, \"pred\"] = model.predict_proba(x_dev)[:,1]\n test_pred += model.predict_proba(test_df)[:,1]/N_SPLITS\n feature_importance_df.loc[:, \"fold_\" + str(m)] = model.get_feature_importance(train_pool)\n m+=1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importance_df[\"mean\"] = feature_importance_df.mean(axis=1)\nfeature_importance_df[\"std\"] = feature_importance_df.std(axis=1)\nfeature_importance_df = feature_importance_df.sort_values(by=\"mean\", ascending=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,10))\nsns.barplot(x=feature_importance_df[\"mean\"].values,\n y=feature_importance_df.index.values, palette=\"Greens_r\");\nplt.title(\"Feature importances\");\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature importances ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"kind=None\nexplainer = shap.TreeExplainer(model)\nshap_values = explainer.shap_values(dev_pool)\nif kind==\"bar\":\n shap.summary_plot(shap_values, x_dev, plot_type=\"bar\")\nelse:\n shap.summary_plot(shap_values, x_dev)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.distplot(test_pred, ax=ax[0], kde=False, color=\"mediumseagreen\")\nax[0].set_title(\"Test predicted probability for having cancer\")\nsns.distplot(oof.pred, ax=ax[1], kde=False, color=\"red\");\nax[1].set_title(\"Oof predicted probability for having cancer\")\nfor n in range(2):\n ax[n].set_xlabel(\"predicted probability of class 1\")\n ax[n].set_ylabel(\"count\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info[\"oof\"] = oof.pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_info[\"gmm_cluster\"] = train_image_stats.gmm_cluster_label\ntrain_info[\"gmm_logL\"] = train_image_stats.gmm_logL","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"spread_of_pos_targets = train_df.groupby(\"gmm_cluster_label\").target.sum() / train_df.target.sum() * 100\nfig, ax = plt.subplots(1,2,figsize=(20,5))\n\nsns.barplot(spread_of_pos_targets.index, spread_of_pos_targets.values, ax=ax[0], palette=\"husl\");\nax[0].set_title(\"How malignant cases are distributed over clusters\");\nax[0].set_ylabel(\"% of 1-targets\");\ng = sns.countplot(train_df.gmm_cluster_label, hue=train_df.target, ax=ax[1], palette=\"Reds\")\ng.set_yscale(\"log\")\nax[1].set_title(\"Class imbalance per cluster\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.violinplot(train_info.gmm_cluster, train_info.oof, palette=\"Blues\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(basepath + \"sample_submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.target = test_pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## What you can try out!\n\nThis is just a sketch how you can find image groups and it would probably lead to some ideas how you like to preprocess the images or to perform augmentations. Furthermore you have obtained more meta_features that could be useful for working with tabular data. I would have many ideas for this competition but I won't be able to try them out. I don't like to share them all as this would spoil the fun of this competition for all having similar ideas. But there are at least 2 very obvious experiments you could try out:\n\n1. Perform a proper preprocessing for all image stats features to obtain clusters of higher quality and play with hyperparameters.\n2. Compute image statistics for external data and do some EDA and clustering with or/and without the competition data to understand why or if it's helpful to be added (not only for better class imbalance). \n\nHave a lot of fun! :-)","execution_count":null}],"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":4,"nbformat_minor":4}