{"cells":[{"metadata":{"_uuid":"9d90940457d5219fcca238343e552f6ff2fc88dd"},"cell_type":"markdown","source":"## Import "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import tensorflow\nimport sys\nimport numpy as np\nimport keras\nfrom keras.utils import Sequence\nfrom PIL import Image\nfrom matplotlib import pyplot as plt\nimport pandas as pd\nfrom tqdm import tqdm\nimport os\nimport imgaug as ia\nfrom imgaug import augmenters as iaa\nimport cv2\nfrom matplotlib import image\n\nBASE_DIR = \"../input\"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ac637dc1eeb6582df9cfdde234c7f997fdd2db01"},"cell_type":"markdown","source":"## Data Exploration"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"def train_data():\n    ''' Gets paths to images for training data and labels '''\n    x_train_path  = os.path.join(BASE_DIR,'train')\n    y_train_path  = os.path.join(BASE_DIR,'train.csv')\n    data          = pd.read_csv(y_train_path)\n    paths         = []\n    labels        = []\n    for example_id, protein_ids in zip(data['Id'], data['Target'].str.split(' ')):\n        y = np.zeros(28)\n        y[np.array(protein_ids,dtype=np.int32)]=1\n        paths.append(os.path.join(x_train_path, example_id))\n        labels.append(y)\n    return np.array(paths), np.array(labels)\n\ndef test_data():\n    ''' Gets paths to images for test data '''\n    x_test_path = os.path.join(BASE_DIR,'test')\n    data_test   = pd.read_csv(os.path.join(BASE_DIR,'sample_submission.csv'))\n    paths       = []\n    y           = []\n    for example_id in data_test.Id:\n        paths.append(os.path.join(x_test_path,example_id))\n        y.append(np.zeros(28))\n    return np.array(paths), np.array(y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"83dea6f760ffd8c7b73b3467f028353004585133"},"cell_type":"code","source":"paths, labels = train_data()\n\nlabel_counts = np.sum(labels,axis=0)\nprotein_id   = np.arange(len(label_counts))\n\nlocalization_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\n# proteins 8,9,10 and 28 are severely underrepresented in the \nplt.figure(1, figsize=(8,5))\nplt.bar(protein_id,label_counts/1000,color='c')\nplt.xticks(protein_id,[localization_names[int(el)] for el in protein_id],rotation='vertical')\nplt.xlabel('Localizations')\nplt.ylabel('Frequency in K')\nplt.title('Frequency of protein localizations')\nplt.box(on=None) \nplt.show()\n\n\nsorted_by_counts = sorted(zip(protein_id,label_counts),key = lambda item: item[1], reverse=False)\nfrequency_counts = pd.DataFrame(data=np.array(sorted_by_counts,dtype=np.int32), columns=['ProteinID','Counts'])\nprint(frequency_counts)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e14226a4d8db37f5d4abdd6a49fa41bd8bad2a76"},"cell_type":"code","source":"# How many localizations there are per example\nlabels_per_image = np.sum(labels,axis=1)\nfrom collections import Counter\ncnt = sorted(Counter(labels_per_image).items(),key=lambda item: item[0])\nlpi    = [int(e[0]) for e in cnt]\ncounts = [int(e[1]) for e in cnt] \n\n\nplt.figure()\nplt.bar(lpi,counts,color='g')\nplt.xticks(lpi, [str(e) for e in lpi])\nplt.xlabel('Labels per image')\nplt.ylabel('Count')\nplt.box(on=None) \nplt.title('Number of labels per data point')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"17dec5b04d46f3645d01dfbb875278b44c2adeb7"},"cell_type":"code","source":"# We want to make split so that train/val/test sets will look the same way, due to \n# significant data imbalances random split may result lead to rarely occuring classes \n# being absent in train or test or validation. In particular we are interested in making \n# sure that examples that contain protein with ID=27 are included.\n\nfrom sklearn.model_selection import train_test_split\n\npaths,labels = train_data()\npaths_train, paths_val, labels_train, labels_val = train_test_split(paths,labels,test_size=0.2,stratify=labels[:,27])\n#paths_val, paths_test, labels_val, labels_test = train_test_split(paths_val_test,labels_val_test,test_size=0.5,stratify=labels_val_test[:,27])\n\n# As you can see from the plot train/val/test have the same distribution\nplt.figure(figsize=(12,6))\nfor series,name,clr in zip([labels_train,labels_val],['train','val'],['ro','bo']):\n     plt.plot(np.arange(series.shape[1]),series.sum(axis=0)/series.sum(),clr)\nplt.xlabel('Protein ID')\nplt.ylabel('Proportion of label')\nplt.title('Comparison of train/val/test')\nplt.show()\n\n# train/val/test contain Protein with ID = 27\nprint(labels_train[:,27].sum())\nprint(labels_val[:,27].sum())\n#print(labels_test[:,27].sum())\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0867b7082477a6b05528e04b5a1d3d45724f15ff"},"cell_type":"code","source":"# Vizualize second image \npath = paths_train[1]\npath_completions = [\"_green.png\",'_yellow.png','_blue.png','_red.png']\ncolor_maps       = ['Greens','Oranges','Blues','Reds']\nplt.figure(figsize=(6,6))\nplt.subplots_adjust(left=0.05, right=0.95)\nfor i, (cmpl,cmp) in enumerate(zip(path_completions,color_maps)):\n    A = image.imread(path+cmpl)\n    plt.subplot(2,2,i+1)\n    plt.imshow(A,cmap=cmp)\n    plt.title(cmp)\n    plt.box(on=None) \nplt.show()   ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"242837adea942f71858d1749b36036b4acc323be"},"cell_type":"code","source":"# loading single image\n\nimage_yellow = image.imread(paths_train[0]+'_yellow.png')\nimage_green  = image.imread(paths_train[0]+'_green.png')\nimage_blue   = image.imread(paths_train[0]+'_blue.png')\nimage_red    = image.imread(paths_train[0]+'_red.png')\n\nimage_general = np.concatenate([np.expand_dims(channel,axis=-1) for channel in (image_yellow,image_green,image_blue,image_red)],axis=-1)\n\ndef image_loader(image_path, width=512, height=512, channels_to_include=['green']):\n    ''' Image loader '''\n    image_matrix = np.zeros([width,height,len(channels_to_include)],dtype=np.float32)\n    channel_index = 0\n    for i,channel_color in enumerate(channels_to_include):\n        channel_extension = \"_\"+channel_color+'.png'\n        image_matrix[:,:,i] = image.imread(image_path+channel_extension)\n    return image_matrix","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"11857fbcc76ea58ccf90127dd4210b01b70e5bd5"},"cell_type":"markdown","source":"## Input Pipeline"},{"metadata":{"trusted":true,"_uuid":"6612a2e21d8fc4a8d87f18f56c49e1ebd1b4a2ae"},"cell_type":"code","source":"from skimage.transform import resize\nfrom itertools import islice\n\ndef train_generator_single_epoch(paths,labels,batch_size=128,row_resize_ratio=4,col_resize_ratio=4,\n                                width=512,height=512,channels_to_include=['green','blue','red','yellow']):\n    ''' Generator for training data that will automatically resize image to a given proportions '''\n    iter_idx = 0\n    while iter_idx*batch_size < labels.shape[0]:\n        batch_paths  = paths[iter_idx*batch_size:(iter_idx+1)*batch_size]\n        batch_images = [image_loader(path,width,height,channels_to_include) for path in batch_paths]\n        batch_labels = labels[iter_idx*batch_size:(iter_idx+1)*batch_size,:]\n        if row_resize_ratio==1 and col_resize_ratio==1:\n            iter_idx += 1\n            yield batch_images,batch_labels\n        else:\n            iter_idx += 1\n            resize_shape = (int(width/row_resize_ratio),int(height/col_resize_ratio),len(channels_to_include))\n            batched_resized_images = [resize(img,resize_shape) for img in batch_images]\n            yield np.array(batched_resized_images),batch_labels\n\n# vizual comparison of resized and original image\ntrain_iterator_original   = train_generator_single_epoch(paths_train,labels_train,batch_size=1,row_resize_ratio=1,col_resize_ratio=1)\ntrain_iterator_resized    = train_generator_single_epoch(paths_train,labels_train,batch_size=1,row_resize_ratio=4,col_resize_ratio=4)\nimage_original,_          = next(train_iterator_original)\nimage_resized,_           = next(train_iterator_resized)\nplt.figure(figsize=(6,6))\nplt.subplots_adjust(left=0.05, right=0.95)\nplt.subplot(1,2,1)\nplt.imshow(np.squeeze(image_resized),cmap='Greens')\nplt.title('Resized Image')\n\n\nplt.subplot(1,2,2)\nplt.imshow(np.squeeze(image_original),cmap='Greens')\nplt.title('Original Image')\nplt.show()\n\n# Comments: it is quite obvious from these images that resizing to the size (128,128) from (512,512) resulted in\n# more blurry image and probably significant loss of information.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c4c5c5efe74897ba9cbef211f29c08bf3681c65c"},"cell_type":"code","source":"def eval_predict_proba_generator(paths,predict_fn,batch_size=128,row_resize_ratio=4,col_resize_ratio=4,\n                                 width=512,height=512,channels_to_include=['green']):\n    iter_idx    = 0\n    accumulator = []\n    resize_shape = (int(width/row_resize_ratio),int(height/col_resize_ratio),len(channels_to_include))\n    while iter_idx*batch_size < paths.shape[0]:\n        batch_paths  = paths[iter_idx*batch_size:(iter_idx+1)*batch_size]\n        batch_images = [image_loader(path,width,height,channels_to_include) for path in batch_paths]\n        if row_resize_ratio==1 and col_resize_ratio==1:\n            iter_idx += 1\n            predicted_probs = predict_fn(np.concatenate([np.expand_dims(im,axis=0) for im in batch_images],axis=0))\n            accumulator.append(np.concatenate(predicted_probs,axis=1))\n        else:\n            iter_idx += 1\n            batched_resized_images = [resize(img,resize_shape) for img in batch_images]\n            predicted_probs = predict_fn(np.concatenate([np.expand_dims(im,axis=0) for im in batched_resized_images],axis=0))\n            accumulator.append(np.concatenate(predicted_probs,axis=1))\n    return np.concatenate(accumulator,axis=0)\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"934962ae9f143589a72b13d5710486a51a43728f"},"cell_type":"markdown","source":"## Build Model"},{"metadata":{"trusted":true,"_uuid":"db864ac73cb610877c9a8e2cdc40a82648bbd974"},"cell_type":"code","source":"\ndef compute_task_specific_gamma(y_true):\n    '''Heuristic to compute values of gamma for different tasks'''\n    props = np.sum(y_true,0) / y_true.shape[0]\n    gammas = [0]*28\n    for i in range(y_true.shape[1]):\n        if props[i] > 0.2:\n            gammas[i] = 2\n        elif props[i] >= 1e-2 and props[i] <= 0.2:\n            gammas[i] = 4\n        else:\n            gammas[i] = 6\n    return gammas\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"502b82bb21a9eccefa6a5e27f76b08acc1160b5e"},"cell_type":"code","source":"# Since our labels are not mutually exclusive using softmax is not an option. As our first baseline we are going to \n# experiment with multi-task learning with hard parameter sharing, where in shared layers we are going to use convolutional \n# layers while task specific layers will be dense ones. Each task will have weighted binary cross-entropy loss function, \n# we will try to use weights to account for severe imbalance in the data.\n# Credit: coursera/convolutional_neural_networks/week_1/home_assignment\nfrom sklearn.metrics import f1_score\nfrom functools import partial\nimport tensorflow as tf\n\ndef create_placeholders(height=128, width=128, channels=4, n_targets=28):\n    ''' Creates the placeholders for the tensorflow session. '''\n    X = tf.placeholder(shape=[None,height,width,channels], dtype=tf.float32, name='X')\n    Y = tf.placeholder(shape=[None,n_targets], dtype=tf.float32, name='Y')    \n    is_train = tf.placeholder(tf.bool, name=\"is_train\")\n    return X, Y, is_train\n\n\n#-------------------------------------------- BASELINE ARCHITECTURE ------------------------------------------------\n\ndef initialize_parameters(channels=4):\n    ''' Initializes weight parameters to build a neural network with tensorflow'''\n    # We attempt to do architecture that closely resembles VGG19 by Karen Simonyan\n    # 1) 2 convolutional layers \n    W1 = tf.get_variable(name='W1',shape=[4,4,channels,24],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W2 = tf.get_variable(name='W2',shape=[4,4,24,32],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    # 2) Another 2 convolutional layers\n    W3 = tf.get_variable(name='W3',shape=[4,4,32,32],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W4 = tf.get_variable(name='W4',shape=[4,4,32,32],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    # 3) Another 2 convolutional layers\n    W5 = tf.get_variable(name='W5',shape=[4,4,32,64],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W6 = tf.get_variable(name='W6',shape=[4,4,64,64],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    return {\"W1\": W1,\"W2\": W2,\"W3\":W3,\"W4\":W4,'W5':W5,'W6':W6}\n\ndef forward_propagation_shared_layers(X, parameters, is_train):\n    \"\"\"\n    Implements forward propagation for shared layers, took inspiration from VGG16 model\n    by Karen Simonyan.\n    \"\"\"    \n    # Retrieve the parameters from the dictionary \"parameters\" \n    W1  = parameters['W1']\n    W2  = parameters['W2']\n    W3  = parameters['W3']\n    W4  = parameters['W4']\n    W5  = parameters['W5']\n    W6  = parameters['W6']\n    # BLOCK 1: CONV2D -> RELU -> CONV2D -> RELU -> MaxPool\n    layer1_Z1    = tf.nn.conv2d(X,W1,strides=[1,1,1,1],padding='SAME')\n    layer1_A1    = tf.nn.relu(layer1_Z1)\n    layer1_Z2    = tf.nn.conv2d(layer1_A1,W2,strides=[1,1,1,1],padding='SAME')\n    layer1_A2    = tf.nn.relu(layer1_Z2)\n    layer1_MP    = tf.nn.max_pool(layer1_A2,ksize=[1,4,4,1],strides=[1,4,4,1],padding='SAME')\n    # BLOCK 2: CONV2D -> RELU -> CONV2D -> RELU -> MaxPool\n    layer2_Z1    = tf.nn.conv2d(layer1_MP,W3,strides=[1,1,1,1],padding='SAME')\n    layer2_A1    = tf.nn.relu(layer2_Z1)\n    layer2_Z2    = tf.nn.conv2d(layer2_A1,W4,strides=[1,1,1,1],padding='SAME')\n    layer2_A2    = tf.nn.relu(layer2_Z2)\n    layer2_MP    = tf.nn.max_pool(layer2_A2,ksize=[1,4,4,1],strides=[1,4,4,1],padding='SAME')\n    # BLOCK 3: CONV2D -> RELU -> CONV2D -> RELU -> MaxPool\n    layer3_Z1    = tf.nn.conv2d(layer2_MP,W5,strides=[1,1,1,1],padding='SAME')\n    layer3_A1    = tf.nn.relu(layer3_Z1)\n    layer3_Z2    = tf.nn.conv2d(layer3_A1,W6,strides=[1,1,1,1],padding='SAME')\n    layer3_A2    = tf.nn.relu(layer3_Z2)\n    layer3_MP    = tf.nn.max_pool(layer3_A2,ksize=[1,4,4,1],strides=[1,4,4,1],padding='SAME')\n    # FLATTEN\n    layer4_FLAT  = tf.contrib.layers.flatten(layer3_MP)\n    # Fully connected layers\n    layer5_FC1   = tf.contrib.layers.fully_connected(layer4_FLAT, 108, activation_fn=tf.nn.relu)\n    layer6_FC2   = tf.contrib.layers.fully_connected(layer5_FC1, 56, activation_fn=tf.nn.relu)\n    return layer6_FC2\n    \ndef forward_propagation_task_specific_layers(shared_last_layer):\n    ''' Forward propagation for task specific layers '''\n    task_specific_layers = [0]*28\n    for i in range(28):\n        task_specific_layers[i] = tf.contrib.layers.fully_connected(shared_last_layer, 1, activation_fn=None)\n    return task_specific_layers\n\n#-------------------------------------------- ARCHITECTURE VERSION 1 ------------------------------------------------\n\ndef initialize_parameters_v1(channels=4):\n    ''' Initializes weight parameters to build a neural network with tensorflow'''\n    # We attempt to do architecture that closely resembles VGG19 by Karen Simonyan\n    # 1) 2 convolutional layers \n    W1 = tf.get_variable(name='W1',shape=[3,3,channels,8],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W2 = tf.get_variable(name='W2',shape=[3,3,8,8],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    # 2) Another 2 convolutional layers\n    W3 = tf.get_variable(name='W3',shape=[3,3,8,16],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W4 = tf.get_variable(name='W4',shape=[3,3,16,16],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    # 3) Another 2 convolutional layers\n    W5 = tf.get_variable(name='W5',shape=[3,3,16,32],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W6 = tf.get_variable(name='W6',shape=[3,3,32,32],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    # 4) Another 2 convolutional layers\n    W7 = tf.get_variable(name='W7',shape=[3,3,32,64],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    W8 = tf.get_variable(name='W8',shape=[3,3,64,128],initializer=tf.contrib.layers.xavier_initializer_conv2d())\n    return {\"W1\": W1,\"W2\": W2,\"W3\":W3,\"W4\":W4,'W5':W5,'W6':W6,'W7':W7,'W8':W8}\n\n\ndef forward_propagation_shared_layers_v1(X, parameters, is_train):\n    \"\"\"\n    Implements forward propagation for shared layers, took inspiration from VGG16 model\n    by Karen Simonyan.\n    \"\"\"    \n    # Retrieve the parameters from the dictionary \"parameters\" \n    W1  = parameters['W1']\n    W2  = parameters['W2']\n    W3  = parameters['W3']\n    W4  = parameters['W4']\n    W5  = parameters['W5']\n    W6  = parameters['W6']\n    W7  = parameters['W7']\n    # BLOCK 1: CONV2D -> RELU -> CONV2D -> BN -> RELU -> MaxPool\n    layer1_Z1    = tf.nn.conv2d(X,W1,strides=[1,1,1,1],padding='SAME')\n    layer1_A1    = tf.nn.relu(layer1_Z1)\n    layer1_Z2    = tf.nn.conv2d(layer1_A1,W2,strides=[1,1,1,1],padding='SAME')\n    layer1_BN    = tf.layers.batch_normalization(layer1_Z2,training=is_train)\n    layer1_A2    = tf.nn.relu(layer1_BN)\n    layer1_MP    = tf.nn.max_pool(layer1_A2,ksize=[1,2,2,1],strides=[1,4,4,1],padding='VALID')\n    # BLOCK 2: CONV2D -> RELU -> CONV2D -> BN -> RELU -> MaxPool\n    layer2_Z1    = tf.nn.conv2d(layer1_MP,W3,strides=[1,1,1,1],padding='SAME')\n    layer2_A1    = tf.nn.relu(layer2_Z1)\n    layer2_Z2    = tf.nn.conv2d(layer2_A1,W4,strides=[1,1,1,1],padding='SAME')\n    layer2_BN    = tf.layers.batch_normalization(layer2_Z2,training=is_train)\n    layer2_A2    = tf.nn.relu(layer2_BN)\n    layer2_MP    = tf.nn.max_pool(layer2_A2,ksize=[1,2,2,1],strides=[1,4,4,1],padding='VALID')\n    # BLOCK 3: CONV2D -> RELU -> CONV2D -> BN-> RELU -> MaxPool\n    layer3_Z1    = tf.nn.conv2d(layer2_MP,W5,strides=[1,1,1,1],padding='SAME')\n    layer3_A1    = tf.nn.relu(layer3_Z1)\n    layer3_Z2    = tf.nn.conv2d(layer3_A1,W6,strides=[1,1,1,1],padding='SAME')\n    layer3_BN    = tf.layers.batch_normalization(layer3_Z2,training=is_train)\n    layer3_A2    = tf.nn.relu(layer3_BN)\n    layer3_MP    = tf.nn.max_pool(layer3_A2,ksize=[1,2,2,1],strides=[1,4,4,1],padding='VALID')\n    # BLOCK 4: CONV2D -> RELU -> CONV2D -> BN-> RELU -> MaxPool\n    layer4_Z1    = tf.nn.conv2d(layer3_MP,W7,strides=[1,1,1,1],padding='SAME')\n    layer4_BN    = tf.layers.batch_normalization(layer4_Z1,training=is_train)\n    layer4_A2    = tf.nn.relu(layer4_BN)\n    layer4_MP    = tf.nn.max_pool(layer4_A2,ksize=[1,2,2,1],strides=[1,2,2,1],padding='VALID')\n    # FLATTEN\n    layer4_FLAT  = tf.contrib.layers.flatten(layer3_MP)\n    # Fully connected layers\n    layer5_FC1   = tf.contrib.layers.fully_connected(layer4_FLAT, 108, activation_fn=tf.nn.relu)\n    layer6_FC2   = tf.contrib.layers.fully_connected(layer5_FC1, 28, activation_fn=None)\n    return layer6_FC2\n\n#--------------------------------------- Loss Fuctions ----------------------------------------------\n\ndef multitask_heads(mlt,weights,labels):\n    ''' Headers for each task, computes predicted probabilities and task specific loss functions'''\n    predict_proba = [0]*28\n    loss = 0\n    for i,w in enumerate(weights):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label       = tf.expand_dims(tf.gather(labels,indices=i,axis=1),axis=1)\n        loss += tf.reduce_mean(tf.pow(tf.constant(w,dtype=tf.float32),labels[:,i])*tf.nn.sigmoid_cross_entropy_with_logits(logits=mlt[i], labels=task_label))\n    return predict_proba, loss\n                                                                \ndef multitask_heads_focal(mlt,weights,labels,gamma=2):\n    ''' Focal Loss '''\n    predict_proba = [0]*28\n    loss          = 0\n    for i in range(28):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label       = tf.expand_dims(tf.gather(labels,indices=i,axis=1),axis=1)\n        focal_loss_multiplier  = tf.where(tf.equal(task_label,tf.constant(1.,dtype=tf.float32)),tf.pow(1-predict_proba[i],gamma),tf.pow(predict_proba[i],gamma))\n        loss += tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=mlt[i], labels=task_label)*focal_loss_multiplier)\n    return predict_proba, loss   \n\ndef multitask_heads_weighted_focal(mlt,weights,labels,gamma=2):\n    ''' Weighted Focal Loss '''\n    predict_proba = [0]*28\n    loss          = 0\n    for i,w in enumerate(weights):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label       = tf.expand_dims(tf.gather(labels,indices=i,axis=1),axis=1)\n        focal_loss_multiplier  = tf.where(tf.equal(task_label,tf.constant(1.,dtype=tf.float32)),tf.pow(1-predict_proba[i],gamma),tf.pow(predict_proba[i],gamma))\n        weighting_multiplier   = tf.pow(tf.constant(w,dtype=tf.float32),labels[:,i])\n        loss += tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=mlt[i], labels=task_label)*focal_loss_multiplier*weighting_multiplier)\n    return predict_proba, loss  \n\ndef multitask_heads_weighted_focal_different_gammas(mlt,weights,labels,gammas):\n    ''' Focal loss with different gammas for different tasks'''\n    predict_proba = [0]*28\n    loss          = 0\n    for i,gamma in enumerate(gammas):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label       = tf.expand_dims(tf.gather(labels,indices=i,axis=1),axis=1)\n        focal_loss_multiplier = tf.where(tf.equal(task_label,tf.constant(1.,dtype=tf.float32)),tf.pow(1-predict_proba[i],gamma),tf.pow(predict_proba[i],gamma))\n        loss += tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=mlt[i], labels=task_label)*focal_loss_multiplier)\n    return predict_proba, loss\n        \ndef multitask_heads_weighted_hinge(mlt,weights,labels):\n    ''' Weighted Hinge Loss'''\n    predict_proba = [0]*28\n    labels_transformed = (labels + 1)/2\n    loss               = 0\n    for i,w in enumerate(weights):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label = tf.expand_dims(tf.gather(labels_transformed,indices=i,axis=1),axis=1)\n        weighting  = tf.pow(tf.constant(w,dtype=tf.float32),task_label)\n        loss      += tf.reduce_mean(weighting*tf.maximum(0., 1 - task_label*mlt[i]))\n    return predict_proba, loss\n\n\ndef multitask_heads_weighted_hinge_squared(mlt,weights,labels):\n    ''' Squared Weighted Hinge Loss '''\n    predict_proba = [0]*28\n    labels_transformed = (labels+1)/2\n    loss          = 0\n    for i,w in enumerate(weights):\n        predict_proba[i] = tf.nn.sigmoid(mlt[i])\n        task_label       = tf.expand_dims(tf.gather(labels_transformed,indices=i,axis=1),axis=1)\n        weighting        = tf.pow(tf.constant(w,dtype=tf.float32),task_label)\n        loss             = tf.reduce_mean(weighting*tf.square(tf.maximum(0.,1.-task_label*mlt[i])))\n    return predict_proba, loss\n\n\ndef _predict_probabilities(x_batch,X,sess,predict_proba,is_train):\n    ''' Helper function that is later wrapped in finctools.partial, for batch evaluation'''\n    return sess.run(predict_proba,feed_dict={X:x_batch,is_train:False})\n\ndef balanced_weights_per_class(targets):\n    ''' Computes balancing weights for imbalanced classes'''\n    pos_weights = np.sum(targets,axis=0)\n    neg_weights = targets.shape[0] - pos_weights\n    return neg_weights/pos_weights","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"18df2bcf64a4e9ae99ced282f8593c40fb49b221","scrolled":true},"cell_type":"code","source":"from functools import partial\nimport imgaug as ia\nfrom imgaug import augmenters as iaa\ntf.reset_default_graph()\n\nBATCH       = 27 # does not include augmentation\nRESIZE_ROWS = 3.2\nRESIZE_COLS = 3.2\nW, H        = 512, 512\nRESIZED_W   = int(W/RESIZE_ROWS)\nRESIZED_H   = int(H/RESIZE_COLS)\nCHANNELS    = ['green','red','blue','yellow']\nLOSSES      = ['WBCE','Focal','WeightedFocal','WeightedHinge','WeightedHingeSquared']\nLOSS_TO_USE = 'Focal'\n# Baseline architecture has more than 7 million weights, taking into account number of examples we have in our \n# dataset we suspect that we are overfitting\nARCHITECTURE = 'Baseline'\nAUGMENT     = True\n\n# Step 0: Compute balancing weights\nbalanced_weights = balanced_weights_per_class(labels_train)\n\n# Step 1: Model definition\nwith tf.device('/device:GPU:0'):\n    #tf.reset_default_graph()   \n    config = tf.ConfigProto()\n    config.gpu_options.allow_growth = True\n    sess = tf.Session(config=config)\n    X, Y, is_train             = create_placeholders(height=RESIZED_H, width=RESIZED_W,channels=len(CHANNELS))\n    \n    # choose architecture\n    if ARCHITECTURE == 'Baseline':\n        parameters             = initialize_parameters(channels=len(CHANNELS))\n        shared_layers          = forward_propagation_shared_layers(X, parameters, is_train)\n        mlt_logits             = forward_propagation_task_specific_layers(shared_layers)\n    elif ARCHITECTURE == 'Version1':\n        parameters             = initialize_parameters_v1(channels=len(CHANNELS))\n        shared_layers          = forward_propagation_shared_layers_v1(X, parameters, is_train)\n        mlt_logits             = forward_propagation_task_specific_layers(shared_layers)\n    elif ARCHITECTURE == 'Version2':\n        # this one is very similar to baseline however does not have task specific layers\n        parameters             = initialize_parameters(channels=len(CHANNELS))    \n        shared_layers          = forward_propagation_shared_layers(X, parameters, is_train)\n        last_layer             = tf.contrib.layers.fully_connected(shared_layers, 28, activation_fn=None)\n        mlt_logits             = [tf.expand_dims(tf.gather(last_layer,indices=i,axis=1),axis=-1) for i in range(28)]\n        \n    # choose loss \n    if LOSS_TO_USE == 'Focal':\n        predict_proba, loss    = multitask_heads_focal(mlt_logits,balanced_weights,Y, gamma=6)\n    if LOSS_TO_USE == 'FocalDifferentGammas':\n        gammas                 = compute_task_specific_gamma(labels)\n        predict_proba, loss    = multitask_heads_weighted_focal_different_gammas(mlt_logits,balanced_weights,Y,gammas)\n    elif LOSS_TO_USE == 'WBCE':\n        predict_proba, loss    = multitask_heads(mlt_logits,balanced_weights,Y)\n    elif LOSS_TO_USE == 'WeightedFocal':\n        predict_proba, loss    = multitask_heads_weighted_focal(mlt_logits,balanced_weights,Y)\n    elif LOSS_TO_USE == 'WeightedHinge':\n        predict_proba, loss    = multitask_heads_weighted_hinge(mlt_logits,balanced_weights,Y)\n    elif LOSS_TO_USE == 'WeightedHingeSquared':\n        predict_proba, loss    = multitask_heads_weighted_hinge_squared(mlt_logits,balanced_weights,Y)\n    \n    predict_probabilities  = partial(_predict_probabilities,X=X,sess=sess,predict_proba=predict_proba,is_train=is_train)\n    #loss                   = tf.reduce_sum(weighted_task_losses)\n    optimizer              = tf.train.AdamOptimizer(learning_rate=0.001)\n    # account for batch normalization layers\n    if ARCHITECTURE == 'Version1':\n        update_ops             = tf.get_collection(tf.GraphKeys.UPDATE_OPS)\n        with tf.control_dependencies(update_ops):\n            opt = optimizer.minimize(loss)\n    else:\n        opt = optimizer.minimize(loss)\n    init                   = tf.global_variables_initializer()\n    sess.run(init)\n\n# Step 2: Run Model for specified \nn_epochs = 15\nloss_values_by_iteration = []\nfor j in range(n_epochs):\n    #out = evaluate(paths_val,labels_val,predict_probabilities,batch_size=BATCH,row_resize_ratio=RESIZE_ROWS,\n    #               col_resize_ratio=RESIZE_COLS,width=W,height=H,channels_to_include=['green'])\n    #print('\\n Epoch eval {0} \\n'.format(out))\n    print('****** Epoch {0} started! *****'.format(j))\n    train_gen = train_generator_single_epoch(paths_train,labels_train,batch_size=BATCH,row_resize_ratio=RESIZE_ROWS,\n                                             col_resize_ratio=RESIZE_COLS,width=W,height=H,channels_to_include=CHANNELS)\n    for i,(batch_images,batch_labels) in enumerate(train_gen): \n        if AUGMENT:\n            augmenter = iaa.Sequential([ iaa.OneOf([ iaa.Affine(rotate=0), iaa.Affine(rotate=90), iaa.Affine(rotate=180),\n                                                     iaa.Affine(rotate=270),     iaa.Fliplr(0.5), iaa.Flipud(0.5)])], random_order=True)\n            batch_images = np.concatenate((batch_images, augmenter.augment_images(batch_images)), 0)\n            batch_labels = np.concatenate((batch_labels, batch_labels), 0)\n        # normalize images \n        batch_images = batch_images/255\n        _, loss_val = sess.run([opt,loss], feed_dict={X:batch_images,Y:batch_labels,is_train:True})\n        if i%100==0:\n            print('Iteration {0}, value of loss function {1}'.format(i,loss_val))\n        loss_values_by_iteration.append(loss_val)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"585a3cd13ba62c71fc2ebc61fadb203941956233","scrolled":true},"cell_type":"code","source":"probs = eval_predict_proba_generator(paths_val,predict_probabilities,batch_size=BATCH,row_resize_ratio=RESIZE_ROWS,\n                                     col_resize_ratio=RESIZE_COLS,width=W,height=H,channels_to_include=CHANNELS)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8b792f81625bbfb2a29d84900eb5361dccda7a79"},"cell_type":"code","source":"def predict_from_probs(predict_proba,cutoffs):\n    y_hat = np.zeros(predict_proba.shape)\n    for j,cutoff in enumerate(cutoffs):\n        y_hat[predict_proba[:,j]>cutoff,j]=1\n    return y_hat\n\ndef evaluate_given_thresholds(probs,thresholds,labels):\n    ''' Evaluate '''\n    y_preds   = np.zeros(probs.shape)\n    for j,thresh in enumerate(thresholds):\n        y_preds[probs[:,j]>thresh,j] = 1\n    f1_macro = f1_score(labels,y_preds,average='macro')\n    return f1_macro\n\n# 1) Strategy 1: Search single constant threshold\nbest_threshold  = None\nbest_macro_f1   = float('-inf')\nn_trials        = 100\nthresholds      = np.linspace(0,1,n_trials)\nfor i in range(n_trials):\n    threshold = thresholds[i]*np.ones(28)\n    f1_new = evaluate_given_thresholds(probs,threshold,labels_val)\n    if f1_new > best_macro_f1:\n        print('New best F-1 {0}'.format(f1_new))\n        best_threshold = threshold\n        best_macro_f1  = f1_new\nprint('Best F-1 macro achieved with first strategy is {0}'.format(best_macro_f1))\nprint(best_threshold)\n    \n# 2) Trying to improve upon startegy 1: Search individual thresholds that corrspond to proportion of labels in train set\ntrials_per_column = 100\ntrue_props      = np.sum(labels_train,axis=0) / labels_train.shape[0]\ncol_indexes     = np.argsort(true_props)\nthreshs = np.linspace(0,1,trials_per_column)\nfor col in col_indexes:\n    best_thresh_val = best_threshold[col]\n    for t in threshs:\n        best_threshold[col] = t\n        f1 = evaluate_given_thresholds(probs,best_threshold,labels_val)\n        if f1 > best_macro_f1:\n            best_thresh_val = t\n            best_macro_f1   = f1\n            print('New Best F-1 is {0}'.format(best_macro_f1))\n    best_threshold[col] = best_thresh_val    \n\n    \nprint('FINAL BEST F-1 is {0}'.format(evaluate_given_thresholds(probs,best_threshold,labels_val)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9d261110b73b852fbe6310ce7674303610ab78ba"},"cell_type":"markdown","source":"## Submission Preparation"},{"metadata":{"trusted":true,"_uuid":"38b4a1b8b36a29bce59b405dce023fe1a7fbad14"},"cell_type":"code","source":"import tqdm\n\nsubmit           = pd.read_csv('../input/sample_submission.csv')\npath_test_data,_ = test_data()\nsubmission       = []\nprobs_test_data  = eval_predict_proba_generator(path_test_data,predict_probabilities,batch_size=BATCH,row_resize_ratio=RESIZE_ROWS,\n                                                 col_resize_ratio=RESIZE_COLS,width=W,height=H,channels_to_include=CHANNELS)\nlabel_predict    = predict_from_probs(probs_test_data,best_threshold)\nfor j in range(label_predict.shape[0]):\n    pos_labels   = np.where(label_predict[j,:]==1)[0]\n    submission.append(' '.join(str(l) for l in pos_labels))\nsubmit['Predicted'] = submission\nnp.save('draw_predict_proba_baseline.npy', probs_test_data)\nsubmit.to_csv('submit_baseline.csv', index=False)\n# record loss function\nloss_function_values = np.asarray(loss_values_by_iteration)\nnp.save('loss_function_values.npy',loss_function_values)\nlv_save = np.asarray(loss_values_by_iteration)\nnp.save('loss_values.npy',lv_save)\nbest_f1 = np.array([best_f1])\nnp.save('val_best_f1.npy',best_f1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4cb05a91a704465a1ee004129e20044a85f78cd4"},"cell_type":"markdown","source":"## Additional Postprocessing"},{"metadata":{"trusted":true,"_uuid":"a5d2ce362ff0abfab0107544caf9658ad3cb8aca"},"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\n\nlabel_predictors = [0]*28\ninclusion_mask = np.array([True]*28)\nmodels = [0]*28\n\nfor j in range(28):\n    model = RandomForestClassifier\n    ","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}