{"cells":[{"metadata":{"_uuid":"2d6da9673d8622bb6edbadb79b98b77263fc1d70"},"cell_type":"raw","source":""},{"metadata":{"_uuid":"16c5e9f6b0cc1e2be83d54adff6ec5ed16dbe0d3"},"cell_type":"markdown","source":"## Our goal\n\n* Predict various protein structures in cellular images\n* there are 28 different target proteins\n* multiple proteins can be present in one image (multilabel classification)"},{"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)\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom PIL import Image\nfrom scipy.misc import imread\n\nimport tensorflow as tf\nsns.set()\n\nimport os\nprint(os.listdir(\"../input\"))\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train_labels = pd.read_csv(\"../input/train.csv\")\ntrain_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"69a6b9c1584ebcf68b5aab128199f524daccd365"},"cell_type":"markdown","source":"How many samples do we have?"},{"metadata":{"trusted":true,"_uuid":"abaa2a7225e485102ceae9bdda12f4cf835c492f","_kg_hide-input":true},"cell_type":"code","source":"train_labels.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"882bba3d7f68f209051fea3899f046bd96dc2915"},"cell_type":"markdown","source":"## Helper code"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"beaac70d0392fab7c5b1d49649a8c994310320e4"},"cell_type":"code","source":"label_names = {\n    0:  \"Nucleoplasm\",  \n    1:  \"Nuclear membrane\",   \n    2:  \"Nucleoli\",   \n    3:  \"Nucleoli fibrillar center\",   \n    4:  \"Nuclear speckles\",\n    5:  \"Nuclear bodies\",   \n    6:  \"Endoplasmic reticulum\",   \n    7:  \"Golgi apparatus\",   \n    8:  \"Peroxisomes\",   \n    9:  \"Endosomes\",   \n    10:  \"Lysosomes\",   \n    11:  \"Intermediate filaments\",   \n    12:  \"Actin filaments\",   \n    13:  \"Focal adhesion sites\",   \n    14:  \"Microtubules\",   \n    15:  \"Microtubule ends\",   \n    16:  \"Cytokinetic bridge\",   \n    17:  \"Mitotic spindle\",   \n    18:  \"Microtubule organizing center\",   \n    19:  \"Centrosome\",   \n    20:  \"Lipid droplets\",   \n    21:  \"Plasma membrane\",   \n    22:  \"Cell junctions\",   \n    23:  \"Mitochondria\",   \n    24:  \"Aggresome\",   \n    25:  \"Cytosol\",   \n    26:  \"Cytoplasmic bodies\",   \n    27:  \"Rods & rings\"\n}\n\nreverse_train_labels = dict((v,k) for k,v in label_names.items())\n\ndef fill_targets(row):\n    row.Target = np.array(row.Target.split(\" \")).astype(np.int)\n    for num in row.Target:\n        name = label_names[int(num)]\n        row.loc[name] = 1\n    return row","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"707280dec8e4b10dc743b5473ba8d1c4492141ff"},"cell_type":"markdown","source":"## Which proteins occur most often in images?"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"13a1f4fd6b594cd4a225f6d79d966096d1c800ae"},"cell_type":"code","source":"for key in label_names.keys():\n    train_labels[label_names[key]] = 0","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"3662893c91d0f74fe666300c214f67cfb03060c7"},"cell_type":"code","source":"train_labels = train_labels.apply(fill_targets, axis=1)\ntrain_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"f516b3a1a297a167264b0f22803de231096c9719"},"cell_type":"code","source":"target_counts = train_labels.drop([\"Id\", \"Target\"],axis=1).sum(axis=0).sort_values(ascending=False)\nplt.figure(figsize=(15,15))\nsns.barplot(y=target_counts.index.values, x=target_counts.values, order=target_counts.index)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af190deeab9f05afcba339f31c7162d841e58b33"},"cell_type":"markdown","source":"### Take-Away\n\n* We can see that most common protein structures belong to coarse grained cellular components like the plasma membrane, the cytosol and the nucleus. \n* In contrast small components like the lipid droplets, peroxisomes, endosomes, lysosomes, microtubule ends, rods and rings are very seldom in our train data. For these classes the prediction will be very difficult as we have only a few examples that may not cover all variabilities and as our model probably will be confused during ins learning process by the major classes. Due to this confusion we will make less accurate predictions on the minor classes.\n* Consequently accuracy is not the right score here to measure your performance and validation strategy should be very fine. "},{"metadata":{"_uuid":"40bff26b2d170e9936c5e6ce4fbef12e7cb4f2ba"},"cell_type":"markdown","source":"## How many targets are most common?"},{"metadata":{"trusted":true,"_uuid":"313f47ac92fa1ff241148af6282c200451edbea7","_kg_hide-input":true},"cell_type":"code","source":"train_labels[\"number_of_targets\"] = train_labels.drop([\"Id\", \"Target\"],axis=1).sum(axis=1)\ncount_perc = np.round(100 * train_labels[\"number_of_targets\"].value_counts() / train_labels.shape[0], 2)\nplt.figure(figsize=(20,5))\nsns.barplot(x=count_perc.index.values, y=count_perc.values, palette=\"Reds\")\nplt.xlabel(\"Number of targets per image\")\nplt.ylabel(\"% of data\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d6e2a72c6f544f4d03a7f0e273a5819982b58108"},"cell_type":"markdown","source":"### Take-away\n\n* Most train images only have 1 or two target labels.\n* More than 3 targets are very seldom!"},{"metadata":{"_uuid":"cb957f819453857d62a96a236febabcfb9ce5a22"},"cell_type":"markdown","source":"## Which targets are correlated?\n\nLet's see if we find some correlations between our targets. This way we may already see that some proteins often come together."},{"metadata":{"trusted":true,"_uuid":"b3c0f5a77e34296416136ca2800b8623f94885d8","_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(15,15))\nsns.heatmap(train_labels[train_labels.number_of_targets>1].drop(\n    [\"Id\", \"Target\", \"number_of_targets\"],axis=1\n).corr(), cmap=\"RdYlBu\", vmin=-1, vmax=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c8083f38bf943a2d38c86544cf322ffe721d8a19"},"cell_type":"markdown","source":"### Take-away\n\n* We can see that many targets only have very slight correlations. \n* In contrast, endosomes and lysosomes often occur together and sometimes seem to be located at the endoplasmatic reticulum. \n* In addition we find that the mitotic spindle often comes together with the cytokinetic bridge. This makes sense as both are participants for cellular division. And in this process microtubules and thier ends are active and participate as well. Consequently we find a positive correlation between these targets."},{"metadata":{"_uuid":"3b7bc009eb7ec4b6d24e4e0951d5d799ac42d289"},"cell_type":"markdown","source":"## How are special and seldom targets grouped?"},{"metadata":{"_uuid":"6f0bf1de2bf7c2550cd7bb320a16ce0d79b25087"},"cell_type":"markdown","source":"### Lysosomes and endosomes\n\nLet's start with these high correlated features!"},{"metadata":{"trusted":true,"_uuid":"86c06ed4c0dbf02c605e447a5ecc40373b35771a","_kg_hide-input":true},"cell_type":"code","source":"def find_counts(special_target, labels):\n    counts = labels[labels[special_target] == 1].drop(\n        [\"Id\", \"Target\", \"number_of_targets\"],axis=1\n    ).sum(axis=0)\n    counts = counts[counts > 0]\n    counts = counts.sort_values()\n    return counts","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e2af724b721a15ebd294771629fb270c9ad9f0ae","_kg_hide-input":true},"cell_type":"code","source":"lyso_endo_counts = find_counts(\"Lysosomes\", train_labels)\n\nplt.figure(figsize=(10,3))\nsns.barplot(x=lyso_endo_counts.index.values, y=lyso_endo_counts.values, palette=\"Blues\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc6dbcea143b60d0ef360e75d14b53b09d36e7be"},"cell_type":"markdown","source":"### Rods and rings"},{"metadata":{"trusted":true,"_uuid":"5f13e4113cfb00ec8b2de6d01fb548e10b794279","_kg_hide-input":true},"cell_type":"code","source":"rod_rings_counts = find_counts(\"Rods & rings\", train_labels)\nplt.figure(figsize=(15,3))\nsns.barplot(x=rod_rings_counts.index.values, y=rod_rings_counts.values, palette=\"Greens\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fea4a1a472bf60fdc1a5317a2e616c98b781f7ab"},"cell_type":"markdown","source":"### Peroxisomes"},{"metadata":{"trusted":true,"_uuid":"11fef619f19eafc3657d693d05f7a50c786e0621","_kg_hide-input":true},"cell_type":"code","source":"peroxi_counts = find_counts(\"Peroxisomes\", train_labels)\n\nplt.figure(figsize=(15,3))\nsns.barplot(x=peroxi_counts.index.values, y=peroxi_counts.values, palette=\"Reds\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"555447f455cac28ec0c84b09b89f32b657371470"},"cell_type":"markdown","source":"### Microtubule ends"},{"metadata":{"trusted":true,"_uuid":"a543faa3994e6bd7b57fda135ea4b5bea40efed3","_kg_hide-input":true},"cell_type":"code","source":"tubeends_counts = find_counts(\"Microtubule ends\", train_labels)\n\nplt.figure(figsize=(15,3))\nsns.barplot(x=tubeends_counts.index.values, y=tubeends_counts.values, palette=\"Purples\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0a04ae9a48f60f02021c0dba920b2ac3d6d65b38"},"cell_type":"markdown","source":"### Nuclear speckles"},{"metadata":{"trusted":true,"_uuid":"85d838037fc3310b9eef31dc4e6c64e1bc35efec","_kg_hide-input":true},"cell_type":"code","source":"nuclear_speckles_counts = find_counts(\"Nuclear speckles\", train_labels)\n\nplt.figure(figsize=(15,3))\nsns.barplot(x=nuclear_speckles_counts.index.values, y=nuclear_speckles_counts.values, palette=\"Oranges\")\nplt.xticks(rotation=\"70\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ba2b973b09d2436f65e8775b214d377d46b29188"},"cell_type":"markdown","source":"### Take-away\n\n* We can see that even with very seldom targets we find some kind of grouping with other targets that reveal where the protein structure seems to be located. \n* For example, we can see that rods and rings have something to do with the nucleus whereas peroxisomes may be located in the nucleus as well as in the cytosol.\n* Perhaps this patterns might help to build a more robust model!  "},{"metadata":{"_uuid":"0714e168ea347df2111124d0f341de3172bac00d"},"cell_type":"markdown","source":"## How do the images look like?\n\n"},{"metadata":{"_uuid":"4af4b12200fd62e7726ff8a481d7566e564947ad"},"cell_type":"markdown","source":"### Peek into the directory\n\nBefore we start loading images, let's have a look into the train directory to get an impression of what we can find there:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"1d11fa62aaa798cfcb8671b7aac0e4cd4fa18776"},"cell_type":"code","source":"from os import listdir\n\nfiles = listdir(\"../input/train\")\nfor n in range(10):\n    print(files[n])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f8804380faaaeb428227dec5ee280f27c57baea"},"cell_type":"markdown","source":"Ah, ok, great! It seems that for one image id, there are different color channels present. Looking into the data description of this competition we can find that:\n\n* Each image is actually splitted into 4 different image files. \n* These 4 files correspond to 4 different filter:\n    * a **green** filter for the **target protein structure** of interest\n    * **blue** landmark filter for the **nucleus**\n    * **red** landmark filter for **microtubules**\n    * **yellow** landmark filter for the **endoplasmatic reticulum**\n* Each image is of size 512 x 512"},{"metadata":{"_uuid":"e29074f7ddba22fbdd2275fdc8782dc2dd5885c6"},"cell_type":"markdown","source":"Let's check if the number of files divided by 4 yields the number of target samples:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"8545dd9afec3e0cdc02873375ba81ff4f4658aab"},"cell_type":"code","source":"len(files) / 4 == train_labels.shape[0]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b89a5ae1b6e358d7100c79e5a2cdc5d77ff88cc2"},"cell_type":"markdown","source":"## How do images of specific targets look like?\n\nWhile looking at examples, we can build an batch loader:"},{"metadata":{"trusted":true,"_uuid":"ce954da724c8a400c12080f8c3f986a537d4a4c4"},"cell_type":"code","source":"train_path = \"../input/train/\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c7e4ffb20ebaee046d888620d0cf016f0cebf2c4","_kg_hide-input":true},"cell_type":"code","source":"def load_image(basepath, image_id):\n    images = np.zeros(shape=(4,512,512))\n    images[0,:,:] = imread(basepath + image_id + \"_green\" + \".png\")\n    images[1,:,:] = imread(basepath + image_id + \"_red\" + \".png\")\n    images[2,:,:] = imread(basepath + image_id + \"_blue\" + \".png\")\n    images[3,:,:] = imread(basepath + image_id + \"_yellow\" + \".png\")\n    return images\n\ndef make_image_row(image, subax, title):\n    subax[0].imshow(image[0], cmap=\"Greens\")\n    subax[1].imshow(image[1], cmap=\"Reds\")\n    subax[1].set_title(\"stained microtubules\")\n    subax[2].imshow(image[2], cmap=\"Blues\")\n    subax[2].set_title(\"stained nucleus\")\n    subax[3].imshow(image[3], cmap=\"Oranges\")\n    subax[3].set_title(\"stained endoplasmatic reticulum\")\n    subax[0].set_title(title)\n    return subax\n\ndef make_title(file_id):\n    file_targets = train_labels.loc[train_labels.Id==file_id, \"Target\"].values[0]\n    title = \" - \"\n    for n in file_targets:\n        title += label_names[n] + \" - \"\n    return title","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"090d1b258d2641e747e58a49376ab3d32b3893b8","_kg_hide-input":true},"cell_type":"code","source":"class TargetGroupIterator:\n    \n    def __init__(self, target_names, batch_size, basepath):\n        self.target_names = target_names\n        self.target_list = [reverse_train_labels[key] for key in target_names]\n        self.batch_shape = (batch_size, 4, 512, 512)\n        self.basepath = basepath\n    \n    def find_matching_data_entries(self):\n        train_labels[\"check_col\"] = train_labels.Target.apply(\n            lambda l: self.check_subset(l)\n        )\n        self.images_identifier = train_labels[train_labels.check_col==1].Id.values\n        train_labels.drop(\"check_col\", axis=1, inplace=True)\n    \n    def check_subset(self, targets):\n        return np.where(set(self.target_list).issuperset(set(targets)), 1, 0)\n    \n    def get_loader(self):\n        filenames = []\n        idx = 0\n        images = np.zeros(self.batch_shape)\n        for image_id in self.images_identifier:\n            images[idx,:,:,:] = load_image(self.basepath, image_id)\n            filenames.append(image_id)\n            idx += 1\n            if idx == self.batch_shape[0]:\n                yield filenames, images\n                filenames = []\n                images = np.zeros(self.batch_shape)\n                idx = 0\n        if idx > 0:\n            yield filenames, images\n            ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"61ba9ebe5b6bfbea57d24e43bc8f4a1c2638d5eb"},"cell_type":"markdown","source":"Let's try to visualize specific target groups. **In this example we will see images that contain the protein structures lysosomes or endosomes**. Set target values of your choice and the target group iterator will collect all images that are subset of your choice:"},{"metadata":{"trusted":true,"_uuid":"445f454dc36c1f4056bc870e8810f219f1c560ea"},"cell_type":"code","source":"your_choice = [\"Lysosomes\", \"Endosomes\"]\nyour_batch_size = 3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6dcd0c8f57651ba21b82420edd37f70f93953215","_kg_hide-input":true},"cell_type":"code","source":"imageloader = TargetGroupIterator(your_choice, your_batch_size, train_path)\nimageloader.find_matching_data_entries()\niterator = imageloader.get_loader()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7d1f440b6ed70150f46a76970b544752f0d194a5"},"cell_type":"markdown","source":"To keep the kernel dense, the target group iterator has a batch size which stands for the number of examples you like to look at once. In this example you can see a maximum amount of 3 images at one iteration.  **To observe the next 3 examples of your target group, just run the cell below again.** This way you can run the cell until you have seen all images of your group without polluting the kernel:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"5f7c98a6e15c081fcbe54c4e678ea2ea5ae4fbd7"},"cell_type":"code","source":"file_ids, images = next(iterator)\n\nfig, ax = plt.subplots(len(file_ids),4,figsize=(20,5*len(file_ids)))\nif ax.shape == (4,):\n    ax = ax.reshape(1,-1)\nfor n in range(len(file_ids)):\n    make_image_row(images[n], ax[n], make_title(file_ids[n]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a66576fbdab4678c42d183017335d1ef27eed655"},"cell_type":"markdown","source":"### Take-Away\n\n* Looking at this few examples we can already obtain some insights:\n    * The staining of target proteins in the green channel was not equally successful. The **images differ in their intensities and the target proteins are not always located the same way**. The first image you can get by the loader shows endosomes that are spread all over the cells and in the second and third you can find endosomes and lysosomes more concetrated around the nucleus. \n    * Especially **in the red channel we can see morphological differences**. It looks like if the cells are of different types. This is just an assumption but perhaps one could use the red channel information to reveal cell types. "},{"metadata":{"_uuid":"ad5c5a22a659540bcb3da81453a6ef2cfaa958cc"},"cell_type":"markdown","source":"## Building a baseline model"},{"metadata":{"_uuid":"8dd65575d8a244f142b060aeb5cca30db9f27576"},"cell_type":"markdown","source":"### K-Fold Cross-Validation"},{"metadata":{"_uuid":"63437e877ed24bd86d34d5c7e8cd9f9e98d29d59"},"cell_type":"markdown","source":"Let's see how many test and train samples we have in this competition:"},{"metadata":{"trusted":true,"_uuid":"b2d74fa662ec90ebc2912105848daf7301845ef8","_kg_hide-input":true},"cell_type":"code","source":"train_files = listdir(\"../input/train\")\ntest_files = listdir(\"../input/test\")\npercentage = np.round(len(test_files) / len(train_files) * 100)\n\nprint(\"The test set size turns out to be {} % compared to the train set.\".format(percentage))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2d9613edd4f1a05b54d79dbe142b78e8cbec7625"},"cell_type":"markdown","source":"To understand the performance of our model we will use **k-fold cross validation**. The train data is splitted into k chunks and each chunk is used once for testing the prediction performance whereas the others are used for training. As our targets show relationships seemed to be grouped somehow the performance per test chunk probably highly depends on the target distribution per test chunk. For example there could be chunks with very seldom targets that may obtain a bad score and some chunks with very common targets and a very good score. To reduce this effect, we will **repeat the K-Fold several times** and look at scoing distributions in the end."},{"metadata":{"_uuid":"d94c71d94925daa6fe0a71a07c1da985925f7d91"},"cell_type":"markdown","source":"As our test data is 38 % of size compared to the train set it makes sense to use 3-Fold cross validation where the test set is 33 % of size compared to the train set. As we are working with neural networks that can be demanding in computational resources, let's only use 2 repetitions. "},{"metadata":{"trusted":true,"_uuid":"29a024cbd71145018c36342a5e8a483b6e010f5b"},"cell_type":"code","source":"from sklearn.model_selection import RepeatedKFold\n\nsplitter = RepeatedKFold(n_splits=3, n_repeats=1, random_state=0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"37fc0f4495c93fc45379aab2eb156fdbd4d3f7c3"},"cell_type":"markdown","source":"This splitter is now a generator. Hence if you call splitters split method it will yield one Fold of the repeated K-Folds. Consequently if we choose n_repeats=2 we will end up with 6 Folds in total: 3 Folds for the first cross validation and again 3 Folds for the repeated cross validation. We will perform the splitting on the image ids. This way we can easily load images and targets given the chunk ids. **Due to performance reasons I will only use one cv-fold to explore results and one repeat!**"},{"metadata":{"trusted":true,"_uuid":"539a35dfd7d735d7068f32277146e60e15729e55"},"cell_type":"code","source":"partitions = []\n\nfor train_idx, test_idx in splitter.split(train_labels.index.values):\n    partition = {}\n    partition[\"train\"] = train_labels.Id.values[train_idx]\n    partition[\"validation\"] = train_labels.Id.values[test_idx]\n    partitions.append(partition)\n    print(\"TRAIN:\", train_idx, \"TEST:\", test_idx)\n    print(\"TRAIN:\", len(train_idx), \"TEST:\", len(test_idx))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"65f4c326d8f32f13872d0dfaaad2bbbfdda257a7"},"cell_type":"code","source":"partitions[0][\"train\"][0:5]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e0c336986b46a35c224dd447190fbb8c553bb64a"},"cell_type":"markdown","source":"### Collecting ideas\n\nNext we need to setup a simple baseline model. This need not be very complex or very good. Its our first attempt to play with and to figure out how to improve. For this purpose let's use the deep learning library [keras](https://keras.io/). This tools makes it easy for us to build and train neural networks. First of all, we should collect some ideas:\n\n* To **stay simple let's use only the green channel image of our images per id**. The competition says that it shows the stained target proteins and consequently it's hopefully the most informative one. The other images are like references showing microtubules, nucleus and endoplasmatic reticulum. We don't acutally now how informative they are and in our current state they would blow up our neural network with a huge amount of network weigths that we might not need.\n* Let's use **generators to only load data images of our batch and not all in once**. Using keras fit_generator, evaluate_generator and predict_generator we can directly connect them to keras without worrying much about how keras does its job. For this purpose I highly follow a descprition of a post in the www for which you will find the link below.\n* It could be advantegous to write a **small class that does simple preprocessing per image.** This way we can easily change something of this phase without producing chaos in the model itself or during data loading.   \n* I'm going to use a **small class that hold parameters that are used or shared between the data loader, the image preprocessor and the baseline model**. Passing an instance of this class to them reduced the risk of setting different parameters and obtaining mismatch errors for example during build & compile of the network layers. \n"},{"metadata":{"_uuid":"5cd298c9360c772a9c8b7946425449db8be0ab8f"},"cell_type":"markdown","source":"### Shared Parameter class"},{"metadata":{"trusted":true,"_uuid":"96df4f1aa1085b95b806b48071c0737d45a5b708","_kg_hide-input":true},"cell_type":"code","source":"class ModelParameter:\n    \n    def __init__(self, num_classes=28,\n                 image_rows=512,\n                 image_cols=512,\n                 batch_size=128,\n                 n_channels=1,\n                 row_scale_factor=4,\n                 col_scale_factor=4,\n                 shuffle=True,\n                 n_epochs=1):\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 = True\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":{"_uuid":"bc20cc4fdebf53b66bda3899974bd8758159d85f"},"cell_type":"markdown","source":"Ok, now we will create an instance of this class and pass it to the DataGenerator, the BaseLineModel and the ImagePreprocessor."},{"metadata":{"trusted":true,"_uuid":"96b80b33c817ad07deac345e36c46c062255c60d"},"cell_type":"code","source":"parameter = ModelParameter()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4497282fcdc0d68554e3878bb82bae8942eda681"},"cell_type":"markdown","source":"### Image Preprocessor\n\nLet's write a simple image preprocessor that handles for example the rescaling of the images. Perhaps we can expand its functionality during improvement of the baseline model. "},{"metadata":{"trusted":true,"_uuid":"2674c1331715aeab66a0194b3ad31afa7cfc3444","_kg_hide-input":true},"cell_type":"code","source":"from skimage.transform import resize\n\nclass ImagePreprocessor:\n    \n    def __init__(self, modelparameter):\n        self.parameter = modelparameter\n        self.scaled_row_dim = self.parameter.scaled_row_dim\n        self.scaled_col_dim = self.parameter.scaled_col_dim\n        self.n_channels = self.parameter.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        ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d93eeb59fa1f35d3e46f24832a40bee0086abb5b"},"cell_type":"markdown","source":"Let's create an instance of this preprocessor and pass it to the data generator."},{"metadata":{"trusted":true,"_uuid":"589121752a0ca1693c7e794655047f947a476b8e"},"cell_type":"code","source":"preprocessor = ImagePreprocessor(parameter)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c83a1888064cf71767aee2d462171227836600f5"},"cell_type":"markdown","source":"#### Looking at a preprocessed example image"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"cc9f0ff3eecc870337c475533b4c1cffed704403"},"cell_type":"code","source":"example = images[0,0]\npreprocessed = preprocessor.preprocess(example)\nprint(example.shape)\nprint(preprocessed.shape)\n\nfig, ax = plt.subplots(1,2,figsize=(20,10))\nax[0].imshow(example, cmap=\"Greens\")\nax[1].imshow(preprocessed.reshape(parameter.scaled_row_dim,parameter.scaled_col_dim), cmap=\"Greens\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"615d734ed90fb1def88151c5365443d61b411536"},"cell_type":"markdown","source":"You can see that we have lost a lot of information by downscaling the image!"},{"metadata":{"_uuid":"2b3cecd1dde678563c35405a7ddde5c35314677b"},"cell_type":"markdown","source":"### Data Generator\n\nI highly build upon the [nice data generator presented by Shervine Amidi.](https://stanford.edu/~shervine/blog/keras-how-to-generate-data-on-the-fly) Thank you! :-) \n\n\n** WORK in progess! - Currently there are missing samples of the last batch < batch_size** "},{"metadata":{"trusted":true,"_uuid":"6e77e06f77164446a8689b55fd56734f943368b1","_kg_hide-input":true},"cell_type":"code","source":"import keras\n\nclass DataGenerator(keras.utils.Sequence):\n    \n    def __init__(self, basepath, list_IDs, labels, modelparameter, imagepreprocessor):\n        self.basepath = basepath\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)\n    \n    def get_targets_per_image(self, identifier):\n        return self.labels.loc[self.labels.Id==identifier].drop(\n                [\"Id\", \"Target\", \"number_of_targets\"], 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 = load_image(self.basepath, identifier)[0]\n            image = self.preprocessor.preprocess(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.ceil(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":{"_uuid":"48e03c679baabe6f4a57cb9da335ccb6d6f9710b"},"cell_type":"markdown","source":"### CNN Baseline model using keras"},{"metadata":{"trusted":true,"_uuid":"76a48f7a606ef5293bfc465871c707ce9b3900a8","_kg_hide-input":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\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        self.model.add(Conv2D(32, (3, 3), activation='relu'))\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        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):\n        return self.model.predict_generator(generator=self.validation_generator,\n                                           use_multiprocessing=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"79ae4ddf40420cb60ed564e515dc6a8fe3fd2d82"},"cell_type":"markdown","source":"### Training the baseline on the first cv-fold"},{"metadata":{"trusted":true,"_uuid":"6bdc3802ae488499de2afacb32e2e93bee5dccf4"},"cell_type":"code","source":"# Datasets\npartition = partitions[0]\nlabels = train_labels\n\n# Generators\ntraining_generator = DataGenerator(train_path, partition['train'], labels, parameter, preprocessor)\nvalidation_generator = DataGenerator(train_path, partition['validation'], labels, parameter, preprocessor)\n\n\nmodel = BaseLineModel(parameter)\nmodel.build_model()\nmodel.compile_model()\n\nmodel.set_generators(training_generator, validation_generator)\nhistory = model.learn()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"98c5206497d690c27fca5fb0cc2866ab31b3f673"},"cell_type":"markdown","source":"Even though this accuracy looks nice it's an illusion! We are far away from a good model. Let's try to understand why..."},{"metadata":{"_uuid":"5e093fa4cd74f6ceda7f395c2017a39a335f2fe9"},"cell_type":"markdown","source":"## What do the results tell us?\n\nLet's have a look at predicted probabilites per target class:"},{"metadata":{"trusted":true,"_uuid":"8ed1e7784c82006e8aa08eb1e18a06cd4af8ecd8"},"cell_type":"code","source":"proba_predictions = model.predict()\nvalidation_labels = train_labels.loc[train_labels.Id.isin(partition[\"validation\"])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22fb79c1ff6d84f7403ad26418b66669e73fab6a"},"cell_type":"code","source":"print(validation_labels.shape)\nprint(proba_predictions.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1bdefea146c0c50bc9d78fb2ddf41c2f9e4b1f14"},"cell_type":"code","source":"proba_predictions = proba_predictions[0:-10]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"851060fb988a7942a3325349e9999f57f2a6225f"},"cell_type":"code","source":"hot_values = validation_labels.drop([\"Id\", \"Target\", \"number_of_targets\"], axis=1).values.flatten()\none_hot = (hot_values.sum()) / hot_values.shape[0] * 100\nzero_hot = (hot_values.shape[0] - hot_values.sum()) / hot_values.shape[0] * 100\n\nfig, ax = plt.subplots(1,2, figsize=(20,5))\nsns.distplot(proba_predictions.flatten() * 100, color=\"DodgerBlue\", ax=ax[0])\nax[0].set_xlabel(\"Probability in %\")\nax[0].set_ylabel(\"Density\")\nax[0].set_title(\"Predicted probabilities\")\nsns.barplot(x=[\"label = 0\", \"label = 1\"], y=[zero_hot, one_hot], ax=ax[1])\nax[1].set_ylim([0,100])\nax[1].set_title(\"True target label count\")\nax[1].set_ylabel(\"Percentage\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4ee24c3b3799e06b86d54c5a24323b6b5902d0ed"},"cell_type":"markdown","source":"### Take-Away\n\n* We can see that our model was always very uncertain to predict the presence of a target protein. All probabilities are close to zero and there are only a few with targets where our model predicted a protein structure with higher than 10 %.\n* If we take a look at the true target label count we can see that most of our targets are filled with zero. This corresponds to an absence of corresponding target proteins. This makes sense: For each image we have a high probability to contain either 1 or 2 target protein structures. Their label values are one whereas all others are zero. \n* Consequently our high accuracy belongs to the high correct prediction of the absence of target proteins. In contrast we weren't able to predict the presence of a target protein which is the most relevant part! \n* Now a bell should ring :-) Have you ever heard about imbalanced classes and model confusion? "},{"metadata":{"_uuid":"6475127409a86d7c4da3544478456d0d17c148eb"},"cell_type":"markdown","source":"### To which targets do the high and small predicted probabilities belong to?"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7174d05d4e824c290df5149e6990827a076755a7"},"cell_type":"code","source":"mean_predictions = np.mean(proba_predictions, axis=0)\nstd_predictions = np.std(proba_predictions, axis=0)\nmean_targets = validation_labels.drop([\"Id\", \"Target\", \"number_of_targets\"], axis=1).mean()\n\nlabels = validation_labels.drop([\"Id\", \"Target\", \"number_of_targets\"], axis=1).columns.values\n\nfig, ax = plt.subplots(1,2,figsize=(20,5))\nsns.barplot(x=labels,\n            y=mean_predictions,\n            ax=ax[0])\nax[0].set_xticklabels(labels=labels,\n                      rotation=90)\nax[0].set_ylabel(\"Mean predicted probability\")\nax[0].set_title(\"Mean predicted probability per class over all samples\")\nsns.barplot(x=labels,\n           y=std_predictions,\n           ax=ax[1])\nax[1].set_xticklabels(labels=labels,\n                      rotation=90)\nax[1].set_ylabel(\"Standard deviation\")\nax[1].set_title(\"Standard deviation of predicted probability per class over all samples\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b058c1f00898dad04957433672e726822426e478"},"cell_type":"code","source":"fig, ax = plt.subplots(1,1,figsize=(20,5))\nsns.barplot(x=labels, y=mean_targets.values, ax=ax)\nax.set_xticklabels(labels=labels,\n                      rotation=90)\nax.set_ylabel(\"Percentage of hot (1)\")\nax.set_title(\"Percentage of hot counts (ones) per target class\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e3717835b45904306c0b00f66949ef8f8f7f86a2"},"cell_type":"markdown","source":"### Take-Away\n\n* Our baseline model seemed to learn something even if this something does not look very nice. Instead of predicting 0.5 probability of each target class if yields the fraction of the number of one-hot-targets devided by the number of of zero-hot targets. Taking a look at the standard deviation we can see that all samples have the same predicted values. There is no deviation, no difference between them. This is of course very bad!\n\nHence, our next goal should be to find out what to tune in such a way that our model really starts learning! "},{"metadata":{"_uuid":"e9c23453f175fc3768f4b94b659c96f8fb72c547"},"cell_type":"markdown","source":"## One-Step-Improvement\n\n\nOk, again let's go one step back and choose the most common target proteins that are present in our data: nucleoplasm, cytosol and plasma membrane. If we are not able to predict them we can go home and stay in bed ;-) .\n\n### The target wish list\n\nTo introduce a target wishlist that we can change whenever we want we need to improve the data generator. For this purpose we're going to extend the class we have already written. Taking a closer look at the base generator you can see that there is just one line code in def data_generation(self, list_IDs_temp) we have to change, namely the part with y[i] = ... inside the for loop over temp list ids (image identifiers of the batch). To make things easier, I added a small method to the DataGenerator we already had:\n\n```\ndef get_targets_per_image(self, identifier):\n        return self.labels.loc[self.labels.Id==identifier].drop(\n                [\"Id\", \"Target\", \"number_of_targets\"], axis=1).values\n```\n\nThis method just avoids the direct pass to y[i], the targets per image in a batch. Now, we can overwrite this method in our ImprovedDataGenerator without loosing functionality:"},{"metadata":{"trusted":true,"_uuid":"24b1f61b38bc5ec85150f0ffbeb1000d2b4c4ee2"},"cell_type":"code","source":"class ImprovedDataGenerator(DataGenerator):\n    \n    # in contrast to the base DataGenerator we add a target wishlist to init\n    def __init__(self,basepath, list_IDs, labels, modelparameter, imagepreprocessor, target_wishlist):\n        super(DataGenerator, self).__init__(basepath, list_IDs, labels, modelparameter, imagepreprocessor)\n        self.target_wishlist = target_wishlist\n    \n    def get_targets_per_image(self, identifier):\n        return self.labels.loc[self.labels.Id==identifier, self.target_wishlist].values","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"66ad84e730aea51fed90d03a40e6615726ea88fb"},"cell_type":"markdown","source":"### Plug and Play with the baseline model\n\n:-) \n\nThis part makes fun! Actually I don't know if our model learns something. But we can try to find out by improving our model adding **new features** and playing with **different parameter settings**. The latter I like to do in a **plug-and-play style**: Setting the flag *improve=True* adds a change to our model whereas *improve=False* uses the old concept we already used in the baseline.  \n\n#### Add scoring metrics\n\nWe have already seen that the accuracy score is an illusion and does not help to figure out how good our predictions are. Let's take a closer look to the competition scoring and alternatives:\n\n* **F1 macro score**: Check out this [nice implementation of Guglielmo Camporese](https://www.kaggle.com/guglielmocamporese/macro-f1-score-keras). Thank you very much! We can easily add it to our model.\n* But even with that score we should be careful! We have 28 different classes that are **very different in their frequency of being present**.  In addition we have to deal with **highly imbalanced classes per single target**. Even for the most common target nucleoplasm there are only 40 % of samples that show it and 60 % not. This imbalance becomes even more dramatic for seldom targets like rods and rings. We should **attach more importance to true positives**. "},{"metadata":{"trusted":true,"_uuid":"b1d64d6e4dbc439fb135cd77d0d0ef2d65706c96","_kg_hide-input":true},"cell_type":"code","source":"import keras.backend as K\n\ndef f1(y_true, y_pred):\n    y_pred = K.round(y_pred)\n    tp = K.sum(K.cast(y_true*y_pred, 'float'), axis=0)\n    tn = K.sum(K.cast((1-y_true)*(1-y_pred), 'float'), axis=0)\n    fp = K.sum(K.cast((1-y_true)*y_pred, 'float'), axis=0)\n    fn = K.sum(K.cast(y_true*(1-y_pred), 'float'), axis=0)\n\n    p = tp / (tp + fp + K.epsilon())\n    r = tp / (tp + fn + K.epsilon())\n\n    f1 = 2*p*r / (p+r+K.epsilon())\n    f1 = tf.where(tf.is_nan(f1), tf.zeros_like(f1), f1)\n    return K.mean(f1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0b763fe315bbc0960599bb4149bce9f600ef45f3"},"cell_type":"code","source":"class ImprovedModel(BaseLineModel):\n    \n    def __init__(self, modelparameter, my_metrics=[f1]):\n        super(BaseLineModel, self).init(modelparameter)\n        self.my_metrics = my_metrics","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9a0fbfd1d2ab88ea2be51687311a5da8b6be7bc1"},"cell_type":"markdown","source":"## Where to go next?"},{"metadata":{"trusted":true,"_uuid":"9d0b2e6b4af7dff12d40abd4711c95ecca43c94f"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b28cf97bcf95f7e77b95589aa8d751a8fcacabea"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}