{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Human Protein Atlas - Single Cell Classification\n\nThis notebook is heavily inspired and sometimes copied from [this](https://www.kaggle.com/allunia/protein-atlas-exploration-and-baseline) and [this](https://www.kaggle.com/dhananjay3/human-protein-atlas-eda-all-you-need-to-know). I wanted to find some starter code so the first mentioned kernel is a must check out if you want to find out what you should about this dataset, as the previous Human Protein Atlas was really similar. Nevertheless, I'll be sharing parts of [that](https://www.kaggle.com/allunia/protein-atlas-exploration-and-baseline) in here as I see fit. "},{"metadata":{"trusted":true},"cell_type":"code","source":"class KernelSettings:\n    \n    def __init__(self, fit_baseline=False):\n        self.fit_baseline = fit_baseline\n        \nkernelsettings = KernelSettings(fit_baseline=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Reading the data\n\nBefore diving into the cell images, I'd rather get comfortable with the data and labels first to figure out what we are dealing with."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\n\nROOT = \"../input/hpa-single-cell-image-classification/\"\ntrain =  pd.read_csv(ROOT+\"train.csv\")\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks like we have a cell ID, where each cell has multiple images but we'll get to that later, and each ID has multiple labels in the Label column. So this is a multi-label classification problem. To make it more simple, you could look at it as if it a multiple binary classification problems.\n\nWe should split these labels into separate columns and then start exploring them."},{"metadata":{"trusted":true},"cell_type":"code","source":"label_names = {\n0: \"Nucleoplasm\",\n1: \"Nuclear membrane\",\n2: \"Nucleoli\",\n3: \"Nucleoli fibrillar center\",\n4: \"Nuclear speckles\",\n5: \"Nuclear bodies\",\n6: \"Endoplasmic reticulum\",\n7: \"Golgi apparatus\",\n8: \"Intermediate filaments\",\n9: \"Actin filaments\",\n10: \"Microtubules\",\n11: \"Mitotic spindle\",\n12: \"Centrosome\",\n13: \"Plasma membrane\",\n14: \"Mitochondria\",\n15: \"Aggresome\",\n16: \"Cytosol\",\n17: \"Vesicles and punctate cytosolic patterns\",\n18: \"Negative\",\n}\n\nname_labels = dict((v,k) for k,v in label_names.items())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_labels_columns(row):\n    for label in row.Label.split('|'):\n        name = label_names[int(label)]\n        row.loc[name] = 1\n    return row\n\nfor label_name in label_names.values():\n    train[label_name] = 0\n    \ntrain = train.apply(make_labels_columns, axis=1)\ntrain_labels = train[label_names.values()]\ntrain_labels.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now each label has it's own column, so let's start exploring.\n\n## What is the frequency of each label?"},{"metadata":{"trusted":true},"cell_type":"code","source":"labels_counts = train_labels.sum(axis=0)\nlabels_counts.sort_values(ascending=True).plot(kind='barh', figsize=(10, 5), title='Frequency of Labels');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"> We can that some labels are barely present in the data like Mitotic spindle, Aggresome, etc. These insufficiency of these labels could result in certain inaccuracies during prediction as a model would not have been exposed enough to it. "},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## What the is the frequency of multiple targets? Do data points with more than 1 target dominate, or is mostly dominated by only 1 target?"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_labels.sum(axis=1).value_counts(ascending=True).plot(kind='barh', figsize=(10, 3));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"> We can see than the majority of the data is single or double targeted."},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Target correlation and co-occurence\n\n**Target correlation could be easily calculated, but co-occurence is calculated using method provided by [this kernel](https://www.kaggle.com/dhananjay3/human-protein-atlas-eda-all-you-need-to-know)**"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt \nimport seaborn as sns\n\n# Target correlation matrix\nplt.figure(figsize=(10, 7))\nsns.heatmap(train_labels[train_labels.sum(axis=1) > 1].corr(), cmap=\"icefire\", vmin=-1, vmax=1);\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"We can see how this color palette exposes that most correlation between targets are negative, except for Plasma membrane's correlation with Intermediate and Actin filaments."},{"metadata":{"trusted":true},"cell_type":"code","source":"u = train_labels\nv = u.T.dot(u)\n\nplt.figure(figsize=(10, 7))\nsns.heatmap((v / np.sum(v, axis=0)).T, cbar=True, annot=False);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"We can see that most co-occurences occur owing to the most frequent target which are Nucleosome and Cytosol, and therefore if ranked the labels according to their frequency, I guess that we would see a fading effect."},{"metadata":{"trusted":true},"cell_type":"code","source":"v = v[labels_counts.sort_values(ascending=False).index]\nv = v.reindex(labels_counts.sort_values(ascending=False).index, axis=0)\n\nplt.figure(figsize=(10, 7))\nsns.heatmap((v / np.sum(v, axis=0)).T, cbar=True, annot=False);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And we can observe the fading effect as we assumed."},{"metadata":{},"cell_type":"markdown","source":"## Helper functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.image as mpimg\n\ntargets = ['Cytosol']\nmax_rows = 5\n\n\nfor target in targets:\n    images = train.loc[np.where(train_labels[target] == 1)]\n\nfig, axes = plt.subplots(max_rows, 4, figsize=(20, 5*max_rows))\naxes = axes.flatten()\n\ncolor_filter = {'green': lambda x: ' - '.join([label_names[int(label)] for label in x.split('|')]),\n                'blue': lambda x: 'Nucleus',\n                'red': lambda x:'Microtubules',\n                'yellow': lambda x: 'Endoplasmic reticulum'}\n\nax_id = 0\nfor row, id_ in enumerate(images.ID):\n    if row == max_rows:\n        break\n        \n    for color, filter_ in color_filter.items():\n        path = f'{ROOT}train/{id_}_{color}.png'\n        img = mpimg.imread(path)\n        axes[ax_id].imshow(img)\n        axes[ax_id].set_title(filter_(images.loc[images.ID == id_, 'Label'].item()))\n        ax_id += 1\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Baseline"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install /kaggle/input/iterative-stratification/iterative-stratification-master/","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from iterstrat.ml_stratifiers import MultilabelStratifiedKFold, MultilabelStratifiedShuffleSplit\n\n\nN_SPLITS = 10\nSEED = 41295\nmskf = MultilabelStratifiedKFold(n_splits=N_SPLITS, random_state=SEED, shuffle=True)\n\npartitions = []\n\nfor train_idx, val_idx in mskf.split(train, train[label_names.values()]):\n    partition = {}\n    partition[\"train\"] = train.ID.values[train_idx]\n    partition[\"validation\"] = train.ID.values[val_idx]\n    partitions.append(partition)\n    print(\"TRAIN:\", train_idx, \"VALIDATION:\", val_idx)\n    print(\"TRAIN:\", len(train_idx), \"VALIDATION:\", len(val_idx))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class ModelParameters:\n    \"\"\"\n    Holds parameters shared between dataloader, model and image processor.\n    \"\"\"\n    \n    def __init__(self, basepath,\n                 num_classes=19,\n                 image_rows=2048,\n                 image_cols=2048,\n                 batch_size=200,\n                 n_channels=1,\n                 row_scale_factor=4,\n                 col_scale_factor=4,\n                 shuffle=False,\n                 n_epochs=1):\n        self.basepath = basepath\n        self.num_classes = num_classes\n        self.image_rows = image_rows\n        self.image_cols = image_cols\n        self.batch_size = batch_size\n        self.n_channels = n_channels\n        self.shuffle = shuffle\n        self.row_scale_factor = row_scale_factor\n        self.col_scale_factor = col_scale_factor\n        self.scaled_row_dim = np.int(self.image_rows / self.row_scale_factor)\n        self.scaled_col_dim = np.int(self.image_cols / self.col_scale_factor)\n        self.n_epochs = n_epochs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# init model parameters class\nparameters = ModelParameters(ROOT)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from skimage.transform import resize\n\nclass ImagePreprocessor:\n    \n    def __init__(self, paramters):\n        self.parameters = parameters\n        self.basepath = self.parameters.basepath\n        self.scaled_row_dim = self.parameters.scaled_row_dim\n        self.scaled_col_dim = self.parameters.scaled_col_dim\n        self.n_channels = self.parameters.n_channels\n        \n    def preprocess(self, image):\n        image = self.resize(image)\n        image = self.reshape(image)\n        image = self.normalize(image)\n        return image\n    \n    def resize(self, image):\n        image = resize(image, (self.scaled_row_dim, self.scaled_col_dim))\n        return image\n    \n    def reshape(self, image):\n        image = np.reshape(image, (image.shape[0], image.shape[1], self.n_channels))\n        return image\n    \n    def normalize(self, image):\n        image /= 255 \n        return image\n    \n    \n    def load_image(self, image_id):\n        path = f'{self.basepath}/{image_id}_green.png'\n        image = mpimg.imread(path)\n        \n        image = np.zeros(shape=(image.shape[0], image.shape[1], 4))\n        image[:,:,0] = mpimg.imread(self.basepath + image_id + \"_green\" + \".png\")\n        image[:,:,1] = mpimg.imread(self.basepath + image_id + \"_blue\" + \".png\")\n        image[:,:,2] = mpimg.imread(self.basepath + image_id + \"_red\" + \".png\")\n        image[:,:,3] = mpimg.imread(self.basepath + image_id + \"_yellow\" + \".png\")\n        return image[:, :, 0:self.n_channels]\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preprocessor = ImagePreprocessor(parameters)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# example of preprocessing an image\nid_ = train.ID[np.random.randint(0, len(train))]\ncolor = 'green'\npath = f'{ROOT}train/{id_}_{color}.png'\nimg = mpimg.imread(path)\n\npp_img = preprocessor.preprocess(img)\n\nfig, axes = plt.subplots(1, 2, figsize=(20, 10))\naxes[0].imshow(img)\naxes[0].set_title('Original Image')\naxes[1].imshow(pp_img.squeeze())\naxes[1].set_title('Preprocessed Image');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import keras\n\nclass DataGenerator(keras.utils.Sequence):\n    \n    def __init__(self, list_IDs, labels, modelparameter, imagepreprocessor):\n        self.current_epoch = 0\n        self.params = modelparameter\n        self.labels = labels\n        self.list_IDs = list_IDs\n        self.dim = (self.params.scaled_row_dim, self.params.scaled_col_dim)\n        self.batch_size = self.params.batch_size\n        self.n_channels = self.params.n_channels\n        self.num_classes = self.params.num_classes\n        self.shuffle = self.params.shuffle\n        self.preprocessor = imagepreprocessor\n        self.on_epoch_end()\n    \n    def on_epoch_end(self):\n        self.indexes = np.arange(len(self.list_IDs))\n        if self.shuffle == True:\n            np.random.shuffle(self.indexes, random_state=self.current_epoch)\n            self.current_epoch += 1\n    \n    def get_targets_per_image(self, identifier):\n        return self.labels.loc[self.labels.ID==identifier].drop(\n                [\"ID\", \"Label\"], axis=1).values\n            \n    def __data_generation(self, list_IDs_temp):\n        'Generates data containing batch_size samples' # X : (n_samples, *dim, n_channels)\n        # Initialization\n        X = np.empty((self.batch_size, *self.dim, self.n_channels))\n        y = np.empty((self.batch_size, self.num_classes), dtype=int)\n        # Generate data\n        for i, identifier in enumerate(list_IDs_temp):\n            # Store sample\n            image = self.preprocessor.load_image(identifier)\n            image = self.preprocessor.preprocess(image)\n            X[i] = image\n            # Store class\n            y[i] = self.get_targets_per_image(identifier)\n        return X, y\n    \n    def __len__(self):\n        'Denotes the number of batches per epoch'\n        return int(np.floor(len(self.list_IDs) / self.batch_size))\n    \n    def __getitem__(self, index):\n        'Generate one batch of data'\n        # Generate indexes of the batch\n        indexes = self.indexes[index*self.batch_size:(index+1)*self.batch_size]\n        # Find list of IDs\n        list_IDs_temp = [self.list_IDs[k] for k in indexes]\n        # Generate data\n        X, y = self.__data_generation(list_IDs_temp)\n        return X, y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class PredictGenerator:\n    \n    def __init__(self, predict_Ids, imagepreprocessor, predict_path):\n        self.preprocessor = imagepreprocessor\n        self.preprocessor.basepath = predict_path\n        self.identifiers = predict_Ids\n    \n    def predict(self, model):\n        y = np.empty(shape=(len(self.identifiers), self.preprocessor.parameter.num_classes))\n        for n in range(len(self.identifiers)):\n            image = self.preprocessor.load_image(self.identifiers[n])\n            image = self.preprocessor.preprocess(image)\n            image = image.reshape((1, *image.shape))\n            y[n] = model.predict(image)\n        return y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Flatten\nfrom keras.layers import Conv2D, MaxPooling2D\nfrom keras.losses import binary_crossentropy\nfrom keras.optimizers import Adadelta\nfrom keras.initializers import VarianceScaling\n\n\nclass BaseLineModel:\n    \n    def __init__(self, modelparameter):\n        self.params = modelparameter\n        self.num_classes = self.params.num_classes\n        self.img_rows = self.params.scaled_row_dim\n        self.img_cols = self.params.scaled_col_dim\n        self.n_channels = self.params.n_channels\n        self.input_shape = (self.img_rows, self.img_cols, self.n_channels)\n        self.my_metrics = ['accuracy']\n    \n    def build_model(self):\n        self.model = Sequential()\n        self.model.add(Conv2D(16, kernel_size=(3, 3), activation='relu', input_shape=self.input_shape,\n                             kernel_initializer=VarianceScaling(seed=0)))\n        self.model.add(Conv2D(32, (3, 3), activation='relu',\n                             kernel_initializer=VarianceScaling(seed=0)))\n        self.model.add(MaxPooling2D(pool_size=(2, 2)))\n        self.model.add(Dropout(0.25))\n        self.model.add(Flatten())\n        self.model.add(Dense(64, activation='relu',\n                            kernel_initializer=VarianceScaling(seed=0),))\n        self.model.add(Dropout(0.5))\n        self.model.add(Dense(self.num_classes, activation='sigmoid'))\n    \n    def compile_model(self):\n        self.model.compile(loss=keras.losses.binary_crossentropy,\n              optimizer=keras.optimizers.Adadelta(),\n              metrics=self.my_metrics)\n    \n    def set_generators(self, train_generator, validation_generator):\n        self.training_generator = train_generator\n        self.validation_generator = validation_generator\n    \n    def learn(self):\n        return self.model.fit_generator(generator=self.training_generator,\n                    validation_data=self.validation_generator,\n                    epochs=self.params.n_epochs, \n                    use_multiprocessing=True,\n                    workers=8)\n    \n    def score(self):\n        return self.model.evaluate_generator(generator=self.validation_generator,\n                                      use_multiprocessing=True, \n                                      workers=8)\n    \n    def predict(self, predict_generator):\n        y = predict_generator.predict(self.model)\n        return y\n    \n    def save(self, modeloutputpath):\n        self.model.save(modeloutputpath)\n    \n    def load(self, modelinputpath):\n        self.model = load_model(modelinputpath)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Datasets\npartition = partitions[0]\nlabels = train\n\nprint(\"Number of samples in train: {}\".format(len(partition[\"train\"])))\nprint(\"Number of samples in validation: {}\".format(len(partition[\"validation\"])))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_generator = DataGenerator(partition['train'], labels, parameters, preprocessor)\nvalidation_generator = DataGenerator(partition['validation'], labels, parameters, preprocessor)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predict_generator = PredictGenerator(partition['validation'], preprocessor, f'{ROOT}train/')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\ntest_names = submission.ID.values\n\ntest_preprocessor = ImagePreprocessor(parameters)\nsubmission_predict_generator = PredictGenerator(test_names, test_preprocessor, f'{ROOT}test/')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_labels = pd.DataFrame(data=test_names, columns=[\"ID\"])\nfor col in train_labels.columns.values:\n    if col != \"ID\":\n        test_labels[col] = 0\ntest_labels.head(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if kernelsettings.fit_baseline == True:\n    model = BaseLineModel(parameters)\n    model.build_model()\n    model.compile_model()\n    model.set_generators(training_generator, validation_generator)\n    history = model.learn()\n    \n    proba_predictions = model.predict(predict_generator)\n    baseline_proba_predictions = pd.DataFrame(index = partition['validation'],\n                                              data=proba_predictions,\n                                              columns=target_names)\n    baseline_proba_predictions.to_csv(\"baseline_predictions.csv\")\n    baseline_losses = pd.DataFrame(history.history[\"loss\"], columns=[\"train_loss\"])\n    baseline_losses[\"val_loss\"] = history.history[\"val_loss\"]\n    baseline_losses.to_csv(\"baseline_losses.csv\")\n    \n    \n    submission_proba_predictions = model.predict(submission_predict_generator)\n    baseline_labels = test_labels.copy()\n    baseline_labels.loc[:, test_labels.drop([\"ID\"].columns.values)] = submission_proba_predictions\n    baseline_labels.to_csv(\"baseline_submission_proba.csv\")\n    \n# If you already have done a baseline fit once, \n# you can load predictions as csv and further fitting is not neccessary:\nelse:\n    baseline_proba_predictions = pd.read_csv(\"../input/protein-atlas-eab-predictions/baseline_predictions.csv\", index_col=0)\n    baseline_losses = pd.read_csv(\"../input/protein-atlas-eab-predictions/baseline_losses.csv\", index_col=0)\n    baseline_labels = pd.read_csv(\"../input/protein-atlas-eab-predictions/baseline_submission_proba.csv\", index_col=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Huge thanks to [Dhananjay Raut](https://www.kaggle.com/dhananjay3) and [Laura Fink](https://www.kaggle.com/allunia) for insipiring this kernel, and more is yet to come, so stay tuned."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}