{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 3D ConvNet with 3D Augmentations\n\nThis is a kernel I was playing around with to see if a 3D convolutional network could produce useful features from scan images preprocessed to the same size using code as published in my previous [kernel](https://www.kaggle.com/donkeys/preprocessing-images-to-normalize-colors-and-sizes). To recap, this includes:\n\n- resizing all images to same size (width, height)\n- rescaling the z-axis of all 3D arrays to the same depth.\n- removing boundaries from images there present.\n\nSince there are not that many images to start with, I also added augmentations to the 3D images, according to my earlier [kernel](https://www.kaggle.com/donkeys/ct-slices-basic-eda):\n\n- 3D gaussian blur\n- 3D flips on x- and y-axis\n- 3D rotation on x-axis\n- 3D shift on x- and y-axis\n- 3D zoom on x- and y-axis\n\nThe 3D in the above list simply refers to applying the augmentations/transformations on the whole z-axis of the 3D numpy array at once.\n\nOne point of this kernel was also to allow faster and separate iterations of trying to build a 3D CNN to use as part of other models (e.g., with the tabular dataset).\n\nThe accuracy of the model in this kernel is low, but it could provide some useful building blocks for other kernels to build on. Or maybe someone will spot some errors and let me know..","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport cv2\n\nimport scipy\nfrom scipy import ndimage\nfrom scipy.ndimage.filters import gaussian_filter\nfrom scipy.ndimage import zoom\nimport random\n\nfrom tqdm.auto import tqdm\ntqdm.pandas()\n\nfrom skimage.io import imread\nfrom skimage.transform import resize\nfrom keras.utils import Sequence\nimport math\n\nimport PIL\nfrom PIL import Image, ImageOps\nimport matplotlib.pylab as plt\n\nfrom sklearn.model_selection import StratifiedKFold, GroupKFold, GroupShuffleSplit # Used to use Kfold to train our model\nfrom sklearn.model_selection import train_test_split\n\nimport tensorflow\nimport tensorflow as tf\n\nfrom tensorflow import keras\nfrom tensorflow.keras.datasets import mnist\nfrom tensorflow.keras.layers import Dense, Flatten, Input\nfrom tensorflow.keras.layers import Conv2D, MaxPooling2D, BatchNormalization, Dropout, Activation\nfrom tensorflow.keras.layers import Conv3D, MaxPooling3D\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import backend as K # The backend give us access to tensorflow operations and allow us to create the Attention class\n#from tensorflow.keras.models import Sequential\nimport matplotlib.pylab as plt\n\nfrom tensorflow.keras.callbacks import (ModelCheckpoint, LearningRateScheduler,\n                             EarlyStopping, ReduceLROnPlateau,CSVLogger)\n                             \nfrom sklearn.model_selection import train_test_split\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# A Few Configuration Variables","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"batch_size = 8\nepochs = 30\n#how many splits to do on the training data, or how many cross-validation rounds to use\nN_SPLITS = 5\n#my_test_pct is a percentage of values left out of test/validation data to compare the final model against known results\nmy_test_pct = 0.05\n\n#there are some different sizes of images in my rescaled dataset, so using some of those here\n#because the 3D arrays are quite large, they tend to take memory and this (192) was a size that did not cause out of memory errors\nimg_size = 192\n#img_size = 256\n#img_size = 512\nimg_depth = 30\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Brief Overview of the Data\n\nThe preprocessed data is in a [dataset](https://www.kaggle.com/donkeys/osic-pulmonary-fibrosispreprocessed) I previously uploaded. This is where it mounts:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls /kaggle/input/osic-pulmonary-fibrosispreprocessed","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this case I used a continous value for the smoking feature, although you can say it is categorical. But 0 for non-smoking, 1 for currently smoking, and 0.5 for used to smoke.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"df_train_orig = pd.read_csv(\"/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv\")\ndf_train = pd.read_csv(\"/kaggle/input/osic-pulmonary-fibrosispreprocessed/dataset/df_train_scaled_continous_smoke.csv\").drop(\"Unnamed: 0\", axis=1)\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train[\"SmokingStatus\"].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test = pd.read_csv(\"/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv\")\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/osic-pulmonary-fibrosispreprocessed/dataset\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train[df_train[\"Patient\"] == \"ID00126637202218610655908\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"List of patients in the training set. This can be useful, for example, to make patient-grouped data-splits.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ids = df_train[\"Patient\"].unique()\npatient_ids.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_id = \"ID00047637202184938901501\"\npatient_fvc = df_train[df_train[\"Patient\"] == patient_id]\npatient_fvc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#this is the original data from kaggle competition, no preprocessing done\n!ls /kaggle/input/osic-pulmonary-fibrosis-progression","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create Train and Test Rows\n\nThis is just a little something I ended up with to format my data after playing with various approaches. Maybe not really necessary but it works for this kernel.. I really just need the patient ID for accessing the image in this kernel, the rest are just someting I used in other kernels.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ids = df_train[\"Patient\"].unique()\ntraining_rows = []\npatient_count = 0\nfor patient_id in tqdm(patient_ids):\n    df_patient = df_train[df_train[\"Patient\"] == patient_id]\n    patient_row_count = 0\n    row = df_patient.iloc[0]\n    row_fvc = row[\"FVC\"]\n    patient_row_count += 1\n    training_row = {}\n    training_row[\"patient_id\"] = row[\"Patient\"]\n    training_row[\"base_week\"] = row[\"Weeks\"]\n    training_row[\"pct\"] = row[\"Percent\"]\n    training_row[\"age\"] = row[\"Age\"]\n    training_row[\"gender_female\"] = row[\"Sex_Female\"]\n    training_row[\"gender_male\"] = row[\"Sex_Male\"]\n    training_row[\"smoking_status\"] = row[\"SmokingStatus\"]\n    training_row[\"target_fvc\"] = row[\"fvc_raw\"]\n    training_rows.append(training_row)\n    patient_count += 1\nprint(f\"processed {patient_count} patients\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test = pd.read_csv(f\"{DATA_DIR}/df_test_scaled_continous_smoke.csv\")\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_ids = df_test[\"Patient\"].unique()\ntest_rows = []\npatient_count = 0\nfor patient_id in tqdm(patient_ids):\n    df_patient = df_test[df_test[\"Patient\"] == patient_id]\n    patient_row_count = 0\n    for idx, row in df_patient.iterrows():\n        row_fvc = row[\"FVC\"]\n        patient_row_count += 1\n        test_row = {}\n        test_row[\"patient_id\"] = row[\"Patient\"]\n        test_row[\"base_fvc\"] = row_fvc\n        test_row[\"base_week\"] = row[\"Weeks\"]\n        test_row[\"pct\"] = row[\"Percent\"]\n        test_row[\"age\"] = row[\"Age\"]\n        test_row[\"gender_female\"] = row[\"Sex_Female\"]\n        test_row[\"gender_male\"] = row[\"Sex_Male\"]\n        test_row[\"smoking_status\"] = row[\"SmokingStatus\"]\n        test_rows.append(test_row)\n    print(f\"created {patient_row_count} instances for patient {patient_id}\")\n    patient_count += 1\nprint(f\"processed {patient_count} patients\")\n    #break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_new_train = pd.DataFrame(training_rows)\ndf_new_train\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_new_train.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train[df_train[\"Patient\"] == \"ID00007637202177411956430\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I have scaled the columns to range 0..1, except the column now named *target_fvc*. I use this column as the target variable to predict here. It is the same as the base fvc given for each patient. Just renamed target here, as it is the prediction target..","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#the only thing this kernel uses from the x_cols is actually the patient id, as that can be used to find the filename of the image.\n#some other kernels i publish use the other colums as well, which is why they are there\nx_cols = [col for col in df_new_train.columns if col != \"target_fvc\"]\ndf_x = df_new_train[x_cols]\ndf_y = df_new_train[\"target_fvc\"]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_x.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_y.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Custom Keras Sequence Generator\n\nHere I define a Keras Sequence generator to build the augmented dataset on the fly.\n\nFirst, a small utility function to shuffle x and y at the end of each training epoch:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def unison_shuffled_copies(a, b):\n    assert len(a) == len(b)\n    p = np.random.permutation(len(a))\n    new_a = a.iloc[p]\n    new_b = b.iloc[p]\n    return new_a, new_b","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A custom Keras generator sequence to produce batches of 3D augmented images:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class MySequence3D(Sequence):\n\n    def __init__(self, x_set, y_set, batch_size, mode=\"train\", augment=True):\n        self.x, self.y = x_set, y_set\n        self.batch_size = batch_size\n        self.mode = mode\n        self.max_idx = math.ceil(len(x_set)/batch_size)\n        self.augment = augment\n\n    def __len__(self):\n        #TODO: check is correct\n        return int(np.ceil(len(self.x) / float(self.batch_size)))\n\n    def __getitem__(self, idx):\n        idx2 = idx % self.max_idx\n        start = idx2 * self.batch_size\n        end = min(start + batch_size, len(self.x))\n        batch_x = self.x.iloc[start : end]\n        batch_y = self.y.iloc[start : end]\n        \n        next_batch_img = []\n        next_batch_num = []\n        for index, row in batch_x.iterrows():\n            #print(row)\n            #nums = [row[col] for col in dense_cols]\n            nums = []\n            nums = np.array(nums)\n            file_name = row[\"patient_id\"]\n            #file_path = row[\"path\"]\n            augmented = img_augment_3d(self.mode, self.x, file_name, self.augment)\n            next_batch_num.append(nums)\n            next_batch_img.append(augmented)\n        np_y = np.array(batch_y)\n        np_x_img = np.array(next_batch_img)\n        np_x_num = np.array(next_batch_num)\n        del next_batch_img\n        del next_batch_num\n        del batch_y\n        #print(f\"loaded shape: {np_x.shape}, batch={idx}\")\n#        result = [np_x_img, np_x_num], np_y\n        result = [np_x_img], np_y\n\n        #print(f\"shapes: {result[0].shape}, {result[1].shape}\")\n        return result\n\n    def on_epoch_end(self):\n        self.x, self.y = unison_shuffled_copies(self.x, self.y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3D Augmentation Helpers\n\nA number of functions to perform augmentations, one per function, on given 3D images. These mostly start with a random chance of whether it will be applied at all or not. Just to avoid applying all augmentations all the time. Then the augmentation itself.\n\n### Rotate\n\nFirst, a helper to rotage a given 3D image by a random angle between given min and max angle:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#https://stackoverflow.com/questions/43922198/how-to-rotate-a-3d-image-by-a-random-angle-in-python\ndef random_rotation_3d(img, min_angle, max_angle):\n    \"\"\" Randomly rotate an image by a random angle (-max_angle, max_angle).\n\n    Arguments:\n    max_angle: `float`. The maximum rotation angle.\n\n    Returns:\n    rotated 3D image\n    \"\"\"\n    if random.randint(1,100) > 30:\n        #with some chance, do not rotate at all\n        return img\n    img_rot = np.zeros(img.shape)\n    angle = random.uniform(min_angle, max_angle)\n    if random.randint(1,100) > 50:\n        #in half the cases, rotate left. in other half, rotate right.\n        angle *= -1\n        # Following lines would rotate on z and y axis as well, but not using them in this kernel\n#        # rotate along z-axis\n#        image2 = scipy.ndimage.interpolation.rotate(image1, angle, mode='nearest', axes=(0, 1), reshape=False)\n#        # rotate along y-axis\n#        image3 = scipy.ndimage.interpolation.rotate(image2, angle, mode='nearest', axes=(0, 2), reshape=False)\n\n    # rotate along x-axis\n    img_rot = scipy.ndimage.interpolation.rotate(img, angle, mode='nearest', axes=(1, 2), reshape=False)\n    return img_rot.reshape(img.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Gaussian Blur\n\nA helper to do a gaussian blur on a 3D array (of image pixels):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#https://stackoverflow.com/questions/29920114/how-to-gauss-filter-blur-a-floating-point-numpy-array\ndef gaussian_blur_3d(img):\n    if random.randint(1,100) > 15:\n        return img\n    sigma = random.uniform(0.1,0.9)\n    blurred = gaussian_filter(img, sigma=sigma)\n    return blurred\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Horizontal Flip\n\nFlip an entire 3D array along x-axis:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#https://stackoverflow.com/questions/7416170/numpy-reverse-multidimensional-array\ndef x_flip(img):\n    if random.randint(1,100) > 50:\n        flipped = img[:, :, ::-1]\n    else:\n        flipped = img\n    return flipped","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Vertical Flip\n\nFlip an entire 3D array along y-axis (some images seem to be upside down as noted in my [previous kernel](https://www.kaggle.com/donkeys/ct-slices-basic-eda)):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def y_flip(img):\n    if random.randint(1,100) > 70:\n        flipped = img[:, ::-1, :]\n    else:\n        flipped = img\n    return flipped","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Shift X-axis\n\nShift a 3D array on x-axis by random number of pixels between given min and max. Left or right, random choice of direction.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def x_shift(img, min_shift, max_shift):\n    if random.randint(1,100) > 30:\n        return img\n    shift_dir = 1\n    if random.randint(1,100) > 50:\n        shift_dir = -1\n    roll_amount = random.randint(min_shift, max_shift)\n    roll_amount *= shift_dir\n    img = np.roll(img, roll_amount, axis=2)\n    #z,y,x?\n    if shift_dir > 0:\n        img[:, :, 0:roll_amount] = 0\n    else:\n        img[:, :, roll_amount:] = 0\n#    print(img)\n    return img\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Shift Y-axis\n\nShift a 3D array on y-axis by random number of pixels between given min and max. Left or right, random choice of direction.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def y_shift(img, min_shift, max_shift):\n    if random.randint(1,100) > 30:\n        return img\n    shift_dir = 1\n    if random.randint(1,100) > 50:\n        shift_dir = -1\n    roll_amount = random.randint(min_shift, max_shift)\n    roll_amount *= shift_dir\n    img = np.roll(img, roll_amount, axis=1)\n    #z,y,x?\n    if shift_dir > 0:\n        img[:, 0:roll_amount, :] = 0\n    else:\n        img[:, roll_amount:, :] = 0\n#    print(img)\n    return img\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Zoom\n\nZoom the 3D image by a random zoom factor (same for both axis) between given min and max on both x- and y-axis:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#https://stackoverflow.com/questions/37119071/scipy-rotate-and-zoom-an-image-without-changing-its-dimensions\ndef zoom_xy(img, min_zoom, max_zoom):\n    if random.randint(1,100) > 20:\n        return img\n    zoom_factor = random.uniform(min_zoom, max_zoom)\n    h, w = img_size, img_size\n\n    # For multichannel images we don't want to apply the zoom factor to the RGB\n    # dimension, so instead we create a tuple of zoom factors, one per array\n    # dimension, with 1's for any trailing dimensions after the width and height.\n    zoom_tuple = (1, zoom_factor, zoom_factor)\n\n    # Zooming out\n    if zoom_factor < 1:\n\n        # Bounding box of the zoomed-out image within the output array\n        zh = int(np.round(h * zoom_factor))\n        zw = int(np.round(w * zoom_factor))\n        top = (h - zh) // 2\n        left = (w - zw) // 2\n\n        # Zero-padding\n        out = np.zeros_like(img)\n        zoomed_img = zoom(img, zoom_tuple, order=0)\n        #print(f\"zoomed shape: {zoomed_img.shape}\")\n        #print(f\"out shape:{out.shape}\")\n        #print(f\"w:{w},h:{h},l:{left},t:{top},zw:{zw}, zh:{zh}\")\n        out[:, top:top+zh, left:left+zw] = zoomed_img\n\n    # Zooming in\n    elif zoom_factor > 1:\n\n        # Bounding box of the zoomed-in region within the input array\n        zh = int(np.ceil(h / zoom_factor))\n        zw = int(np.ceil(w / zoom_factor))\n        top = (h - zh) // 2\n        left = (w - zw) // 2\n\n        #out_template = np.zeros_like(img)\n        out = zoom(img[:, top:top+zh, left:left+zw], zoom_tuple, order=0)\n        #print(f\"out shape:{out.shape}\")\n        #print(f\"w:{w},h:{h},l:{left},t:{top},zw:{zw}, zh:{zh}\")\n\n        # `out` might still be slightly larger than `img` due to rounding, so\n        # trim off any extra pixels at the edges\n        trim_top = ((out.shape[1] - h) // 2)\n        trim_left = ((out.shape[2] - w) // 2)\n        #print(f\"out shape before:{out.shape}\")\n        out = out[:, trim_top:trim_top+h, trim_left:trim_left+w]\n        #print(f\"out shape after:{out.shape}\")\n        #print(f\"w:{w},h:{h},l:{left},trimtop:{trim_top},trimleft:{trim_left}\")\n\n    # If zoom_factor == 1, just return the input array\n    else:\n        out = img\n    #print(out.shape)\n    return out","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Failed Attempt at Zoom\n\nThis was an attempt at using CV2 library for zooming, keeping it here for posterity:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#https://stackoverflow.com/questions/37119071/scipy-rotate-and-zoom-an-image-without-changing-its-dimensions\n#open cv does not seem to support 3d image resizing so cannot do that..\ndef cv2_zoom_xy(img, min_zoom, max_zoom):\n    if random.randint(1,100) > 20:\n        return img\n    zoom_factor = random.uniform(min_zoom, max_zoom)\n    zoom_factor = 2.0\n    print(f\"zf:{zoom_factor}\")\n    height, width = img_size, img_size\n    new_height, new_width = int(height * zoom_factor), int(width * zoom_factor)\n\n    ### Crop only the part that will remain in the result (more efficient)\n    # Centered bbox of the final desired size in resized (larger/smaller) image coordinates\n    y1, x1 = max(0, new_height - height) // 2, max(0, new_width - width) // 2\n    y2, x2 = y1 + height, x1 + width\n    bbox = np.array([y1,x1,y2,x2])\n    # Map back to original image coordinates\n    bbox = (bbox / zoom_factor).astype(np.int)\n    y1, x1, y2, x2 = bbox\n    cropped_img = img[:, y1:y2, x1:x2]\n\n    # Handle padding when downscaling\n    resize_height, resize_width = min(new_height, height), min(new_width, width)\n    pad_height1, pad_width1 = (height - resize_height) // 2, (width - resize_width) //2\n    pad_height2, pad_width2 = (height - resize_height) - pad_height1, (width - resize_width) - pad_width1\n    pad_spec = [(0,0), (pad_height1, pad_height2), (pad_width1, pad_width2)] + [(0,0)] * (img.ndim - 2)\n\n    result = cv2.resize(cropped_img, (img.shape[0], resize_width, resize_height))\n    result = np.pad(result, pad_spec, mode='constant')\n    assert result.shape[0] == height and result.shape[1] == width\n    return result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The One Function to Bind them Together\n\nThe actual augmentation function that aggregates all of the above into one:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def img_augment_3d(df_name, df, patient_id, do_augment):\n    filename = f\"{DATA_DIR}/scaled_png/{df_name}_{img_depth}/{patient_id}/full_3d_{img_size}.npy\"\n    img = np.load(filename)\n    img = img / 255.0\n    if do_augment:\n        #give it 15% chance of not doing any augmentation\n        if random.randint(1,100) > 15:\n            #you can just comment all but one of below calls to just see effects of a single one later\n            #in which case, you might want to modify above if the > 0 to make it always true and see changes\n            img = gaussian_blur_3d(img)\n            img = x_flip(img)\n            img = y_flip(img)\n            img = random_rotation_3d(img, 1, 5)\n            img = x_shift(img, 5, 15)\n            img = y_shift(img, 5, 15)\n            img = zoom_xy(img, 0.9, 1.1)\n            #pass\n    return img\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Smoke Test\n\nSometimes the custom sequence generator has issues when some combination of augmentations hits, so this allows testing it with various combinations to see if it crashes (since it loops the batch generator with augmentation, it could also be used as a performance test):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"experiment_gen = MySequence3D(df_x, df_y, batch_size, augment=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for x in tqdm(range(10)):\n    batch = experiment_gen.__getitem__(0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualizing Augmentations to See How They Work\n\nA helper to plot the 3D augmented images the above custom generator produces (for debugging/testing):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_gen_batch(generator, idx):\n    # configure batch size and retrieve one batch of images\n    plt.clf() #clears matplotlib data and axes\n    #for batch in train_generator:\n    rows = (batch_size / 3)+1\n    plt.figure(figsize=[30,10*rows])\n    batch = generator.__getitem__(idx)\n    print(f\"showing {len(batch[0])} images\")\n    #have to use len(batch[0] here, as batch size can vary if it is the last part of the images (truncated to dataset length)\n    for x in range(0, len(batch[0][0])):\n    #    print(train_generator.filenames[x])\n    \n        plt.subplot(rows, 3, x+1)\n        #batch[0] is x, which is [imgs, nums]\n        #batch[0][0] is imgs\n        #so the following line just takes the first slice of each 3D image in the batch\n        img_2d_plane = batch[0][0][x][0]\n        plt.imshow(img_2d_plane, interpolation='nearest')\n\n        num = \"disabled\"\n        y = batch[1][x]\n        print(f\"num: {num}, y: {y}, img min: {img_2d_plane.min()} max: {img_2d_plane.max()}\")\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#this gives length 2, since it is an array of (x,y) batch items\nlen(batch)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#with batch size 8 we get 8 images\nlen(batch[0][0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#with img_depth=30 and img_size=192, we get 8 images in a batch, each of size depth(z)=30, height(y)=192, width(x)=192\nbatch[0][0].shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot Some Example Augmented Batches\n\nEvery time `__getitem__()` is called on the generator, it should produce differenct augmentations on the fly. So here we call it twice to see that it works as expected:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_gen_batch(experiment_gen, 0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_gen_batch(experiment_gen, 0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Just a little cleanup as we are resource-wasting paranoid:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"del batch\ndel experiment_gen","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create the Model\n\nHere, I build the 3D CNN to play with. First basic config for image sizes etc:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"h = img_size\nw = img_size\ninput_shape = (h, w, img_depth, 1) #the image has just one color channel\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The actual `create_model()` function further below will call `create_cnn_model()` to create the CNN. This is because some other kernels I reused the model I build in this kernel and this allows me to combine them easier. I built this here separately since it allowed me to focus on the CNN only, and run it faster with a smaller dataset. To iterate model architectures faster and with less use of the limited GPU resources on Kaggle.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_cnn_model():\n    img_input = Input(shape=(input_shape), name=\"img_input\")\n    cnn = Conv3D(32, kernel_size=(15), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_1\")(img_input)\n    #cnn = Dropout(0.5)(cnn)\n    #strides=(1,1)?\n    #batchnorm seems to really wreck havoc on the results if added anywhere in the cnn\n    #cnn = BatchNormalization()(cnn)\n    #cnn = MaxPooling3D(pool_size=(2,2,2))(cnn) #strides=(2,2)?\n    #cnn = Dropout(0.6)(cnn)\n    cnn = Conv3D(32, kernel_size=(7), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_2\")(img_input)\n    #cnn = Dropout(0.5)(cnn)\n    #strides=(1,1)?\n    #cnn = BatchNormalization()(cnn)\n    cnn = MaxPooling3D(pool_size=(2,2,2), strides=(2, 2, 2))(cnn) #strides=(2,2)?\n    #cnn = Dropout(0.45)(cnn)\n    cnn = Conv3D(64, kernel_size=(3), strides=(1), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_3\")(cnn)\n    cnn = Conv3D(64, kernel_size=(3), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_4\")(cnn)\n    #cnn = BatchNormalization()(cnn)\n    #cnn = MaxPooling3D(pool_size=(2,2,2), strides=(2, 2, 2))(cnn) #strides=(2,2)?\n    #cnn = Dropout(0.45)(cnn)\n    #cnn = Conv3D(128, kernel_size=(5), strides=(1), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_5\")(cnn)\n    #cnn = Conv3D(128, kernel_size=(3), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_6\")(cnn)\n    #cnn = BatchNormalization()(cnn)\n    cnn = MaxPooling3D(pool_size=(2,2,2), strides=(2, 2, 2))(cnn) #strides=(2,2)?\n    #cnn = Dropout(0.3)(cnn)\n    flatten = Flatten()(cnn)\n    final_cnn_dense = Dense(100, activation='relu')(flatten)\n    model = keras.Model(\n        inputs=[img_input],\n        outputs=[final_cnn_dense],\n    )\n    return model\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Another alternative I tried, just for illustration:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#tried a few options, just rename the one to try as \"create_cnn_model\" and run the thing..\ndef create_cnn_model_small():\n    img_input = Input(shape=(input_shape), name=\"img_input\")\n    cnn = Conv3D(32, kernel_size=(5), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_1\")(img_input)\n    cnn = MaxPooling3D(pool_size=(2,2,2), strides=(2))(cnn)\n    cnn = Conv3D(64, kernel_size=(3), strides=(2), padding='same', activation='relu', kernel_initializer='he_uniform', name=\"conv_3d_2\")(img_input)\n    cnn = MaxPooling3D(pool_size=(2,2,2), strides=(2, 2, 2))(cnn)\n    flatten = Flatten()(cnn)\n    final_cnn_dense = Dense(100, activation='relu')(flatten)\n    model = keras.Model(\n        inputs=[img_input],\n        outputs=[final_cnn_dense],\n    )\n    return model\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This creates the actual model for training folds:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.optimizers import Adam\n\ndef create_model():\n    cnn_model = create_cnn_model()\n    x = Dropout(0.5)(cnn_model.output)\n    x = Dense(200, activation=\"relu\")(x)\n    x = Dropout(0.5)(x)\n    x = Dense(1, activation=\"linear\", name=\"final_dense\")(x)\n    model = keras.Model(inputs=[cnn_model.input], outputs=x)\n    adam = tensorflow.keras.optimizers.Adam(lr=0.01)\n#    adam = tensorflow.keras.optimizers.Adam(lr=0.001)\n#    adam = tensorflow.keras.optimizers.Adam(lr=1e-2, decay=1e-2/epochs)\n\n    model.compile(loss='mean_squared_error',\n                  optimizer=adam,  #keras.optimizers.SGD(lr=0.01),\n                  metrics=[tf.keras.metrics.RootMeanSquaredError(name='rmse')])\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And our regular Keras callbacks to save the best model during training, lower training rate on plateuau, and stop a bit earlier if no gains:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_callbacks(idx):\n    checkpoint = ModelCheckpoint(f'../working/weights_best_{idx}.h5', monitor='val_loss', verbose=1, \n                                 save_best_only=True, mode='min', save_weights_only = True)\n    reduceLROnPlat = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=3, min_lr=1e-6,\n                                       verbose=1, mode='auto', epsilon=0.0001)\n    early = EarlyStopping(monitor=\"val_loss\", \n                          mode=\"min\", \n                          patience=22) #OK, patience is silly / too high for this number of epochs but whatever :)\n\n    csv_logger = CSVLogger(filename='../working/training_log.csv',\n                           separator=',',\n                           append=True)\n\n    callbacks_list = [checkpoint, reduceLROnPlat, csv_logger, early]\n    return callbacks_list","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Function to Train the Model on Given Data\n\nThis function runs the training for a single fold with given data:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def fit_model(model, callbacks_list, df_x_train, df_y_train, df_x_val, df_y_val):\n    train_gen = MySequence3D(df_x_train, df_y_train, batch_size, augment=True)\n    valid_gen = MySequence3D(df_x_val, df_y_val, batch_size, augment=False)\n\n    #the total number of images we have:\n    train_size = df_x_train.shape[0]\n    #train_steps is how many steps per epoch Keras runs the genrator. One step is batch_size*images\n    train_steps = train_size/batch_size\n    train_steps = int(train_steps)\n    #same for the validation set\n    valid_size = df_x_val.shape[0]\n    valid_steps = valid_size/batch_size\n    valid_steps = int(valid_steps)\n\n    fit_history = model.fit_generator(\n            train_gen,\n            steps_per_epoch=train_steps,\n            epochs = epochs,\n            validation_data=valid_gen,\n            validation_steps=valid_steps,\n            callbacks=callbacks_list,\n        use_multiprocessing=True,\n        workers=2,\n        verbose = 1\n    )\n    return fit_history","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Make Predictions\n\nMake predictions for the Kaggle test set, using target week numbers -12...133. This is just to see how the prediction runs here, did not use this for actual submission, since I made this kernel more to explore the CNN architecture. It does not use any time related data, or the tabular data, so its not really very good for the actual final prediction over all weeks. But keeping these here anyway.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_predictions_dict(model, test_rows):\n    if isinstance(test_rows, pd.DataFrame):\n        test_rows = [row.to_dict() for (idx, row) in test_rows.iterrows()]\n    predictions = {}\n    col_names = []\n    print(f\"predicting {test_rows}\")\n    for target_week in tqdm(range(-12,134)):\n        for idx, row in enumerate(test_rows):\n            row[\"target_week\"] = (target_week+12)/(133+12)\n            patient_id = row[\"patient_id\"]\n            img = img_augment_3d(\"train\", None, patient_id, False)\n            img = np.array([img])\n            #print(img)\n            #print(img.shape)\n            pred = model.predict([img])\n            col_name = f\"{idx+1}_{target_week}\"\n            col_names.append(col_name)\n            predictions[col_name] = pred.flatten()[0]\n            #print(f\"target week: {row['target_week']}, pred: {pred}\")\n        \n    return predictions, col_names","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Predict also for my own test set that was put aside from the training data. So I can compare later how well/bad it did with the final trained model (since training set has known target FVC values):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_predictions_my_test(model, test_rows):\n    test_rows = [row.to_dict() for (idx, row) in test_rows.iterrows()]\n    predictions = []\n    for idx, row in tqdm(enumerate(test_rows), total=len(test_rows)):\n        patient_id = row[\"patient_id\"]\n        img = img_augment_3d(\"train\", None, patient_id, False)\n        img = np.array([img])\n        pred = model.predict([img])\n        predictions.append(pred.flatten()[0])\n        #print(f\"target week: {row['target_week']}, pred: {pred}\")\n        \n    return predictions","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Split into N Groups to Train N Models","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"indices = np.arange(df_x.shape[0])\n#https://github.com/scikit-learn/scikit-learn/issues/9193\ntrain_indices, my_test_indices = next(GroupShuffleSplit(test_size=my_test_pct, random_state=8).split(indices, groups=df_x[\"patient_id\"]))\nmy_test_X = df_x.iloc[my_test_indices]\nmy_test_y = df_y.iloc[my_test_indices]\n\nfull_indices = indices\nindices = train_indices","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Run the Training on Model","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\n#split the training data by patient id, so if patient has multiple rows we put them in the same group\nsplits = list(GroupKFold(n_splits=N_SPLITS).split(indices, groups=df_x.iloc[indices][\"patient_id\"]))\npreds_test = []\npreds_my_test = []\ncol_names = []\nfit_histories = []\n\nfor idx, (train_idx, val_idx) in enumerate(splits):\n    K.clear_session() # start Keras from clean state in each iteration\n    print(\"Beginning fold {}\".format(idx+1))\n    # use the indexes to extract the folds in the train and validation data\n    train_X, train_y, val_X, val_y = df_x.iloc[train_idx], df_y.iloc[train_idx], df_x.iloc[val_idx], df_y.iloc[val_idx]\n    # instantiate the model for this fold\n    model = create_model()\n    callbacks = create_callbacks(idx)\n    #train the model on this fold\n    history = fit_model(model, callbacks, train_X, train_y, val_X, val_y)\n    fit_histories.append(history)\n    # loads the best weights saved by the checkpoint\n    model.load_weights(f'weights_best_{idx}.h5')\n    pred_test, col_names = make_predictions_dict(model, test_rows)\n    preds_test.append(pred_test)\n    pred_test = make_predictions_my_test(model, my_test_X)\n    preds_my_test.append(pred_test)\n    del model\n    del callbacks\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"my_test_X\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = create_model()\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.utils import plot_model\nplot_model(model)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot Training Losses per Fold","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_loss_and_accuracy(fit_history, n=0):\n    plt.clf()\n    plt.plot(fit_history.history['rmse'][n:])\n    plt.plot(fit_history.history['val_rmse'][n:])\n    plt.title('model rmse')\n    plt.ylabel('rmse')\n    plt.xlabel('epoch')\n    plt.legend(['train', 'test'], loc='upper left')\n    plt.show()\n    plt.clf()\n    #since the loss is just mse here, the graph is practically identical to RMSE (metric above). Save space and plot it just once..\n    # summarize history for loss\n    #plt.plot(fit_history.history['loss'][n:])\n    #plt.plot(fit_history.history['val_loss'][n:])\n    #plt.title('model loss')\n    #plt.ylabel('loss')\n    #plt.xlabel('epoch')\n    #plt.legend(['train', 'test'], loc='upper left')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Full loss history for all folds:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"for idx, fit_history in enumerate(fit_histories):\n    print(f\"fold {idx}\")\n    plot_loss_and_accuracy(fit_history)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The above charts often have a few versions where the first few epochs have very high loss, and then these drop down significatnly. But because the start is high, variance in the later epochs is hard to spot in the chart. To avoid this, here we plot loss history starting at epoch 4 for all folds:\n    ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"for idx, fit_history in enumerate(fit_histories):\n    print(f\"fold {idx}\")\n    plot_loss_and_accuracy(fit_history, 4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Look at Custom Test Set Predictions vs Actual Values","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"len(preds_my_test[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds_my_test_mean = np.mean(preds_my_test, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\n\nmean_squared_error(my_test_y, preds_my_test_mean, squared=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Convert our test set predictions into a dataframe for easier manipulation:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions = pd.DataFrame()\nfor idx, mtp in enumerate(preds_my_test):\n    df_mytest_predictions[f\"{idx+1}\"] = mtp\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions.T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions.T.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions.T.describe().T.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Check Diff per Patient in Test Set\n\nBy looking at how much the FVC predictions differ from the actual FVC value per patient, maybe we could identify some trends.\n\nLets see.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions_diff = pd.DataFrame()\nfor idx, mtp in enumerate(preds_my_test):\n    df_mytest_predictions_diff[f\"{idx+1}\"] = np.abs(np.array(mtp) - my_test_y.values)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions_diff.T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Above shows one row per fold (5 folds = 5 rows), where each column is the prediction difference vs actual value for a patient (0-8 patients, for a total of 9 patients in this test set).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions_diff.T.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Above shows the mean, str, etc over the 5 folds for all patients in the test set.\n\nBy using describe() on describe(), we can get aggregated statistics on the overall dataset (so mean over all patients means, etc):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_mytest_predictions_diff.T.describe().T.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So on average, when I ran this, we were off across all folds and patients by 566 units. The results on this posted version should be little bit different but not too much, due to randomness of the process. When the range of units is roughly about 2000-3000, I'd say thats a pretty bad result.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Predictions for Kaggle Test Set\n\nThis would be the test set Kaggle provices, and the one that would be providing the actual competition results if we submitted this. But not going there with these poor results. Lets play anyway.\n\nPartly my idea was to look at the differences above, and compare with these values to see if I could find some trend to apply. Of course not.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test_predictions = pd.DataFrame(preds_test)\ndf_test_predictions = df_test_predictions[col_names] #to get sorted order","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test_predictions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test_predictions.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"descriptions = df_test_predictions.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"descriptions.max(axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"descriptions.T.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}