{"cells":[{"metadata":{"_uuid":"1677fddbb95f7545b6540e9201f3339a0fdbfc5d","_cell_guid":"e61ef2d8-f315-4f7f-b07e-1de0f4e8441a"},"cell_type":"markdown","source":"**Author:** Raoul Malm\n\n**Abstract:** \n\nThe 2018 Data Science Bowl \"Find the nuclei in divergent images to advance medical discovery\" provides in its first stage a training and test data set consisting of 670 and 65 microscopic images of varying size showing ensembles of cells and their nuclei. For the training images the nuclei are segmented by humans such that we know their number and location within each image. The goal is to find the correct number and location of all nuclei shown in the test images. The performance of an algorithm is evaluated on the mean average precision at different intersection over union (IoU) thresholds, which will be referred to as the score in the following. \n\nFor the task we implement a deep neural network of the U-Net type consisting of several convolutional and max-pooling layers. The network is written in TensorFlow. Each model can be saved/loaded and the training process can be visualized with TensorBoard. The input of the network are images of shape (height, width, channels) while the output are corresponding binary masks of shape (height, width, 1). \n\nPrior to training the network we resize, normalize and transform the images. We use 10% of the training data for validation. Furthermore, we implement data augmentation by making use of translations, rotations, horizontal/vertical flipping and zoom. Choosing an image size of 384x384 pixels the network requires roughly 30 training epochs before the training seems to converge. The network achieves a score of 0.56/0.352 on the validation/test set. A major reason for the score discrepancy can be explained by overlapping/touching nuclei that are identified as a single nucleous by the current implementation. Furthermore, we have not tuned the hyperparameters, so there is still a lot of room for improvement. \n\nBeing constrained on kaggel hardware and running time, we can choose an image size of 256x256 pixels and train for 4 epochs, which takes roughly 40 minutes. The network can achieve a score of 0.45 on the validation set.\n\n**Outline:**\n\n1. [Modules and global settings](#1-bullet)\n2. [Analyse data](#2-bullet)\n3. [Manipulate data](#3-bullet)\n4. [Score Metric](#4-bullet)\n5. [Implement the Neural Network class](#5-bullet)\n6. [Train the Neural Network](#6-bullet)\n7. [Validate the Neural Network](#7-bullet)\n8. [Make Test Prediction](#8-bullet)\n9. [Submit](#9-bullet)\n\n**Reference:**\n\n[U-Net: Convolutional Networks for Biomedical Image Segmentation](https://arxiv.org/abs/1505.04597)\n","outputs":[],"execution_count":null},{"metadata":{"_uuid":"180bee0860395d3ce516f38d62814d8ebc5f061c","_cell_guid":"b992f278-0b93-4b68-93ae-9532f5eabb67"},"cell_type":"markdown","source":"# 1. Modules and global settings  <a class=\"anchor\" id=\"1-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"5c38504af3a84bee68c66d3cde74443c58df422f","_cell_guid":"c332549b-8d23-4bb5-8497-e7a8eb8b21d2","trusted":true},"cell_type":"code","source":"# Import necessary modules and set global constants and variables. \n\nimport tensorflow as tf            \nimport pandas as pd                 \nimport numpy as np                                       \nimport sklearn.model_selection     # For using KFold\nimport keras.preprocessing.image   # For using image generation\nimport datetime                    # To measure running time \nimport skimage.transform           # For resizing images\nimport skimage.morphology          # For using image labeling\nimport cv2                         # To read and manipulate images\nimport os                          # For filepath, directory handling\nimport sys                         # System-specific parameters and functions\nimport tqdm                        # Use smart progress meter\nimport seaborn as sns              # For pairplots\nimport matplotlib.pyplot as plt    # Python 2D plotting library\nimport matplotlib.cm as cm         # Color map\n%matplotlib inline                  \n\n# Global constants.\nIMG_WIDTH = 256       # Default image width\nIMG_HEIGHT = 256      # Default image height\nIMG_CHANNELS = 3      # Default number of channels\nCW_DIR = os.getcwd()  \nTRAIN_DIR = os.path.join(os.path.dirname(CW_DIR), 'input', 'stage1_train')\nTEST_DIR = os.path.join(os.path.dirname(CW_DIR), 'input', 'stage1_test')\nIMG_TYPE = '.png'         # Image type\nIMG_DIR_NAME = 'images'   # Folder name including the image\nMASK_DIR_NAME = 'masks'   # Folder name including the masks\nLOGS_DIR_NAME = 'logs'    # Folder name for TensorBoard summaries \nSAVES_DIR_NAME = 'saves'  # Folder name for storing network parameters\nSEED = 123                # Random seed for splitting train/validation sets\n    \n# Global variables.\nmin_object_size = 1       # Minimal nucleous size in pixels\nx_train = []\ny_train = []\nx_test = []\ny_test_pred_proba = {}\ny_test_pred = {}\n\n# Display working/train/test directories.\nprint('CW_DIR = {}'.format(CW_DIR))\nprint('TRAIN_DIR = {}'.format(TRAIN_DIR))\nprint('TEST_DIR = {}'.format(TEST_DIR))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8fcf5470931c0d3cfc744ff09f50c0aa1c1c1b2a","_cell_guid":"b878becd-e3d4-42e3-8f3d-b33e4f152e48"},"cell_type":"markdown","source":"# 2. Analyse data  <a class=\"anchor\" id=\"2-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"6f01ddea27c26f0eec0564512093664612ff31eb","collapsed":true,"_cell_guid":"7d3fdaed-e6a9-4b60-9a02-d1c6324d8f89","trusted":true},"cell_type":"code","source":"# Collection of methods for data operations. Implemented are functions to read  \n# images/masks from files and to read basic properties of the train/test data sets.\n\ndef read_image(filepath, color_mode=cv2.IMREAD_COLOR, target_size=None):\n    \"\"\"Read an image from a file and resize it.\"\"\"\n    img = cv2.imread(filepath, color_mode)\n    if target_size: \n        img = cv2.resize(img, target_size, interpolation = cv2.INTER_AREA)\n    return img\n\ndef read_mask(directory, target_size=None):\n    \"\"\"Read and resize masks contained in a given directory.\"\"\"\n    for i,filename in enumerate(next(os.walk(directory))[2]):\n        mask_path = os.path.join(directory, filename)\n        mask_tmp = read_image(mask_path, cv2.IMREAD_GRAYSCALE, target_size)\n        if not i: mask = mask_tmp\n        else: mask = np.maximum(mask, mask_tmp)\n    return mask \n\ndef read_train_data_properties(train_dir, img_dir_name, mask_dir_name):\n    \"\"\"Read basic properties of training images and masks\"\"\"\n    tmp = []\n    for i,dir_name in enumerate(next(os.walk(train_dir))[1]):\n\n        img_dir = os.path.join(train_dir, dir_name, img_dir_name)\n        mask_dir = os.path.join(train_dir, dir_name, mask_dir_name) \n        num_masks = len(next(os.walk(mask_dir))[2])\n        img_name = next(os.walk(img_dir))[2][0]\n        img_name_id = os.path.splitext(img_name)[0]\n        img_path = os.path.join(img_dir, img_name)\n        img_shape = read_image(img_path).shape\n        tmp.append(['{}'.format(img_name_id), img_shape[0], img_shape[1],\n                    img_shape[0]/img_shape[1], img_shape[2], num_masks,\n                    img_path, mask_dir])\n\n    train_df = pd.DataFrame(tmp, columns = ['img_id', 'img_height', 'img_width',\n                                            'img_ratio', 'num_channels', \n                                            'num_masks', 'image_path', 'mask_dir'])\n    return train_df\n\ndef read_test_data_properties(test_dir, img_dir_name):\n    \"\"\"Read basic properties of test images.\"\"\"\n    tmp = []\n    for i,dir_name in enumerate(next(os.walk(test_dir))[1]):\n\n        img_dir = os.path.join(test_dir, dir_name, img_dir_name)\n        img_name = next(os.walk(img_dir))[2][0]\n        img_name_id = os.path.splitext(img_name)[0]\n        img_path = os.path.join(img_dir, img_name)\n        img_shape = read_image(img_path).shape\n        tmp.append(['{}'.format(img_name_id), img_shape[0], img_shape[1],\n                    img_shape[0]/img_shape[1], img_shape[2], img_path])\n\n    test_df = pd.DataFrame(tmp, columns = ['img_id', 'img_height', 'img_width',\n                                           'img_ratio', 'num_channels', 'image_path'])\n    return test_df\n\ndef imshow_args(x):\n    \"\"\"Matplotlib imshow arguments for plotting.\"\"\"\n    if len(x.shape)==2: return x, cm.gray\n    if x.shape[2]==1: return x[:,:,0], cm.gray\n    return x, None\n\ndef load_raw_data(image_size=(IMG_HEIGHT, IMG_WIDTH)):\n    \"\"\"Load raw data.\"\"\"\n    # Python lists to store the training images/masks and test images.\n    x_train, y_train, x_test = [],[],[]\n\n    # Read and resize train images/masks. \n    print('Loading and resizing train images and masks ...')\n    sys.stdout.flush()\n    for i, filename in tqdm.tqdm(enumerate(train_df['image_path']), total=len(train_df)):\n        img = read_image(train_df['image_path'].loc[i], target_size=image_size)\n        mask = read_mask(train_df['mask_dir'].loc[i], target_size=image_size)\n        x_train.append(img)\n        y_train.append(mask)\n\n    # Read and resize test images. \n    print('Loading and resizing test images ...')\n    sys.stdout.flush()\n    for i, filename in tqdm.tqdm(enumerate(test_df['image_path']), total=len(test_df)):\n        img = read_image(test_df['image_path'].loc[i], target_size=image_size)\n        x_test.append(img)\n\n    # Transform lists into 4-dim numpy arrays.\n    x_train = np.array(x_train)\n    y_train = np.expand_dims(np.array(y_train), axis=4)\n    x_test = np.array(x_test)\n\n    print('x_train.shape: {} of dtype {}'.format(x_train.shape, x_train.dtype))\n    print('y_train.shape: {} of dtype {}'.format(y_train.shape, x_train.dtype))\n    print('x_test.shape: {} of dtype {}'.format(x_test.shape, x_test.dtype))\n    \n    return x_train, y_train, x_test\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"45dbdf0f65e54707543059026edcde2679154300","_cell_guid":"2d0fcfce-732d-4a6c-a1ed-1ca8fff50e68","trusted":true},"cell_type":"code","source":"# Basic properties of images/masks. \ntrain_df = read_train_data_properties(TRAIN_DIR, IMG_DIR_NAME, MASK_DIR_NAME)\ntest_df = read_test_data_properties(TEST_DIR, IMG_DIR_NAME)\nprint('train_df:')\nprint(train_df.describe())\nprint('')\nprint('test_df:')\nprint(test_df.describe())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e458bd1eaff5cd0179f240fb94f0a11d4b2daad5","_cell_guid":"dce96692-18bb-4eee-99e8-bfa304f83ee6","trusted":false,"collapsed":true},"cell_type":"code","source":"# Counting unique image shapes.\ndf = pd.DataFrame([[x] for x in zip(train_df['img_height'], train_df['img_width'])])\ndf[0].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1d3fd7275aa001854f509041af956035aa9866e","_cell_guid":"3e6a4997-cae2-406e-8b01-9d568cd4cf13","trusted":false,"collapsed":true},"cell_type":"code","source":"# Overview of train images/masks. There is a lot of variation concerning\n# the form/size/number of nuclei and the darkness/lightness/colorfulness of \n# the images. \nfig, axs = plt.subplots(4,4,figsize=(20,20))\nfor i in range(4):\n    for j in range(2):\n        n = np.random.randint(0,len(train_df))\n        axs[i,j*2].imshow(read_image(train_df['image_path'].loc[n]))\n        axs[i,j*2].set_title('{}. image'.format(n))\n        axs[i,j*2+1].imshow(read_mask(train_df['mask_dir'].loc[n]), cmap='gray') \n        axs[i,j*2+1].set_title('{}. mask'.format(n))       ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7794ecc749aefa2e6d74b1b204bb09b7cb687789","_cell_guid":"a13030a7-118b-44ee-856d-2205c0edbb66","trusted":false,"collapsed":true},"cell_type":"code","source":"# Read images/masks from files and resize them. Each image and mask \n# is stored as a 3-dim array where the number of channels is 3 and 1, respectively.\nx_train, y_train, x_test = load_raw_data()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c5df940637b8be98a0edb08d073633bbefbceb80","_cell_guid":"f1b38887-bf0a-46ed-8d04-c6051db5e8f8","trusted":false,"collapsed":true},"cell_type":"code","source":"# Study the pixel intensity. On average the red, green and blue channels have similar\n# intensities for all images. It should be noted that the background can be dark \n# (black) as  as well as light (white). \ndef img_intensity_pairplot(x):\n    \"\"\"Plot intensity distributions of color channels.\"\"\"\n    df = pd.DataFrame()\n    df['Gray'] = np.mean(x[:,:,:,:], axis=(1,2,3))\n    if x.shape[3]==3:\n        df['Red'] = np.mean(x[:,:,:,0], axis=(1,2))\n        df['Blue'] = np.mean(x[:,:,:,1], axis=(1,2))\n        df['Green'] = np.mean(x[:,:,:,2], axis=(1,2))\n    return df\n\ncolor_df = img_intensity_pairplot(np.concatenate([x_train, x_test]))\ncolor_df['images'] = ['train']*len(x_train) + ['test']*len(x_test)\nsns.pairplot(color_df, hue = 'images');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e12809a76b5d4bd7c52d0c7359d768bb1dbce71f","_cell_guid":"68c9e4e8-9486-4d32-bd48-aa3ef7a3ee39"},"cell_type":"markdown","source":"# 3. Manipulate data <a class=\"anchor\" id=\"3-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"67bdcee8a2eddfa76da62313dbd01fa23a8fa594","collapsed":true,"_cell_guid":"ed05dea5-07aa-4757-a686-83c028485273","trusted":false},"cell_type":"code","source":"# Collection of methods for basic data manipulation like normalizing, inverting, \n# color transformation and generating new images/masks\n\ndef normalize_imgs(data):\n    \"\"\"Normalize images.\"\"\"\n    return normalize(data, type_=1)\n\ndef normalize_masks(data):\n    \"\"\"Normalize masks.\"\"\"\n    return normalize(data, type_=1)\n    \ndef normalize(data, type_=1): \n    \"\"\"Normalize data.\"\"\"\n    if type_==0:\n        # Convert pixel values from [0:255] to [0:1] by global factor\n        data = data.astype(np.float32) / data.max()\n    if type_==1:\n        # Convert pixel values from [0:255] to [0:1] by local factor\n        div = data.max(axis=tuple(np.arange(1,len(data.shape))), keepdims=True) \n        div[div < 0.01*data.mean()] = 1. # protect against too small pixel intensities\n        data = data.astype(np.float32)/div\n    if type_==2:\n        # Standardisation of each image \n        data = data.astype(np.float32) / data.max() \n        mean = data.mean(axis=tuple(np.arange(1,len(data.shape))), keepdims=True) \n        std = data.std(axis=tuple(np.arange(1,len(data.shape))), keepdims=True) \n        data = (data-mean)/std\n\n    return data\n\ndef trsf_proba_to_binary(y_data):\n    \"\"\"Transform propabilities into binary values 0 or 1.\"\"\"  \n    return np.greater(y_data,.5).astype(np.uint8)\n\ndef invert_imgs(imgs, cutoff=.5):\n    '''Invert image if mean value is greater than cutoff.'''\n    imgs = np.array(list(map(lambda x: 1.-x if np.mean(x)>cutoff else x, imgs)))\n    return normalize_imgs(imgs)\n\ndef imgs_to_grayscale(imgs):\n    '''Transform RGB images into grayscale spectrum.''' \n    if imgs.shape[3]==3:\n        imgs = normalize_imgs(np.expand_dims(np.mean(imgs, axis=3), axis=3))\n    return imgs\n\ndef generate_images(imgs, seed=None):\n    \"\"\"Generate new images.\"\"\"\n    # Transformations.\n    image_generator = keras.preprocessing.image.ImageDataGenerator(\n        rotation_range = 90., width_shift_range = 0.02 , height_shift_range = 0.02,\n        zoom_range = 0.10, horizontal_flip=True, vertical_flip=True)\n    \n    # Generate new set of images\n    imgs = image_generator.flow(imgs, np.zeros(len(imgs)), batch_size=len(imgs),\n                                shuffle = False, seed=seed).next()    \n    return imgs[0]\n\ndef generate_images_and_masks(imgs, masks):\n    \"\"\"Generate new images and masks.\"\"\"\n    seed = np.random.randint(10000) \n    imgs = generate_images(imgs, seed=seed)\n    masks = trsf_proba_to_binary(generate_images(masks, seed=seed))\n    return imgs, masks\n\ndef preprocess_raw_data(x_train, y_train, x_test, grayscale=False, invert=False):\n    \"\"\"Preprocessing of images and masks.\"\"\"\n    # Normalize images and masks\n    x_train = normalize_imgs(x_train)\n    y_train = trsf_proba_to_binary(normalize_masks(y_train))\n    x_test = normalize_imgs(x_test)\n    print('Images normalized.')\n \n    if grayscale:\n        # Remove color and transform images into grayscale spectrum.\n        x_train = imgs_to_grayscale(x_train)\n        x_test = imgs_to_grayscale(x_test)\n        print('Images transformed into grayscale spectrum.')\n\n    if invert:\n        # Invert images, such that each image has a dark background.\n        x_train = invert_imgs(x_train)\n        x_test = invert_imgs(x_test)\n        print('Images inverted to remove light backgrounds.')\n\n    return x_train, y_train, x_test\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b95c8b74cc603b9426fe277ff0c92485bce4828f","_cell_guid":"dc5dc1e5-38a2-4908-8042-aa3b9055d8d9","trusted":false,"collapsed":true},"cell_type":"code","source":"# Normalize all images and masks. There is the possibility to transform images \n# into the grayscale sepctrum and to invert images which have a very \n# light background.\nx_train, y_train, x_test = preprocess_raw_data(x_train, y_train, x_test, invert=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4c3b9b1076eabca43e31149cc35e9678985d1079","_cell_guid":"7263f364-58a1-4e31-95e1-78cc0c4e66da","trusted":false,"collapsed":true},"cell_type":"code","source":"# Analyze nuclei sizes.\ndef get_nuclei_sizes():\n    nuclei_sizes = []\n    mask_idx = []\n    for i in range(len(y_train)):\n        mask = y_train[i].reshape(y_train.shape[1], y_train.shape[2])\n        lab_mask = skimage.morphology.label(mask > .5)\n        (mask_labels, mask_sizes) = np.unique(lab_mask, return_counts=True)\n        nuclei_sizes.extend(mask_sizes[1:])\n        mask_idx.extend([i]*len(mask_sizes[1:]))\n    return mask_idx, nuclei_sizes\n\nmask_idx, nuclei_sizes = get_nuclei_sizes()\nnuclei_sizes_df = pd.DataFrame()   \nnuclei_sizes_df['mask_index'] = mask_idx\nnuclei_sizes_df['nucleous_size'] = nuclei_sizes\n\nprint(nuclei_sizes_df.describe())\nnuclei_sizes_df.sort_values(by='nucleous_size', ascending=True).head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"abd0470203eea27e58679140e0cfe08f6e0147fc","_cell_guid":"721b38ad-ae2c-45ac-b4af-aa23cfcae0cd","trusted":false,"collapsed":true},"cell_type":"code","source":"color_df = img_intensity_pairplot(np.concatenate([x_train, x_test]))\ncolor_df['images'] = ['train']*len(x_train) + ['test']*len(x_test)\nsns.pairplot(color_df, hue = 'images');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6ed2c5ee02f0ae007eca90b0adebac6717a1b6ef","_cell_guid":"2344f4e1-84a3-45cc-80ab-aa4a50a30272","trusted":false,"collapsed":true},"cell_type":"code","source":"# Check the image transformation procedure (resizing, normalizing, inverting) \n# by looking at a sample. \ndef img_comparison_plot(n):\n    \"\"\"Plot the original and transformed images/masks.\"\"\"\n    fig, axs = plt.subplots(1,4,figsize=(20,20))\n    axs[0].imshow(read_image(train_df['image_path'].loc[n]))\n    axs[0].set_title('{}.) original image'.format(n))\n    img, img_type = imshow_args(x_train[n])\n    axs[1].imshow(img, img_type)\n    axs[1].set_title('{}.) transformed image'.format(n))\n    axs[2].imshow(read_mask(train_df['mask_dir'].loc[n]), cm.gray) \n    axs[2].set_title('{}.) original mask'.format(n))\n    axs[3].imshow(y_train[n,:,:,0], cm.gray)\n    axs[3].set_title('{}.) transformed mask'.format(n));\n\nn = 617 # np.random.randint(0, len(x_train))\nimg_comparison_plot(n)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"66a47fdfe4dec54198c99237c4ddbd2f1b372800","_cell_guid":"e3a084e0-d092-4896-a6aa-8971a7c767ec","trusted":false,"collapsed":true},"cell_type":"code","source":"# Generate new images/masks via transformations applied on the original \n# images/maks. Data augmentations can be used for regularization.\ndef plot_generated_image_mask(n):\n    fig, axs = plt.subplots(1,4,figsize=(20,20))\n    img_new, mask_new = generate_images_and_masks(x_train[n:n+1], y_train[n:n+1])\n    img, img_type = imshow_args(x_train[n])\n    axs[0].imshow(img, img_type)\n    axs[0].set_title('{}. original image'.format(n))\n    img, img_type = imshow_args(img_new[0])\n    axs[1].imshow(img, img_type)\n    axs[1].set_title('{}. generated image'.format(n))\n    axs[2].imshow(y_train[n,:,:,0], cmap='gray')\n    axs[2].set_title('{}. original mask'.format(n))\n    axs[3].imshow(mask_new[0,:,:,0], cmap='gray')\n    axs[3].set_title('{}. generated mask'.format(n));\n\nn = 166 # np.random.randint(len(x_train))\nplot_generated_image_mask(n)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c5c0eec3d5f75d6f5bd57a3a80ca0b2c9a821559","_cell_guid":"cfdf240c-a2fd-4a16-8cfc-98edee5fa420"},"cell_type":"markdown","source":"# 4. Score Metric <a class=\"anchor\" id=\"4-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"7db4e41c6c4a36707b9b24cfcf3ffdac48334d13","collapsed":true,"_cell_guid":"e6a31484-94c3-4efc-8399-cf275cdce49a","trusted":false},"cell_type":"code","source":"\"\"\" Collection of methods to compute the score.\n\n1. We start with a true and predicted mask, corresponding to one train image.\n\n2. The true mask is segmented into different objects. Here lies a main source \nof error. Overlapping or touching nuclei are not separated but are labeled as \none object. This means that the target mask can contain less objects than \nthose that have been originally identified by humans.\n\n3. In the same manner the predicted mask is segmented into different objects.\n\n4. We compute all intersections between the objects of the true and predicted \nmasks. Starting with the largest intersection area we assign true objects to \npredicted ones, until there are no true/pred objects left that overlap. \nWe then compute for each true/pred object pair their corresponding intersection \nover union (iou) ratio. \n\n5. Given some threshold t we count the object pairs that have an iou > t, which\nyields the number of true positives: tp(t). True objects that have no partner are \ncounted as false positives: fp(t). Likewise, predicted objects without a counterpart\na counted as false negatives: fn(t).\n\n6. Now, we compute the precision tp(t)/(tp(t)+fp(t)+fn(t)) for t=0.5,0.55,0.60,...,0.95\nand take the mean value as the final precision (score).\n\"\"\"\n\ndef get_labeled_mask(mask, cutoff=.5):\n    \"\"\"Object segmentation by labeling the mask.\"\"\"\n    mask = mask.reshape(mask.shape[0], mask.shape[1])\n    lab_mask = skimage.morphology.label(mask > cutoff) \n    \n    # Keep only objects that are large enough.\n    (mask_labels, mask_sizes) = np.unique(lab_mask, return_counts=True)\n    if (mask_sizes < min_object_size).any():\n        mask_labels = mask_labels[mask_sizes < min_object_size]\n        for n in mask_labels:\n            lab_mask[lab_mask == n] = 0\n        lab_mask = skimage.morphology.label(lab_mask > cutoff) \n    \n    return lab_mask  \n\ndef get_iou(y_true_labeled, y_pred_labeled):\n    \"\"\"Compute non-zero intersections over unions.\"\"\"\n    # Array of different objects and occupied area.\n    (true_labels, true_areas) = np.unique(y_true_labeled, return_counts=True)\n    (pred_labels, pred_areas) = np.unique(y_pred_labeled, return_counts=True)\n\n    # Number of different labels.\n    n_true_labels = len(true_labels)\n    n_pred_labels = len(pred_labels)\n\n    # Each mask has at least one identified object.\n    if (n_true_labels > 1) and (n_pred_labels > 1):\n        \n        # Compute all intersections between the objects.\n        all_intersections = np.zeros((n_true_labels, n_pred_labels))\n        for i in range(y_true_labeled.shape[0]):\n            for j in range(y_true_labeled.shape[1]):\n                m = y_true_labeled[i,j]\n                n = y_pred_labeled[i,j]\n                all_intersections[m,n] += 1 \n\n        # Assign predicted to true background.\n        assigned = [[0,0]]\n        tmp = all_intersections.copy()\n        tmp[0,:] = -1\n        tmp[:,0] = -1\n\n        # Assign predicted to true objects if they have any overlap.\n        for i in range(1, np.min([n_true_labels, n_pred_labels])):\n            mn = list(np.unravel_index(np.argmax(tmp), (n_true_labels, n_pred_labels)))\n            if all_intersections[mn[0], mn[1]] > 0:\n                assigned.append(mn)\n            tmp[mn[0],:] = -1\n            tmp[:,mn[1]] = -1\n        assigned = np.array(assigned)\n\n        # Intersections over unions.\n        intersection = np.array([all_intersections[m,n] for m,n in assigned])\n        union = np.array([(true_areas[m] + pred_areas[n] - all_intersections[m,n]) \n                           for m,n in assigned])\n        iou = intersection / union\n\n        # Remove background.\n        iou = iou[1:]\n        assigned = assigned[1:]\n        true_labels = true_labels[1:]\n        pred_labels = pred_labels[1:]\n\n        # Labels that are not assigned.\n        true_not_assigned = np.setdiff1d(true_labels, assigned[:,0])\n        pred_not_assigned = np.setdiff1d(pred_labels, assigned[:,1])\n        \n    else:\n        # in case that no object is identified in one of the masks\n        iou = np.array([])\n        assigned = np.array([])\n        true_labels = true_labels[1:]\n        pred_labels = pred_labels[1:]\n        true_not_assigned = true_labels\n        pred_not_assigned = pred_labels\n        \n    # Returning parameters.\n    params = {'iou': iou, 'assigned': assigned, 'true_not_assigned': true_not_assigned,\n             'pred_not_assigned': pred_not_assigned, 'true_labels': true_labels,\n             'pred_labels': pred_labels}\n    return params\n\ndef get_score_summary(y_true, y_pred):\n    \"\"\"Compute the score for a single sample including a detailed summary.\"\"\"\n    \n    y_true_labeled = get_labeled_mask(y_true)  \n    y_pred_labeled = get_labeled_mask(y_pred)  \n    \n    params = get_iou(y_true_labeled, y_pred_labeled)\n    iou = params['iou']\n    assigned = params['assigned']\n    true_not_assigned = params['true_not_assigned']\n    pred_not_assigned = params['pred_not_assigned']\n    true_labels = params['true_labels']\n    pred_labels = params['pred_labels']\n    n_true_labels = len(true_labels)\n    n_pred_labels = len(pred_labels)\n\n    summary = []\n    for i,threshold in enumerate(np.arange(0.5, 1.0, 0.05)):\n        tp = np.sum(iou > threshold)\n        fn = n_true_labels - tp\n        fp = n_pred_labels - tp\n        if (tp+fp+fn)>0: \n            prec = tp/(tp+fp+fn)\n        else: \n            prec = 0\n        summary.append([threshold, prec, tp, fp, fn])\n\n    summary = np.array(summary)\n    score = np.mean(summary[:,1]) # Final score.\n    params_dict = {'summary': summary, 'iou': iou, 'assigned': assigned, \n                   'true_not_assigned': true_not_assigned, \n                   'pred_not_assigned': pred_not_assigned, 'true_labels': true_labels,\n                   'pred_labels': pred_labels, 'y_true_labeled': y_true_labeled,\n                   'y_pred_labeled': y_pred_labeled}\n    \n    return score, params_dict\n\ndef get_score(y_true, y_pred):\n    \"\"\"Compute the score for a batch of samples.\"\"\"\n    scores = []\n    for i in range(len(y_true)):\n        score,_ = get_score_summary(y_true[i], y_pred[i])\n        scores.append(score)\n    return np.array(scores)\n\ndef plot_score_summary(y_true, y_pred):\n    \"\"\"Plot score summary for a single sample.\"\"\"\n    # Compute score and assign parameters.\n    score, params_dict = get_score_summary(y_true, y_pred)\n    \n    assigned = params_dict['assigned']\n    true_not_assigned = params_dict['true_not_assigned']\n    pred_not_assigned = params_dict['pred_not_assigned']\n    true_labels = params_dict['true_labels']\n    pred_labels = params_dict['pred_labels']\n    y_true_labeled = params_dict['y_true_labeled']\n    y_pred_labeled = params_dict['y_pred_labeled']\n    summary = params_dict['summary']\n\n    n_assigned = len(assigned)\n    n_true_not_assigned = len(true_not_assigned)\n    n_pred_not_assigned = len(pred_not_assigned)\n    n_true_labels = len(true_labels)\n    n_pred_labels = len(pred_labels)\n\n    # Summary dataframe.\n    summary_df = pd.DataFrame(summary,columns=['threshold','precision','tp','fp','fn'])\n    print('Final score:', score)\n    print(summary_df)\n\n    # Plots.\n    fig, axs = plt.subplots(2,3,figsize=(20,13))\n\n    # True mask with true objects.\n    img = y_true\n    axs[0,0].imshow(img, cmap=cm.gray)\n    axs[0,0].set_title('{}.) true mask: {} true objects'.format(n,train_df['num_masks'][n]))\n    \n    # True mask with identified objects.\n    #img = np.zeros(y_true.shape)\n    #img[y_true_labeled > 0.5] = 255\n    img, img_type = imshow_args(y_true_labeled)\n    axs[0,1].imshow(img, img_type)\n    axs[0,1].set_title('{}.) true mask: {} objects identified'.format(n, n_true_labels))\n    \n    # Predicted mask with identified objects.\n    #img = np.zeros(y_true.shape)\n    #img[y_pred_labeled > 0.5] = 255\n    img, img_type = imshow_args(y_pred_labeled)\n    axs[0,2].imshow(img, img_type)\n    axs[0,2].set_title('{}.) predicted mask: {} objects identified'.format(\n        n, n_pred_labels))\n\n    # Prediction overlap with true mask.\n    img = np.zeros(y_true.shape)\n    img[y_true > 0.5] = 100\n    for i,j in assigned: img[(y_true_labeled == i) & (y_pred_labeled == j)] = 255\n    axs[1,0].set_title('{}.) {} pred. overlaps (white) with true objects (gray)'.format(\n        n,len(assigned)))\n    axs[1,0].imshow(img, cmap='gray', norm=None)\n\n    # Intersection over union.\n    img = np.zeros(y_true.shape)\n    img[(y_pred_labeled > 0) & (y_pred_labeled < 100)] = 100\n    img[(y_true_labeled > 0) & (y_true_labeled < 100)] = 100\n    for i,j in assigned: img[(y_true_labeled == i) & (y_pred_labeled == j)] = 255\n    axs[1,1].set_title('{}.) {} intersections (white) over unions (gray)'.format(\n        n, n_assigned))\n    axs[1,1].imshow(img, cmap='gray');\n\n    # False positives and false negatives.\n    img = np.zeros(y_true.shape)\n    for i in pred_not_assigned: img[(y_pred_labeled == i)] = 255\n    for i in true_not_assigned: img[(y_true_labeled == i)] = 100\n    axs[1,2].set_title('{}.) no threshold: {} fp (white), {} fn (gray)'.format(\n        n, n_pred_not_assigned, n_true_not_assigned))\n    axs[1,2].imshow(img, cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"daab06a90fdae02067485444cd385bbf644b56ba","_cell_guid":"77927f19-9c2c-492d-b7c7-d32360064723","trusted":false,"collapsed":true},"cell_type":"code","source":"# Check the score metric for one sample. The predicted mask is simulated\n# and can be modified in order to check the correct implementation of\n# the score metric.\nn = 166 # np.random.randint(len(x_train))\ntrue_mask = y_train[n,:,:,0].copy()\nlab_true_mask = get_labeled_mask(true_mask)\npred_mask = true_mask.copy() # Create predicted mask from true mask.\ntrue_mask[lab_true_mask == 7] = 0 # Remove one object => false postive\npred_mask[lab_true_mask == 10] = 0 # Remove one object => false negative\noffset = 5  # Offset.\npred_mask = pred_mask[offset:, offset:]\npred_mask = np.pad(pred_mask, ((0, offset), (0, offset)), mode=\"constant\")\nplot_score_summary(true_mask, pred_mask) ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"41de2084d9266ff8c3ab682d8786b457cea504c0","_cell_guid":"0e8bfc9a-9ea7-4d1d-9b8e-54e6e65620a1","trusted":false,"collapsed":true},"cell_type":"code","source":"# Study how many objects in the masks can be identified. This is a limiting factor\n# for the overall performance.\nmin_pixels_per_object = 1\nsummary = []\nfor n in range(len(y_train)):\n    img = y_train[n,:,:,0]\n    lab_img=get_labeled_mask(img)\n    img_labels, img_area = np.unique(lab_img, return_counts=True)\n    img_labels = img_labels[img_area>=min_pixels_per_object]\n    img_area = img_area[img_area>=min_pixels_per_object]\n    n_true_labels = train_df['num_masks'][n]\n    n_ident_labels = len(img_labels)\n    diff = np.abs(n_ident_labels-n_true_labels)\n    summary.append([n_true_labels, n_ident_labels, diff])\n\nsum_df = pd.DataFrame(summary, columns=(['true_objects', 'identified_objects', 'subtraction']))\nsum_df.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"56d8d63223a7eb3c818a691f2431896ca7c5d967","_cell_guid":"910e2b9b-f16a-4d89-a20e-d4d7c9327d59"},"cell_type":"markdown","source":"# 5. Implement the Neural Network Class  <a class=\"anchor\" id=\"5-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"f1b41a418790f009b4f66921431ad3d734b04024","collapsed":true,"_cell_guid":"74f13201-4ce9-44c9-bed8-aa723e9c6ad2","trusted":false},"cell_type":"code","source":"class NeuralNetwork():\n    \"\"\" Implements a neural network.\n        \n        TensorFlow is used to implement the U-Net, which consists of convolutional\n        and max pooling layers. Input and output shapes coincide. Methods are\n        implemented to train the model, to save/load the complete session and to \n        attach summaries for visualization with TensorBoard. \n    \"\"\"\n\n    def __init__(self, nn_name='tmp', nn_type='UNet', log_step=0.2, keep_prob=0.33, \n                 mb_size=16, input_shape=[IMG_HEIGHT,IMG_WIDTH,IMG_CHANNELS], \n                 output_shape=[IMG_HEIGHT,IMG_WIDTH,1]):\n        \"\"\"Instance constructor.\"\"\"\n        \n        # Tunable hyperparameters for training.\n        self.mb_size = mb_size       # Mini batch size\n        self.keep_prob = keep_prob   # Keeping probability with dropout regularization \n        self.learn_rate_step = 3     # Step size in terms of epochs\n        self.learn_rate_alpha = 0.25 # Reduction of learn rate for each step \n        self.learn_rate_0 = 0.001    # Starting learning rate \n        self.dropout_proba = 0.1     # == 1-keep_probability\n        \n        # Set helper variables.\n        self.input_shape = input_shape\n        self.output_shape = output_shape\n        self.nn_type = nn_type                # Type of neural network\n        self.nn_name = nn_name                # Name of neural network\n        self.params = {}                      # For storing parameters\n        self.learn_rate_pos = 0                \n        self.learn_rate = self.learn_rate_0\n        self.index_in_epoch = 0 \n        self.epoch = 0. \n        self.log_step = log_step              # Log results in terms of epochs\n        self.n_log_step = 0                   # Count number of mini batches  \n        self.train_on_augmented_data = False  # True = use augmented data \n        self.use_tb_summary = False           # True = use TensorBoard summaries\n        self.use_tf_saver = False             # True = save the session\n        \n        # Parameters that should be stored.\n        self.params['train_loss']=[]\n        self.params['valid_loss']=[]\n        self.params['train_score']=[]\n        self.params['valid_score']=[]\n        \n    def get_learn_rate(self):\n        \"\"\"Compute the current learning rate.\"\"\"\n        if False:\n            # Fixed learnrate\n            learn_rate = self.learn_rate_0\n        else:\n            # Decreasing learnrate each step by factor 1-alpha\n            learn_rate = self.learn_rate_0*(1.-self.learn_rate_alpha)**self.learn_rate_pos\n        return learn_rate\n\n    def next_mini_batch(self):\n        \"\"\"Get the next mini batch.\"\"\"\n        start = self.index_in_epoch\n        self.index_in_epoch += self.mb_size           \n        self.epoch += self.mb_size/len(self.x_train)\n        \n        # At the start of the epoch.\n        if start == 0:\n            np.random.shuffle(self.perm_array) # Shuffle permutation array.\n   \n        # In case the current index is larger than one epoch.\n        if self.index_in_epoch > len(self.x_train):\n            self.index_in_epoch = 0\n            self.epoch -= self.mb_size/len(self.x_train) \n            return self.next_mini_batch() # Recursive use of function.\n        \n        end = self.index_in_epoch\n        \n        # Original data.\n        x_tr = self.x_train[self.perm_array[start:end]]\n        y_tr = self.y_train[self.perm_array[start:end]]\n        \n        # Use augmented data.\n        if self.train_on_augmented_data:\n            x_tr, y_tr = generate_images_and_masks(x_tr, y_tr)\n            y_tr = trsf_proba_to_binary(y_tr)\n        \n        return x_tr, y_tr\n \n    def weight_variable(self, shape, name=None):\n        \"\"\" Weight initialization \"\"\"\n        #initializer = tf.truncated_normal(shape, stddev=0.1)\n        initializer = tf.contrib.layers.xavier_initializer()\n        #initializer = tf.contrib.layers.variance_scaling_initializer()\n        return tf.get_variable(name, shape=shape, initializer=initializer)\n\n    def bias_variable(self, shape, name=None):\n        \"\"\"Bias initialization.\"\"\"\n        #initializer = tf.constant(0.1, shape=shape)  \n        initializer = tf.contrib.layers.xavier_initializer()\n        #initializer = tf.contrib.layers.variance_scaling_initializer()\n        return tf.get_variable(name, shape=shape, initializer=initializer)\n     \n    def conv2d(self, x, W, name=None):\n        \"\"\" 2D convolution. \"\"\"\n        return tf.nn.conv2d(x, W, strides=[1,1,1,1], padding='SAME', name=name)\n\n    def max_pool_2x2(self, x, name=None):\n        \"\"\" Max Pooling 2x2. \"\"\"\n        return tf.nn.max_pool(x, ksize=[1,2,2,1], strides=[1,2,2,1], padding='SAME',\n                              name=name)\n    \n    def conv2d_transpose(self, x, filters, name=None):\n        \"\"\" Transposed 2d convolution. \"\"\"\n        return tf.layers.conv2d_transpose(x, filters=filters, kernel_size=2, \n                                          strides=2, padding='SAME') \n    \n    def leaky_relu(self, z, name=None):\n        \"\"\"Leaky ReLU.\"\"\"\n        return tf.maximum(0.01 * z, z, name=name)\n    \n    def activation(self, x, name=None):\n        \"\"\" Activation function. \"\"\"\n        a = tf.nn.elu(x, name=name)\n        #a = self.leaky_relu(x, name=name)\n        #a = tf.nn.relu(x, name=name)\n        return a \n    \n    def loss_tensor(self):\n        \"\"\"Loss tensor.\"\"\"\n        if True:\n            # Dice loss based on Jaccard dice score coefficent.\n            axis=np.arange(1,len(self.output_shape)+1)\n            offset = 1e-5\n            corr = tf.reduce_sum(self.y_data_tf * self.y_pred_tf, axis=axis)\n            l2_pred = tf.reduce_sum(tf.square(self.y_pred_tf), axis=axis)\n            l2_true = tf.reduce_sum(tf.square(self.y_data_tf), axis=axis)\n            dice_coeff = (2. * corr + 1e-5) / (l2_true + l2_pred + 1e-5)\n            # Second version: 2-class variant of dice loss\n            #corr_inv = tf.reduce_sum((1.-self.y_data_tf) * (1.-self.y_pred_tf), axis=axis)\n            #l2_pred_inv = tf.reduce_sum(tf.square(1.-self.y_pred_tf), axis=axis)\n            #l2_true_inv = tf.reduce_sum(tf.square(1.-self.y_data_tf), axis=axis)\n            #dice_coeff = ((corr + offset) / (l2_true + l2_pred + offset) +\n            #             (corr_inv + offset) / (l2_pred_inv + l2_true_inv + offset))\n            loss = tf.subtract(1., tf.reduce_mean(dice_coeff))\n        if False:\n            # Sigmoid cross entropy. \n            loss = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(\n                    labels=self.y_data_tf, logits=self.z_pred_tf))\n        return loss \n    \n    def optimizer_tensor(self):\n        \"\"\"Optimization tensor.\"\"\"\n        # Adam Optimizer (adaptive moment estimation). \n        optimizer = tf.train.AdamOptimizer(self.learn_rate_tf).minimize(\n                    self.loss_tf, name='train_step_tf')\n        return optimizer\n   \n    def batch_norm_layer(self, x, name=None):\n        \"\"\"Batch normalization layer.\"\"\"\n        if False:\n            layer = tf.layers.batch_normalization(x, training=self.training_tf, \n                                                  momentum=0.9, name=name)\n        else: \n            layer = x\n        return layer\n    \n    def dropout_layer(self, x, name=None):\n        \"\"\"Dropout layer.\"\"\"\n        if False:\n            layer = tf.layers.dropout(x, self.dropout_proba, training=self.training_tf,\n                                     name=name)\n        else:\n            layer = x\n        return layer\n\n    def num_of_weights(self,tensors):\n        \"\"\"Compute the number of weights.\"\"\"\n        sum_=0\n        for i in range(len(tensors)):\n            m = 1\n            for j in range(len(tensors[i].shape)):\n              m *= int(tensors[i].shape[j])\n            sum_+=m\n        return sum_\n\n    def build_UNet_graph(self):\n        \"\"\" Create the UNet graph in TensorFlow. \"\"\"\n        # 1. unit \n        with tf.name_scope('1.unit'):\n            W1_1 = self.weight_variable([3,3,self.input_shape[2],16], 'W1_1')\n            b1_1 = self.bias_variable([16], 'b1_1')\n            Z1 = self.conv2d(self.x_data_tf, W1_1, 'Z1') + b1_1\n            A1 = self.activation(self.batch_norm_layer(Z1)) # (.,128,128,16)\n            A1_drop = self.dropout_layer(A1)\n            W1_2 = self.weight_variable([3,3,16,16], 'W1_2')\n            b1_2 = self.bias_variable([16], 'b1_2')\n            Z2 = self.conv2d(A1_drop, W1_2, 'Z2') + b1_2\n            A2 = self.activation(self.batch_norm_layer(Z2)) # (.,128,128,16)\n            P1 = self.max_pool_2x2(A2, 'P1') # (.,64,64,16)\n        # 2. unit \n        with tf.name_scope('2.unit'):\n            W2_1 = self.weight_variable([3,3,16,32], \"W2_1\")\n            b2_1 = self.bias_variable([32], 'b2_1')\n            Z3 = self.conv2d(P1, W2_1) + b2_1\n            A3 = self.activation(self.batch_norm_layer(Z3)) # (.,64,64,32)\n            A3_drop = self.dropout_layer(A3)\n            W2_2 = self.weight_variable([3,3,32,32], \"W2_2\")\n            b2_2 = self.bias_variable([32], 'b2_2')\n            Z4 = self.conv2d(A3_drop, W2_2) + b2_2\n            A4 = self.activation(self.batch_norm_layer(Z4)) # (.,64,64,32)\n            P2 = self.max_pool_2x2(A4) # (.,32,32,32)\n        # 3. unit\n        with tf.name_scope('3.unit'):\n            W3_1 = self.weight_variable([3,3,32,64], \"W3_1\")\n            b3_1 = self.bias_variable([64], 'b3_1')\n            Z5 = self.conv2d(P2, W3_1) + b3_1\n            A5 = self.activation(self.batch_norm_layer(Z5)) # (.,32,32,64)\n            A5_drop = self.dropout_layer(A5)\n            W3_2 = self.weight_variable([3,3,64,64], \"W3_2\")\n            b3_2 = self.bias_variable([64], 'b3_2')\n            Z6 = self.conv2d(A5_drop, W3_2) + b3_2\n            A6 = self.activation(self.batch_norm_layer(Z6)) # (.,32,32,64)\n            P3 = self.max_pool_2x2(A6) # (.,16,16,64)\n        # 4. unit\n        with tf.name_scope('4.unit'):\n            W4_1 = self.weight_variable([3,3,64,128], \"W4_1\")\n            b4_1 = self.bias_variable([128], 'b4_1')\n            Z7 = self.conv2d(P3, W4_1) + b4_1\n            A7 = self.activation(self.batch_norm_layer(Z7)) # (.,16,16,128)\n            A7_drop = self.dropout_layer(A7)\n            W4_2 = self.weight_variable([3,3,128,128], \"W4_2\")\n            b4_2 = self.bias_variable([128], 'b4_2')\n            Z8 = self.conv2d(A7_drop, W4_2) + b4_2\n            A8 = self.activation(self.batch_norm_layer(Z8)) # (.,16,16,128)\n            P4 = self.max_pool_2x2(A8) # (.,8,8,128)\n        # 5. unit \n        with tf.name_scope('5.unit'):\n            W5_1 = self.weight_variable([3,3,128,256], \"W5_1\")\n            b5_1 = self.bias_variable([256], 'b5_1')\n            Z9 = self.conv2d(P4, W5_1) + b5_1\n            A9 = self.activation(self.batch_norm_layer(Z9)) # (.,8,8,256)\n            A9_drop = self.dropout_layer(A9)\n            W5_2 = self.weight_variable([3,3,256,256], \"W5_2\")\n            b5_2 = self.bias_variable([256], 'b5_2')\n            Z10 = self.conv2d(A9_drop, W5_2) + b5_2\n            A10 = self.activation(self.batch_norm_layer(Z10)) # (.,8,8,256)\n        # 6. unit\n        with tf.name_scope('6.unit'):\n            W6_1 = self.weight_variable([3,3,256,128], \"W6_1\")\n            b6_1 = self.bias_variable([128], 'b6_1')\n            U1 = self.conv2d_transpose(A10, 128) # (.,16,16,128)\n            U1 = tf.concat([U1, A8], 3) # (.,16,16,256)\n            Z11 = self.conv2d(U1, W6_1) + b6_1\n            A11 = self.activation(self.batch_norm_layer(Z11)) # (.,16,16,128)\n            A11_drop = self.dropout_layer(A11)\n            W6_2 = self.weight_variable([3,3,128,128], \"W6_2\")\n            b6_2 = self.bias_variable([128], 'b6_2')\n            Z12 = self.conv2d(A11_drop, W6_2) + b6_2\n            A12 = self.activation(self.batch_norm_layer(Z12)) # (.,16,16,128)\n        # 7. unit \n        with tf.name_scope('7.unit'):\n            W7_1 = self.weight_variable([3,3,128,64], \"W7_1\")\n            b7_1 = self.bias_variable([64], 'b7_1')\n            U2 = self.conv2d_transpose(A12, 64) # (.,32,32,64)\n            U2 = tf.concat([U2, A6],3) # (.,32,32,128)\n            Z13 = self.conv2d(U2, W7_1) + b7_1\n            A13 = self.activation(self.batch_norm_layer(Z13)) # (.,32,32,64)\n            A13_drop = self.dropout_layer(A13)\n            W7_2 = self.weight_variable([3,3,64,64], \"W7_2\")\n            b7_2 = self.bias_variable([64], 'b7_2')\n            Z14 = self.conv2d(A13_drop, W7_2) + b7_2\n            A14 = self.activation(self.batch_norm_layer(Z14)) # (.,32,32,64)\n        # 8. unit\n        with tf.name_scope('8.unit'):\n            W8_1 = self.weight_variable([3,3,64,32], \"W8_1\")\n            b8_1 = self.bias_variable([32], 'b8_1')\n            U3 = self.conv2d_transpose(A14, 32) # (.,64,64,32)\n            U3 = tf.concat([U3, A4],3) # (.,64,64,64)\n            Z15 = self.conv2d(U3, W8_1) + b8_1\n            A15 = self.activation(self.batch_norm_layer(Z15)) # (.,64,64,32)\n            A15_drop = self.dropout_layer(A15)\n            W8_2 = self.weight_variable([3,3,32,32], \"W8_2\")\n            b8_2 = self.bias_variable([32], 'b8_2')\n            Z16 = self.conv2d(A15_drop, W8_2) + b8_2\n            A16 = self.activation(self.batch_norm_layer(Z16)) # (.,64,64,32)\n        # 9. unit \n        with tf.name_scope('9.unit'):\n            W9_1 = self.weight_variable([3,3,32,16], \"W9_1\")\n            b9_1 = self.bias_variable([16], 'b9_1')\n            U4 = self.conv2d_transpose(A16, 16) # (.,128,128,16)\n            U4 = tf.concat([U4, A2],3) # (.,128,128,32)\n            Z17 = self.conv2d(U4, W9_1) + b9_1\n            A17 = self.activation(self.batch_norm_layer(Z17)) # (.,128,128,16)\n            A17_drop = self.dropout_layer(A17)\n            W9_2 = self.weight_variable([3,3,16,16], \"W9_2\")\n            b9_2 = self.bias_variable([16], 'b9_2')\n            Z18 = self.conv2d(A17_drop, W9_2) + b9_2\n            A18 = self.activation(self.batch_norm_layer(Z18)) # (.,128,128,16)\n        # 10. unit: output layer\n        with tf.name_scope('10.unit'):\n            W10 = self.weight_variable([1,1,16,1], \"W10\")\n            b10 = self.bias_variable([1], 'b10')\n            Z19 = self.conv2d(A18, W10) + b10\n            A19 = tf.nn.sigmoid(self.batch_norm_layer(Z19)) # (.,128,128,1)\n        \n        self.z_pred_tf = tf.identity(Z19, name='z_pred_tf') # (.,128,128,1)\n        self.y_pred_tf = tf.identity(A19, name='y_pred_tf') # (.,128,128,1)\n        \n        print('Build UNet Graph: 10 layers, {} trainable weights'.format(\n            self.num_of_weights([W1_1,b1_1,W1_2,b1_2,W2_1,b2_1,W2_2,b2_2,\n                                 W3_1,b3_1,W3_2,b3_2,W4_1,b4_1,W4_2,b4_2,\n                                 W5_1,b5_1,W5_2,b5_2,W6_1,b6_1,W6_2,b6_2,\n                                 W7_1,b7_1,W7_2,b7_2,W8_1,b8_1,W8_2,b8_2,\n                                 W9_1,b9_1,W9_2,b9_2,W10,b10])))\n    \n    def build_graph(self):\n        \"\"\" Build the complete graph in TensorFlow. \"\"\"\n        tf.reset_default_graph()  \n        self.graph = tf.Graph()\n\n        with self.graph.as_default():\n            \n            # Input tensor.\n            shape = [None]\n            shape = shape.extend(self.input_shape)\n            self.x_data_tf = tf.placeholder(dtype=tf.float32, shape=shape, \n                                            name='x_data_tf') # (.,128,128,3)\n            \n            # Generic tensors.\n            self.keep_prob_tf = tf.placeholder_with_default(1.0, shape=(), \n                                                            name='keep_prob_tf') \n            self.learn_rate_tf = tf.placeholder(dtype=tf.float32,\n                                                name=\"learn_rate_tf\")\n            self.training_tf = tf.placeholder_with_default(False, shape=(),\n                                                           name='training_tf')\n            # Build U-Net graph.\n            self.build_UNet_graph()\n\n            # Target tensor.\n            shape = [None]\n            shape = shape.extend(self.output_shape)\n            self.y_data_tf = tf.placeholder(dtype=tf.float32, shape=shape, \n                                            name='y_data_tf') # (.,128,128,1)\n            # Loss tensor\n            self.loss_tf = tf.identity(self.loss_tensor(), name='loss_tf')\n\n            # Optimisation tensor.\n            self.train_step_tf = self.optimizer_tensor()\n            \n            # Extra operations required for batch normalization.\n            self.extra_update_ops_tf = tf.get_collection(tf.GraphKeys.UPDATE_OPS)\n\n            \n    def train_graph(self, sess, x_train, y_train, x_valid, y_valid, n_epoch=1, \n                    train_on_augmented_data=False):\n        \"\"\" Train the graph of the corresponding neural network. \"\"\"\n        # Set training and validation sets.\n        self.x_train = x_train\n        self.y_train = y_train\n        self.x_valid = x_valid\n        self.y_valid = y_valid\n        \n        # Parameters.\n        self.perm_array = np.arange(len(self.x_train))\n        self.train_on_augmented_data = train_on_augmented_data\n        mb_per_epoch = self.x_train.shape[0]/self.mb_size\n        train_loss, train_score, valid_loss, valid_score = [],[],[],[]\n        \n        # Start timer.\n        start = datetime.datetime.now()\n        print('Training the Neural Network')\n        print('\\tnn_name = {}, n_epoch = {}, mb_size = {}, learnrate = {:.7f}'.format(\n               self.nn_name, n_epoch, self.mb_size, self.learn_rate))\n        print('\\tinput_shape = {}, output_shape = {}'.format(\n            self.input_shape, self.output_shape))\n        print('\\tlearn_rate = {:.10f}, learn_rate_0 = {:.10f}, learn_rate_alpha = {}'.format(\n            self.learn_rate, self.learn_rate_0, self.learn_rate_alpha))\n        print('\\tlearn_rate_step = {}, learn_rate_pos = {}, dropout_proba = {}'.format(\n            self.learn_rate_step, self.learn_rate_pos, self.dropout_proba))\n        print('\\tx_train = {}, x_valid = {}'.format(x_train.shape, x_valid.shape))\n        print('\\ty_train = {}, y_valid = {}'.format(y_train.shape, y_valid.shape))\n        print('Training started: {}'.format(datetime.datetime.now().strftime(\n                                     '%d-%m-%Y %H:%M:%S')))\n        \n        # Looping over mini batches.\n        for i in range(int(n_epoch*mb_per_epoch)+1):\n\n            # Adapt the learning rate.\n            if not self.learn_rate_pos == int(self.epoch // self.learn_rate_step):\n                self.learn_rate_pos = int(self.epoch // self.learn_rate_step)\n                self.learn_rate = self.get_learn_rate()\n                print('Update learning rate to {:.10f}. Running time: {}'.format(\n                    self.learn_rate, datetime.datetime.now()-start))\n            \n            # Train the graph.\n            x_batch, y_batch = self.next_mini_batch() # next mini batch\n            sess.run([self.train_step_tf, self.extra_update_ops_tf], \n                     feed_dict={self.x_data_tf: x_batch, self.y_data_tf: y_batch, \n                                self.keep_prob_tf: self.keep_prob, \n                                self.learn_rate_tf: self.learn_rate,\n                                self.training_tf: True})\n            \n            # Store losses and scores.\n            if i%int(self.log_step*mb_per_epoch) == 0:\n             \n                self.n_log_step += 1 # Current number of log steps.\n                \n                # Train data used for evaluation.\n                ids = np.arange(len(self.x_train))\n                np.random.shuffle(ids)\n                ids = ids[:len(x_valid)] # len(x_batch)\n                x_trn = self.x_train[ids]\n                y_trn = self.y_train[ids]\n                \n                # Valid data used for evaluation.\n                ids = np.arange(len(self.x_valid))\n                np.random.shuffle(ids)\n                ids = ids[:len(x_valid)] # len(x_batch)\n                x_vld = self.x_valid[ids]\n                y_vld = self.y_valid[ids]\n                \n                feed_dict_train = {self.x_data_tf: x_trn, self.y_data_tf: y_trn, \n                                   self.keep_prob_tf: 1.0}\n                feed_dict_valid = {self.x_data_tf: x_vld, self.y_data_tf: y_vld, \n                                   self.keep_prob_tf: 1.0}\n                \n                # Evaluate current loss and score\n                train_loss, y_train_pred = sess.run([self.loss_tf, self.y_pred_tf], \n                                                   feed_dict = feed_dict_train)\n                valid_loss, y_valid_pred = sess.run([self.loss_tf, self.y_pred_tf], \n                                                   feed_dict = feed_dict_valid)\n                train_score = np.mean(get_score(y_trn, y_train_pred))\n                valid_score = np.mean(get_score(y_vld, y_valid_pred))\n                \n                print(('{:.2f} epoch: train/valid loss = {:.4f}/{:.4f} ' + \n                       'train/valid score = {:.4f}/{:.4f}').format(\n                        self.epoch, train_loss, valid_loss,  train_score, valid_score))\n\n                # Store losses and scores.\n                self.params['train_loss'].extend([train_loss])\n                self.params['valid_loss'].extend([valid_loss])\n                self.params['train_score'].extend([train_score])\n                self.params['valid_score'].extend([valid_score])\n        \n                # Save summaries for TensorBoard.\n                if self.use_tb_summary:\n                    train_summary = sess.run(self.merged, feed_dict = feed_dict_train)\n                    valid_summary = sess.run(self.merged, feed_dict = feed_dict_valid)\n                    self.train_writer.add_summary(train_summary, self.n_log_step)\n                    self.valid_writer.add_summary(valid_summary, self.n_log_step)\n                \n        # Store parameters.\n        self.params['learn_rate'] = self.learn_rate\n        self.params['learn_rate_step'] = self.learn_rate_step\n        self.params['learn_rate_pos'] = self.learn_rate_pos\n        self.params['learn_rate_alpha'] = self.learn_rate_alpha\n        self.params['learn_rate_0'] = self.learn_rate_0\n        self.params['keep_prob'] = self.keep_prob\n        self.params['epoch'] = self.epoch\n        self.params['n_log_step'] = self.n_log_step\n        self.params['log_step'] = self.log_step\n        self.params['input_shape'] = self.input_shape\n        self.params['output_shape'] = self.output_shape\n        self.params['mb_size'] = self.mb_size\n        self.params['dropout_proba'] = self.dropout_proba\n        \n        print('Training ended. Running time: {}'.format(datetime.datetime.now()-start))\n    \n    def summary_variable(self, var, var_name):\n        \"\"\" Attach summaries to a tensor for TensorBoard visualization. \"\"\"\n        with tf.name_scope(var_name):\n            mean = tf.reduce_mean(var)\n            stddev = tf.sqrt(tf.reduce_mean(tf.square(var - mean)))\n            tf.summary.scalar('mean', mean)\n            tf.summary.scalar('stddev', stddev)\n            tf.summary.scalar('max', tf.reduce_max(var))\n            tf.summary.scalar('min', tf.reduce_min(var))\n            tf.summary.histogram('histogram', var)\n\n    def attach_summary(self, sess):\n        \"\"\" Attach TensorBoard summaries to certain tensors. \"\"\"\n        self.use_tb_summary = True\n        \n        # Create summary tensors for TensorBoard.\n        tf.summary.scalar('loss_tf', self.loss_tf)\n\n        # Merge all summaries.\n        self.merged = tf.summary.merge_all()\n\n        # Initialize summary writer.\n        timestamp = datetime.datetime.now().strftime('%d-%m-%Y_%H-%M-%S')\n        filepath = os.path.join(os.getcwd(), LOGS_DIR_NAME, (self.nn_name+'_'+timestamp))\n        self.train_writer = tf.summary.FileWriter(os.path.join(filepath,'train'), sess.graph)\n        self.valid_writer = tf.summary.FileWriter(os.path.join(filepath,'valid'), sess.graph)\n\n    def attach_saver(self):\n        \"\"\" Initialize TensorFlow saver. \"\"\"\n        with self.graph.as_default():\n            self.use_tf_saver = True\n            self.saver_tf = tf.train.Saver()\n\n    def save_model(self, sess):\n        \"\"\" Save parameters, tensors and summaries. \"\"\"\n        if not os.path.isdir(os.path.join(CW_DIR, SAVES_DIR_NAME)):\n            os.mkdir(SAVES_DIR_NAME)\n        filepath = os.path.join(os.getcwd(), SAVES_DIR_NAME , self.nn_name+'_params.npy')\n        np.save(filepath, self.params) # save parameters of the network\n\n        # TensorFlow saver\n        if self.use_tf_saver:\n            filepath = os.path.join(os.getcwd(),  self.nn_name)\n            self.saver_tf.save(sess, filepath)\n\n        # TensorBoard summaries\n        if self.use_tb_summary:\n            self.train_writer.close()\n            self.valid_writer.close()\n        \n    def load_session_from_file(self, filename):\n        \"\"\" Load session from a file, restore the graph, and load the tensors. \"\"\"\n        tf.reset_default_graph()\n        filepath = os.path.join(os.getcwd(), filename + '.meta')\n        saver = tf.train.import_meta_graph(filepath)\n        sess = tf.Session() # default session\n        saver.restore(sess, filename) # restore session\n        self.graph = tf.get_default_graph() # save default graph\n        self.load_parameters(filename) # load parameters\n        self.load_tensors(self.graph) # define relevant tensors as variables \n        return sess\n    \n    def load_parameters(self, filename):\n        '''Load helper and tunable parameters.'''\n        filepath = os.path.join(os.getcwd(), SAVES_DIR_NAME, filename+'_params.npy')\n        self.params = np.load(filepath).item() # load parameters of network\n        \n        self.nn_name = filename\n        self.learn_rate = self.params['learn_rate']\n        self.learn_rate_0 = self.params['learn_rate_0']\n        self.learn_rate_step = self.params['learn_rate_step']\n        self.learn_rate_alpha = self.params['learn_rate_alpha']\n        self.learn_rate_pos = self.params['learn_rate_pos']\n        self.keep_prob = self.params['keep_prob']\n        self.epoch = self.params['epoch'] \n        self.n_log_step = self.params['n_log_step']\n        self.log_step = self.params['log_step']\n        self.input_shape = self.params['input_shape']\n        self.output_shape = self.params['output_shape'] \n        self.mb_size = self.params['mb_size']   \n        self.dropout_proba = self.params['dropout_proba']\n        \n        print('Parameters of the loaded neural network')\n        print('\\tnn_name = {}, epoch = {:.2f}, mb_size = {}'.format(\n            self.nn_name, self.epoch, self.mb_size))\n        print('\\tinput_shape = {}, output_shape = {}'.format(\n            self.input_shape, self.output_shape))\n        print('\\tlearn_rate = {:.10f}, learn_rate_0 = {:.10f}, dropout_proba = {}'.format(\n            self.learn_rate, self.learn_rate_0, self.dropout_proba))\n        print('\\tlearn_rate_step = {}, learn_rate_pos = {}, learn_rate_alpha = {}'.format(\n            self.learn_rate_step, self.learn_rate_pos, self.learn_rate_alpha))\n\n    def load_tensors(self, graph):\n        \"\"\" Load tensors from a graph. \"\"\"\n        # Input tensors\n        self.x_data_tf = graph.get_tensor_by_name(\"x_data_tf:0\")\n        self.y_data_tf = graph.get_tensor_by_name(\"y_data_tf:0\")\n\n        # Tensors for training and prediction.\n        self.learn_rate_tf = graph.get_tensor_by_name(\"learn_rate_tf:0\")\n        self.keep_prob_tf = graph.get_tensor_by_name(\"keep_prob_tf:0\")\n        self.loss_tf = graph.get_tensor_by_name('loss_tf:0')\n        self.train_step_tf = graph.get_operation_by_name('train_step_tf')\n        self.z_pred_tf = graph.get_tensor_by_name('z_pred_tf:0')\n        self.y_pred_tf = graph.get_tensor_by_name(\"y_pred_tf:0\")\n        self.training_tf = graph.get_tensor_by_name(\"training_tf:0\")\n        self.extra_update_ops_tf = tf.get_collection(tf.GraphKeys.UPDATE_OPS)\n\n    def get_prediction(self, sess, x_data, keep_prob=1.0):\n        \"\"\" Prediction of the neural network graph. \"\"\"\n        return sess.run(self.y_pred_tf, feed_dict={self.x_data_tf: x_data,\n                                                     self.keep_prob_tf: keep_prob})\n       \n    def get_loss(self, sess, x_data, y_data, keep_prob=1.0):\n        \"\"\" Compute the loss. \"\"\"\n        return sess.run(self.loss_tf, feed_dict={self.x_data_tf: x_data, \n                                                 self.y_data_tf: y_data,\n                                                 self.keep_prob_tf: keep_prob})\n ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1e7847e8b5e85af1610da427b3935794b49060e","_cell_guid":"43920015-f35f-4546-929f-208acfb525a0"},"cell_type":"markdown","source":"# 6. Train the Neural Network <a class=\"anchor\" id=\"6-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"4c250877056775e86b2ff91fd6bc2f814b7037a5","collapsed":true,"_cell_guid":"7a7a3deb-702b-476e-8e92-16e0e39c9025","trusted":false},"cell_type":"code","source":"# In case you want to reload and preprocess the raw data.\nif False:\n    x_train, y_train, x_test = load_raw_data((128,128))\n    x_train, y_train, x_test = preprocess_raw_data(x_train, y_train, x_test, grayscale=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8cd4261d15154345cc0bb5e77d5ac9a1d3f176f3","_cell_guid":"b1269b86-40b1-44b2-9564-f3d141133e40","scrolled":true,"trusted":false,"collapsed":true},"cell_type":"code","source":"# Create and start training of a new neural network, or continue training of\n# a pretrained model.\n\n#nn_name = ['nn0','nn1','nn2','nn3','nn4','nn5','nn6','nn7','nn8','nn9']\n#nn_name = ['nn0_128_128_3']  \n#nn_name = ['nn0_384_384_3']  \n#nn_name = ['nn0_512_512_3']  \nnn_name = ['nn0_256_256_3']  \n\n# Implement cross validations\ncv_num = 10 \nkfold = sklearn.model_selection.KFold(cv_num, shuffle=True, random_state=SEED)\n\nfor i,(train_index, valid_index) in enumerate(kfold.split(x_train)):\n\n    # Start timer\n    start = datetime.datetime.now();\n\n    # Split into train and validation \n    x_trn = x_train[train_index]\n    y_trn = y_train[train_index]\n    x_vld = x_train[valid_index]\n    y_vld = y_train[valid_index]\n    \n    # Choose a certain fold.\n    if i==0:\n        # Create and start training of a new model.\n        if True:\n            u_net = NeuralNetwork(nn_name=nn_name[0], log_step=1.0, input_shape=(256,256,3), \n                                  output_shape=(256,256,1)) # Create instance of neural network.\n            u_net.build_graph() # Build graph.\n\n            # Start tensorflow session.\n            with tf.Session(graph=u_net.graph) as sess: \n                u_net.attach_saver() # Attach saver tensor.\n                u_net.attach_summary(sess) # Attach summaries.\n                sess.run(tf.global_variables_initializer()) # Variable initialization.\n\n                # Training on original data.\n                u_net.train_graph(sess, x_trn, y_trn, x_vld, y_vld, n_epoch=1.)\n\n                for _ in range(1):\n                    # Training on augmented data.\n                    u_net.train_graph(sess, x_trn, y_trn, x_vld, y_vld, n_epoch=3.,\n                                      train_on_augmented_data=True)\n                    u_net.save_model(sess) # Save parameters, tensors, summaries.\n\n        # Continue training of a pretrained model.\n        if False:\n            u_net = NeuralNetwork() \n            sess = u_net.load_session_from_file(nn_name[0])  \n            u_net.attach_saver() \n            u_net.attach_summary(sess) \n\n            # Training on original data.\n            #u_net.train_graph(sess, x_trn, y_trn, x_vld, y_vld, n_epoch=0.1) \n\n            for _ in range(9):\n                # Training on augmented data.\n                u_net.train_graph(sess, x_trn, y_trn, x_vld, y_vld, n_epoch=5.,\n                                  train_on_augmented_data = True)\n                u_net.save_model(sess) # Save parameters, tensors, summaries.\n\nprint('Total running time: ', datetime.datetime.now() - start)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"288f1e68f32103babc19febfbcf8ef8249d216c0","collapsed":true,"_cell_guid":"31bb8903-3270-49d3-9034-15f44e6b9df1","trusted":false},"cell_type":"code","source":"# Start TensorBoard visualization. All summaries are written into the\n# logs directory which is contained in the current working directory.\nif False:\n    !tensorboard --logdir=./logs","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f9d02a878056ae3dac6de883d3e44195564c2a9d","_cell_guid":"90181849-927e-4e8d-92a5-4d1dc7b24c1d","trusted":false,"collapsed":true},"cell_type":"code","source":"# Show intermediate losses and scores during the training session.\n#mn = 'nn0_512_512_3' \nmn = 'nn0_256_256_3' \nu_net = NeuralNetwork()\nsess = u_net.load_session_from_file(mn)\nsess.close()\ntrain_loss = u_net.params['train_loss']\nvalid_loss = u_net.params['valid_loss']\ntrain_score = u_net.params['train_score']\nvalid_score = u_net.params['valid_score']\n\nprint('final train/valid loss = {:.4f}/{:.4f}'.format(train_loss[-1], valid_loss[-1]))\nprint('final train/valid score = {:.4f}/{:.4f}'.format(train_score[-1], valid_score[-1]))\nplt.figure(figsize=(10, 5));\nplt.subplot(1,2,1)\nplt.plot(np.arange(0,len(train_loss)), train_loss,'-b', label='Training')\nplt.plot(np.arange(0,len(valid_loss)), valid_loss,'-g', label='Validation')\nplt.legend(loc='lower right', frameon=False)\nplt.ylim(ymax = 0.5, ymin = 0.0)\nplt.ylabel('loss')\nplt.xlabel('steps');\n\nplt.subplot(1,2,2);\nplt.plot(np.arange(0,len(train_score)), train_score,'-b', label='Training')\nplt.plot(np.arange(0,len(valid_score)), valid_score,'-g', label='Validation')\nplt.legend(loc='lower right', frameon=False)\nplt.ylim(ymax = 0.7, ymin = 0.0)\nplt.ylabel('score')\nplt.xlabel('steps');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0b4cba240ed99085d6adf2b33b181d6578b13c95","_cell_guid":"75e26192-41b4-4a32-b108-1a20eecd95d9"},"cell_type":"markdown","source":"# 7. Validate the neural network <a class=\"anchor\" id=\"7-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"719e5fee15550c38283287e91ca4a0acc1dbe2f3","collapsed":true,"_cell_guid":"3c3d59a6-85aa-44e5-b93d-06c5e981b1b8","trusted":false},"cell_type":"code","source":"# Split training and validation data in the same way \n# as was done for training the neural network.\nif True:\n    cv_num = 10 \n    kfold = sklearn.model_selection.KFold(cv_num, shuffle=True, random_state=SEED)\n    for i,(train_index, valid_index) in enumerate(kfold.split(x_train)):\n        if i==1:\n            x_trn = x_train[train_index]\n            y_trn = y_train[train_index]\n            x_vld = x_train[valid_index]\n            y_vld = y_train[valid_index]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b4a1a36a8b832e922e8638513ef447814b85f75c","_cell_guid":"953ae4df-b0b7-4d67-89de-ff0bf6226773","trusted":false,"collapsed":true},"cell_type":"code","source":"# Summary of scores for training and validations sets. Note that the score is\n# better than the true score, since overlapping/touching nuclei can not be\n# separately identified in this version.  \n#mn = 'nn0_512_512_3'\nmn = 'nn0_256_256_3'\nu_net = NeuralNetwork()\nsess = u_net.load_session_from_file(mn)\n\nif False:\n    # Overall score on train set.\n    y_train_pred = trsf_proba_to_binary(u_net.get_prediction(sess, x_trn))\n    train_score = get_score(y_trn, y_train_pred)\n    tmp = np.concatenate([train_index.reshape(-1,1), train_score.reshape(-1,1)], axis=1)\n    train_score_df = pd.DataFrame(tmp, columns=(['train_index','train_score']))\n    print(train_score_df.describe())\n    print('')\n    print(train_score_df.sort_values(by='train_score', ascending=True).head())\n\nif True:\n    # Overall score on validation set.\n    y_valid_pred_proba = u_net.get_prediction(sess, x_vld)\n    y_valid_pred = trsf_proba_to_binary(y_valid_pred_proba)\n    valid_score = get_score(y_vld, y_valid_pred)\n    tmp = np.concatenate([np.arange(len(valid_index)).reshape(-1,1),\n                          valid_index.reshape(-1,1),\n                          valid_score.reshape(-1,1)], axis=1)\n    valid_score_df = pd.DataFrame(tmp, columns=(['index', 'valid_index','valid_score']))\n    print('\\n', valid_score_df.describe())\n    print('\\n', valid_score_df.sort_values(by='valid_score', ascending=True).head())\n\n    # Plot the worst 4 predictions.\n    fig, axs = plt.subplots(4,4,figsize=(20,20))\n    list_ = valid_score_df.sort_values(by='valid_score', ascending=True)[:4]['index'].values.astype(np.int)\n    #list_ = [valid_score_df['valid_score'].idxmin(),valid_score_df['valid_score'].idxmax()]\n    for i,n in enumerate(list_):\n        img, img_type = imshow_args(x_vld[n])\n        axs[i,0].imshow(img, img_type)\n        axs[i,0].set_title('{}.) input image'.format(n))\n        axs[i,1].imshow(y_vld[n,:,:,0], cm.gray)\n        axs[i,1].set_title('{}.) true mask'.format(n))\n        axs[i,2].imshow(y_valid_pred_proba[n,:,:,0], cmap='gray') \n        axs[i,2].set_title('{}.) predicted mask probabilities'.format(n))\n        axs[i,3].imshow(y_valid_pred[n,:,:,0], cmap='gray') \n        axs[i,3].set_title('{}.) predicted mask'.format(n));\n\nsess.close()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5754747c8e470447402032722095e28f74e1e062","scrolled":true,"_cell_guid":"2956f19c-9478-4c91-96c3-2c477a9cd68a","trusted":false,"collapsed":true},"cell_type":"code","source":"# Tune minimal object size for prediction \nif True:\n    #mn = 'nn0_512_512_3'\n    mn = 'nn0_256_256_3'\n    u_net = NeuralNetwork()\n    sess = u_net.load_session_from_file(mn)\n    y_valid_pred_proba = u_net.get_prediction(sess, x_vld)\n    y_valid_pred = trsf_proba_to_binary(y_valid_pred_proba)\n    sess.close()        \n    \n    tmp = min_object_size\n    min_object_sizes = [1,3,5,7,9,20,30,40,50,60,70,80,90,100,110,120,130,140,150,200,300,400,500]\n    for mos in min_object_sizes:\n        min_object_size = mos\n        valid_score = get_score(y_vld, y_valid_pred)\n        print('min_object_size = {}: valid_score min/mean/std/max = {:.3f}/{:.3f}/{:.3f}/{:.3f}'.format(mos, \n                np.min(valid_score), np.mean(valid_score), np.std(valid_score), np.max(valid_score)))\n    min_object_size = tmp","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1965108dfa981b7e8b726c2a13655d1cfcf8120d","scrolled":true,"_cell_guid":"b0510aad-86da-4d8e-a4e9-510fd8196a14","trusted":false,"collapsed":true},"cell_type":"code","source":"# Check one sample prediction in more detail.\n#mn = 'nn0_512_512_3'\nmn = 'nn0_256_256_3'\nu_net = NeuralNetwork()\nsess = u_net.load_session_from_file(mn)\nn = 40 # np.random.randint(len(x_train))\nx_true = x_vld[n]\ny_true = y_vld[n,:,:,0]\ny_pred_proba = (u_net.get_prediction(sess, np.expand_dims(x_true, axis=0)))[0,:,:,0]\ny_pred = trsf_proba_to_binary(y_pred_proba)\nsess.close()\n\nfig, axs = plt.subplots(1,3,figsize=(20,13))\nimg, img_type = imshow_args(x_true)\naxs[0].imshow(img, img_type)\naxs[0].set_title('{}.) input image'.format(n))\naxs[1].imshow(y_pred_proba, cmap='gray') \naxs[1].set_title('{}.) predicted mask probabilities'.format(n));\naxs[2].imshow(y_pred, cmap='gray') \naxs[2].set_title('{}.) predicted mask'.format(n));\nplot_score_summary(y_true, y_pred)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9cd3d6fbd034498973c1928cefe589c5ba6bc8c4","_cell_guid":"93f27a93-6edf-4849-8a0e-a087e5b35b50"},"cell_type":"markdown","source":"# 8. Make Test Prediction <a class=\"anchor\" id=\"8-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"5f2b136b26be120d9020ccd750b374d4a7c1c141","collapsed":true,"_cell_guid":"6714b081-b4fc-400c-9ad1-686b2904952b","trusted":false},"cell_type":"code","source":"# Collection of methods for run length encoding. \n# For example, '1 3 10 5' implies pixels 1,2,3,10,11,12,13,14 are to be included \n# in the mask. The pixels are one-indexed and numbered from top to bottom, \n# then left to right: 1 is pixel (1,1), 2 is pixel (2,1), etc.\n\ndef rle_of_binary(x):\n    \"\"\" Run length encoding of a binary 2D array. \"\"\"\n    dots = np.where(x.T.flatten() == 1)[0] # indices from top to down\n    run_lengths = []\n    prev = -2\n    for b in dots:\n        if (b>prev+1): run_lengths.extend((b + 1, 0))\n        run_lengths[-1] += 1\n        prev = b\n    return run_lengths\n\ndef mask_to_rle(mask, cutoff=.5, min_object_size=1.):\n    \"\"\" Return run length encoding of mask. \"\"\"\n    # segment image and label different objects\n    lab_mask = skimage.morphology.label(mask > cutoff)\n    \n    # Keep only objects that are large enough.\n    (mask_labels, mask_sizes) = np.unique(lab_mask, return_counts=True)\n    if (mask_sizes < min_object_size).any():\n        mask_labels = mask_labels[mask_sizes < min_object_size]\n        for n in mask_labels:\n            lab_mask[lab_mask == n] = 0\n        lab_mask = skimage.morphology.label(lab_mask > cutoff) \n        \n    # Loop over each object excluding the background labeled by 0.\n    for i in range(1, lab_mask.max() + 1):\n        yield rle_of_binary(lab_mask == i)\n        \ndef rle_to_mask(rle, img_shape):\n    ''' Return mask from run length encoding.'''\n    mask_rec = np.zeros(img_shape).flatten()\n    for n in range(len(rle)):\n        for i in range(0,len(rle[n]),2):\n            for j in range(rle[n][i+1]): \n                mask_rec[rle[n][i]-1+j] = 1\n    return mask_rec.reshape(img_shape[1], img_shape[0]).T\n        ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"08e540f6c97c5c6976b8e42b44704c431909c83e","_cell_guid":"7023c9d3-47d0-4c1e-8938-b405c493feb9","trusted":false,"collapsed":true},"cell_type":"code","source":"# Load neural network, make prediction for test masks, resize predicted\n# masks to original image size and apply run length encoding for the\n# submission file. \n\n# Load neural network and make prediction for masks.\n#nn_name = ['nn0_512_512_3']\nnn_name = ['nn0_256_256_3']\n\n# Soft voting majority.\nfor i,mn in enumerate(nn_name):\n    u_net = NeuralNetwork()\n    sess = u_net.load_session_from_file(mn)\n    if i==0: \n        y_test_pred_proba = u_net.get_prediction(sess, x_test)/len(nn_name)\n    else:\n        y_test_pred_proba += u_net.get_prediction(sess, x_test)/len(nn_name)\n    sess.close()\n\ny_test_pred = trsf_proba_to_binary(y_test_pred_proba)\nprint('y_test_pred.shape = {}'.format(y_test_pred.shape))\n\n# Resize predicted masks to original image size.\ny_test_pred_original_size = []\nfor i in range(len(y_test_pred)):\n    res_mask = trsf_proba_to_binary(skimage.transform.resize(np.squeeze(y_test_pred[i]),\n        (test_df.loc[i,'img_height'], test_df.loc[i,'img_width']), \n        mode='constant', preserve_range=True))\n    y_test_pred_original_size.append(res_mask)\ny_test_pred_original_size = np.array(y_test_pred_original_size)\n\nprint('y_test_pred_original_size.shape = {}'.format(y_test_pred_original_size.shape))\n   \n# Run length encoding of predicted test masks.\ntest_pred_rle = []\ntest_pred_ids = []\nfor n, id_ in enumerate(test_df['img_id']):\n    min_object_size = 20*test_df.loc[n,'img_height']*test_df.loc[n,'img_width']/(256*256)\n    rle = list(mask_to_rle(y_test_pred_original_size[n], min_object_size=min_object_size))\n    test_pred_rle.extend(rle)\n    test_pred_ids.extend([id_]*len(rle))\n\nprint('test_pred_ids.shape = {}'.format(np.array(test_pred_ids).shape))\nprint('test_pred_rle.shape = {}'.format(np.array(test_pred_rle).shape))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b8c28b4c33d755eece4c5bfb953b121124f7cec7","_cell_guid":"ed91bb1b-b453-41cc-94dd-5d104f2d5c73","trusted":false,"collapsed":true},"cell_type":"code","source":"# Inspect a test prediction and check run length encoding.\nn = 8 # np.random.randint(len(x_test))\nmask = y_test_pred_original_size[n]\nrle = list(mask_to_rle(mask))\nmask_rec = rle_to_mask(rle, mask.shape)\nprint('Run length encoding: {} matches, {} misses'.format(\n    (mask_rec == mask).sum(),(mask_rec != mask).sum()))\n\nfig, axs = plt.subplots(2,3,figsize=(20,13))\naxs[0,0].imshow(read_image(test_df['image_path'].loc[n]))\naxs[0,0].set_title('{}.) original test image'.format(n))\nimg, img_type = imshow_args(x_test[n])\naxs[0,1].imshow(img, img_type)\naxs[0,1].set_title('{}.) transformed test image'.format(n))\naxs[0,2].imshow(y_test_pred_proba[n][:,:,0], cm.gray) \naxs[0,2].set_title('{}.) predicted test mask probabilities'.format(n));\naxs[1,0].imshow(y_test_pred_proba[n][:,:,0], cm.gray) \naxs[1,0].set_title('{}.) predicted test mask'.format(n));\naxs[1,1].imshow(y_test_pred_original_size[n], cm.gray);\naxs[1,1].set_title('{}.) predicted final test mask in original size'.format(n));\naxs[1,2].imshow(mask_rec[:,:], cm.gray);\naxs[1,2].set_title('{}.) final mask recovered from run length encoding'.format(n));","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e2f9f76bd02069a7bc07e5f825e7d05dcedc8cd9","_cell_guid":"7f7f91be-3592-4ea4-a2aa-cf3c8aa5a8ef"},"cell_type":"markdown","source":"# 9. Submit <a class=\"anchor\" id=\"9-bullet\"></a> ","outputs":[],"execution_count":null},{"metadata":{"_uuid":"f51cb0fb633631b491cf29ebeff624b178518acc","_cell_guid":"5829b715-def7-4e50-bbee-3309fbc4dae2","trusted":false,"collapsed":true},"cell_type":"code","source":"# Create submission file\nsub = pd.DataFrame()\nsub['ImageId'] = test_pred_ids\nsub['EncodedPixels'] = pd.Series(test_pred_rle).apply(lambda x: ' '.join(str(y) for y in x))\nsub.to_csv('sub-dsbowl2018-1.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c29ede08bae7ae8b07186141f647ab2e038bb9b2","collapsed":true,"_cell_guid":"105994cb-8a91-4558-9ad4-3a4f842f1cbc","trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}