{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport seaborn as sns\nimport pymc3 as pm\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold\nfrom sklearn import datasets, linear_model\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\nimport statistics\nimport warnings\nimport copy\nfrom datetime import timedelta, datetime\nimport imageio\nfrom matplotlib import cm\n# import multiprocessing\nimport os\nfrom pathlib import Path\nimport pydicom\nimport pytest\nimport scipy.ndimage as ndimage\nfrom scipy.ndimage.interpolation import zoom\nfrom skimage import measure, morphology, segmentation\nfrom time import time, sleep\nfrom tqdm import trange, tqdm\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, random_split, DistributedSampler, DataLoader\nfrom torch.utils.tensorboard import SummaryWriter\nfrom torchvision import transforms\n\nfrom scipy.stats import skew, kurtosis\nfrom skimage.measure import label, regionprops\nfrom skimage.segmentation import clear_border\nimport multiprocessing as mp\nfrom functools import partial\nimport sys\n\nwarnings.filterwarnings('ignore')\npd.set_option('display.max_rows', 176)\npd.set_option('display.max_columns', 35)\n\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        break\n#print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(sys.version)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Autoencoder code:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# !ls ../","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# root_dir = '/kaggle/input/osic-cached-dataset'\n# train_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/train'\n# test_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/test'\n# model_file = '/kaggle/working/diophantus.pt'\n# resize_dims = (40, 256, 256)\n# clip_bounds = (-1000, 200)\n# watershed_iterations = 1\n# pre_calculated_mean = 0.02865046213070556\n# latent_features = 2\n# batch_size = 16\n# learning_rate = 3e-5\n# num_epochs = 10\n# val_size = 0\n# tensorboard_dir = '/kaggle/working/runs'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# class CTScansDataset(Dataset):\n#     def __init__(self, root_dir, transform=None):\n#         self.root_dir = Path(root_dir)\n#         self.patients = [p for p in self.root_dir.glob('*') if p.is_dir()]\n#         self.transform = transform\n\n#     def __len__(self):\n#         return len(self.patients)\n\n#     def __getitem__(self, idx):\n#         if torch.is_tensor(idx):\n#             idx = idx.tolist()\n\n#         image, metadata = self.load_scan(self.patients[idx])\n#         sample = {'image': image, 'metadata': metadata}\n#         if self.transform:\n#             sample = self.transform(sample)\n\n#         return sample\n\n#     def save(self, path):\n#         t0 = time()\n#         Path(path).mkdir(exist_ok=True, parents=True)\n#         print('Saving pre-processed dataset to disk')\n#         sleep(1)\n#         cum = 0\n\n#         bar = trange(len(self))\n#         for i in bar:\n#             sample = self[i]\n#             image, data = sample['image'], sample['metadata']\n#             cum += torch.mean(image).item()\n\n#             bar.set_description(f'Saving CT scan {data.PatientID}')\n#             fname = Path(path) / f'{data.PatientID}.pt'\n#             torch.save(image, fname)\n\n#         sleep(1)\n#         bar.close()\n#         print(f'Done! Time {timedelta(seconds=time() - t0)}\\n'\n#               f'Mean value: {cum / len(self)}')\n\n#     def get_patient(self, patient_id):\n#         patient_ids = [str(p.stem) for p in self.patients]\n#         return self.__getitem__(patient_ids.index(patient_id))\n\n#     @staticmethod\n#     def load_scan(path):\n#         slices = [pydicom.read_file(p) for p in path.glob('*.dcm')]\n#         try:\n#             slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n#         except AttributeError:\n#             warnings.warn(f'Patient {slices[0].PatientID} CT scan does not '\n#                           f'have \"ImagePositionPatient\". Assuming filenames '\n#                           f'in the right scan order.')\n\n#         image = np.stack([s.pixel_array.astype(float) for s in slices])\n#         return image, slices[0]\n\n# class CropBoundingBox:\n#     @staticmethod\n#     def bounding_box(img3d: np.array):\n#         mid_img = img3d[int(img3d.shape[0] / 2)]\n#         same_first_row = (mid_img[0, :] == mid_img[0, 0]).all()\n#         same_first_col = (mid_img[:, 0] == mid_img[0, 0]).all()\n#         if same_first_col and same_first_row:\n#             return True\n#         else:\n#             return False\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n#         if not self.bounding_box(image):\n#             return sample\n\n#         mid_img = image[int(image.shape[0] / 2)]\n#         r_min, r_max = None, None\n#         c_min, c_max = None, None\n#         for row in range(mid_img.shape[0]):\n#             if not (mid_img[row, :] == mid_img[0, 0]).all() and r_min is None:\n#                 r_min = row\n#             if (mid_img[row, :] == mid_img[0, 0]).all() and r_max is None \\\n#                     and r_min is not None:\n#                 r_max = row\n#                 break\n\n#         for col in range(mid_img.shape[1]):\n#             if not (mid_img[:, col] == mid_img[0, 0]).all() and c_min is None:\n#                 c_min = col\n#             if (mid_img[:, col] == mid_img[0, 0]).all() and c_max is None \\\n#                     and c_min is not None:\n#                 c_max = col\n#                 break\n\n#         image = image[:, r_min:r_max, c_min:c_max]\n#         return {'image': image, 'metadata': data}\n\n# class ConvertToHU:\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n\n#         img_type = data.ImageType\n#         is_hu = img_type[0] == 'ORIGINAL' and not (img_type[2] == 'LOCALIZER')\n#         # if not is_hu:\n#         #     warnings.warn(f'Patient {data.PatientID} CT Scan not cannot be'\n#         #                   f'converted to Hounsfield Units (HU).')\n\n#         intercept = data.RescaleIntercept\n#         slope = data.RescaleSlope\n#         image = (image * slope + intercept).astype(np.int16)\n#         return {'image': image, 'metadata': data}\n    \n# class Resize:\n#     def __init__(self, output_size):\n#         assert isinstance(output_size, tuple)\n#         self.output_size = output_size\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n#         resize_factor = np.array(self.output_size) / np.array(image.shape)\n#         image = zoom(image, resize_factor, mode='nearest')\n#         return {'image': image, 'metadata': data}\n    \n# class Clip:\n#     def __init__(self, bounds=(-1000, 500)):\n#         self.min = min(bounds)\n#         self.max = max(bounds)\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n#         image[image < self.min] = self.min\n#         image[image > self.max] = self.max\n#         return {'image': image, 'metadata': data}\n    \n# class MaskWatershed:\n#     def __init__(self, min_hu, iterations, show_tqdm):\n#         self.min_hu = min_hu\n#         self.iterations = iterations\n#         self.show_tqdm = show_tqdm\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n\n#         stack = []\n#         if self.show_tqdm:\n#             bar = trange(image.shape[0])\n#             bar.set_description(f'Masking CT scan {data.PatientID}')\n#         else:\n#             bar = range(image.shape[0])\n#         for slice_idx in bar:\n#             sliced = image[slice_idx]\n#             stack.append(self.seperate_lungs(sliced, self.min_hu,\n#                                              self.iterations))\n\n#         return {\n#             'image': np.stack(stack),\n#             'metadata': sample['metadata']\n#         }\n\n#     @staticmethod\n#     def seperate_lungs(image, min_hu, iterations):\n#         h, w = image.shape[0], image.shape[1]\n\n#         marker_internal, marker_external, marker_watershed = MaskWatershed.generate_markers(image)\n\n#         # Sobel-Gradient\n#         sobel_filtered_dx = ndimage.sobel(image, 1)\n#         sobel_filtered_dy = ndimage.sobel(image, 0)\n#         sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n#         sobel_gradient *= 255.0 / np.max(sobel_gradient)\n\n#         watershed = morphology.watershed(sobel_gradient, marker_watershed)\n\n#         outline = ndimage.morphological_gradient(watershed, size=(3,3))\n#         outline = outline.astype(bool)\n\n#         # Structuring element used for the filter\n#         blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n#                            [0, 1, 1, 1, 1, 1, 0],\n#                            [1, 1, 1, 1, 1, 1, 1],\n#                            [1, 1, 1, 1, 1, 1, 1],\n#                            [1, 1, 1, 1, 1, 1, 1],\n#                            [0, 1, 1, 1, 1, 1, 0],\n#                            [0, 0, 1, 1, 1, 0, 0]]\n\n#         blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n\n#         # Perform Black Top-hat filter\n#         outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n\n#         lungfilter = np.bitwise_or(marker_internal, outline)\n#         lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n\n#         segmented = np.where(lungfilter == 1, image, min_hu * np.ones((h, w)))\n\n#         return segmented  #, lungfilter, outline, watershed, sobel_gradient\n\n#     @staticmethod\n#     def generate_markers(image, threshold=-400):\n#         h, w = image.shape[0], image.shape[1]\n\n#         marker_internal = image < threshold\n#         marker_internal = segmentation.clear_border(marker_internal)\n#         marker_internal_labels = measure.label(marker_internal)\n\n#         areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n#         areas.sort()\n\n#         if len(areas) > 2:\n#             for region in measure.regionprops(marker_internal_labels):\n#                 if region.area < areas[-2]:\n#                     for coordinates in region.coords:\n#                         marker_internal_labels[coordinates[0], coordinates[1]] = 0\n\n#         marker_internal = marker_internal_labels > 0\n\n#         # Creation of the External Marker\n#         external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n#         external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n#         marker_external = external_b ^ external_a\n\n#         # Creation of the Watershed Marker\n#         marker_watershed = np.zeros((h, w), dtype=np.int)\n#         marker_watershed += marker_internal * 255\n#         marker_watershed += marker_external * 128\n\n#         return marker_internal, marker_external, marker_watershed\n    \n# class Normalize:\n#     def __init__(self, bounds=(-1000, 500)):\n#         self.min = min(bounds)\n#         self.max = max(bounds)\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n#         image = image.astype(np.float)\n#         image = (image - self.min) / (self.max - self.min)\n#         return {'image': image, 'metadata': data}\n    \n\n# class ToTensor:\n#     def __init__(self, add_channel=True):\n#         self.add_channel = add_channel\n\n#     def __call__(self, sample):\n#         image, data = sample['image'], sample['metadata']\n#         if self.add_channel:\n#             image = np.expand_dims(image, axis=0)\n\n#         return {'image': torch.from_numpy(image), 'metadata': data}\n    \n    \n# class ZeroCenter:\n#     def __init__(self, pre_calculated_mean):\n#         self.pre_calculated_mean = pre_calculated_mean\n\n#     def __call__(self, tensor):\n#         return tensor - self.pre_calculated_mean\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# def show(list_imgs, cmap=cm.bone):\n#     list_slices = []\n#     for img3d in list_imgs:\n#         slc = int(img3d.shape[0] / 2)\n#         img = img3d[slc]\n#         list_slices.append(img)\n    \n#     fig, axs = plt.subplots(1, 5, figsize=(15, 7))\n#     for i, img in enumerate(list_slices):\n#         axs[i].imshow(img, cmap=cmap)\n#         axs[i].axis('off')\n        \n#     plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# # Show some images\n# test_img = CTScansDataset(\n#     root_dir=test_dir,\n#     transform=transforms.Compose([\n#         CropBoundingBox(),\n#         ConvertToHU(),\n#         Resize(resize_dims),\n#         Clip(bounds=clip_bounds),\n#         MaskWatershed(min_hu=min(clip_bounds), iterations=1, show_tqdm=True),\n#         Normalize(bounds=clip_bounds)\n#     ]))\n\n# list_imgs = [test_img[i]['image'] for i in range(len(test_img))]\n# show(list_imgs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# Only uncomment if pre-processing train data from scratch.\n# train_img = CTScansDataset(\n#     root_dir=train_dir,\n#     transform=transforms.Compose([\n#         CropBoundingBox(),\n#         ConvertToHU(),\n#         Resize(resize_dims),\n#         Clip(bounds=clip_bounds),\n#         MaskWatershed(min_hu=min(clip_bounds), iterations=1, show_tqdm=True),\n#         Normalize(bounds=clip_bounds)\n#     ]))\n\n# list_imgs = [train_img[i]['image'] for i in range(len(train_img))]\n# show(list_imgs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# class CTTensorsDataset(Dataset):\n#     def __init__(self, root_dir, transform=None):\n#         self.root_dir = Path(root_dir)\n#         self.tensor_files = sorted([f for f in self.root_dir.glob('*.pt')])\n#         self.transform = transform\n\n#     def __len__(self):\n#         return len(self.tensor_files)\n\n#     def __getitem__(self, item):\n#         if torch.is_tensor(item):\n#             item = item.tolist()\n\n#         image = torch.load(self.tensor_files[item])\n#         if self.transform:\n#             image = self.transform(image)\n\n#         return {\n#             'patient_id': self.tensor_files[item].stem,\n#             'image': image\n#         }\n\n#     def mean(self):\n#         cum = 0\n#         for i in range(len(self)):\n#             sample = self[i]['image']\n#             cum += torch.mean(sample).item()\n\n#         return cum / len(self)\n\n#     def random_split(self, val_size: float):\n#         num_val = int(val_size * len(self))\n#         num_train = len(self) - num_val\n#         return random_split(self, [num_train, num_val])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# torch.mean(data[0]['image']).item()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# train_img = CTTensorsDataset(\n#     root_dir=root_dir,\n#     transform=ZeroCenter(pre_calculated_mean=pre_calculated_mean)\n# )\n# cum = 0\n# for i in range(len(train_img)):\n#     sample = train_img[i]['image']\n#     cum += torch.mean(sample).item()\n\n# assert cum / len(train_img) == pytest.approx(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# class AutoEncoder(nn.Module):\n#     def __init__(self, latent_features=latent_features):\n#         super(AutoEncoder, self).__init__()\n#         # Encoder\n#         self.conv1 = nn.Conv3d(1, 16, 3)\n#         self.conv2 = nn.Conv3d(16, 32, 3)\n#         self.conv3 = nn.Conv3d(32, 96, 2)\n#         self.conv4 = nn.Conv3d(96, 1, 1)\n#         self.pool1 = nn.MaxPool3d(kernel_size=2, stride=2, return_indices=True)\n#         self.pool2 = nn.MaxPool3d(kernel_size=3, stride=3, return_indices=True)\n#         self.pool3 = nn.MaxPool3d(kernel_size=2, stride=2, return_indices=True)\n#         self.pool4 = nn.MaxPool3d(kernel_size=2, stride=2, return_indices=True)\n#         self.fc1 = nn.Linear(10 * 10, latent_features)\n#         # Decoder\n#         self.fc2 = nn.Linear(latent_features, 10 * 10)\n#         self.deconv0 = nn.ConvTranspose3d(1, 96, 1)\n#         self.deconv1 = nn.ConvTranspose3d(96, 32, 2)\n#         self.deconv2 = nn.ConvTranspose3d(32, 16, 3)\n#         self.deconv3 = nn.ConvTranspose3d(16, 1, 3)\n#         self.unpool0 = nn.MaxUnpool3d(kernel_size=2, stride=2)\n#         self.unpool1 = nn.MaxUnpool3d(kernel_size=2, stride=2)\n#         self.unpool2 = nn.MaxUnpool3d(kernel_size=3, stride=3)\n#         self.unpool3 = nn.MaxUnpool3d(kernel_size=2, stride=2)\n\n#     def encode(self, x, return_partials=True):\n#         # Encoder\n#         x = self.conv1(x)\n#         up3out_shape = x.shape\n#         x, i1 = self.pool1(x)\n\n#         x = self.conv2(x)\n#         up2out_shape = x.shape\n#         x, i2 = self.pool2(x)\n\n#         x = self.conv3(x)\n#         up1out_shape = x.shape\n#         x, i3 = self.pool3(x)\n\n#         x = self.conv4(x)\n#         up0out_shape = x.shape\n#         x, i4 = self.pool4(x)\n\n#         x = x.view(-1, 10 * 10)\n#         x = F.relu(self.fc1(x))\n\n#         if return_partials:\n#             return x, up3out_shape, i1, up2out_shape, i2, up1out_shape, i3, \\\n#                    up0out_shape, i4\n\n#         else:\n#             return x\n\n#     def forward(self, x):\n#         x, up3out_shape, i1, up2out_shape, i2, \\\n#         up1out_shape, i3, up0out_shape, i4 = self.encode(x)\n\n#         # Decoder\n#         x = F.relu(self.fc2(x))\n#         x = x.view(-1, 1, 1, 10, 10)\n#         x = self.unpool0(x, output_size=up0out_shape, indices=i4)\n#         x = self.deconv0(x)\n#         x = self.unpool1(x, output_size=up1out_shape, indices=i3)\n#         x = self.deconv1(x)\n#         x = self.unpool2(x, output_size=up2out_shape, indices=i2)\n#         x = self.deconv2(x)\n#         x = self.unpool3(x, output_size=up3out_shape, indices=i1)\n#         x = self.deconv3(x)\n\n#         return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# t0 = time()\n\n# # Load the data\n# data = CTTensorsDataset(\n#     root_dir=root_dir,\n#     transform=ZeroCenter(pre_calculated_mean=pre_calculated_mean)\n# )\n# # train_set, val_set = data.random_split(val_size)\n# train_set = data\n# val_set = None\n# # datasets = {'train': train_set, 'val': val_set}\n# datasets = {'train': train_set}\n# dataloaders = {\n#     x: DataLoader(\n#         datasets[x],\n#         batch_size=batch_size,\n# #         shuffle=(x == 'train'),\n#         shuffle=(False),\n#         num_workers=2\n# #     ) for x in ['train', 'val']}\n#     ) for x in ['train']}\n\n# # dataset_sizes = {x: len(datasets[x]) for x in ['train', 'val']}\n# dataset_sizes = {x: len(datasets[x]) for x in ['train']}\n\n# # Prepare for training\n# device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n# model = AutoEncoder(latent_features=latent_features).to(device)\n# criterion = torch.nn.MSELoss()\n# optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\n# best_model_wts = None\n# best_loss = np.inf\n\n# date_time = datetime.now().strftime(\"%Y%m%d-%H%M\")\n# log_dir = Path(tensorboard_dir) / f'{date_time}'\n# writer = SummaryWriter(log_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# # Training loop\n# for epoch in range(num_epochs):\n\n#     # Each epoch has a training and validation phase\n# #     for phase in ['train', 'val']:\n#     for phase in ['train']:\n#         if phase == 'train':\n#             model.train()  # Set model to training mode\n#         else:\n#             model.eval()   # Set model to evaluate mode\n\n#         running_loss = 0.0\n#         running_preds = 0\n\n#         # Iterate over data.\n#         bar = tqdm(dataloaders[phase])\n#         for inputs in bar:\n#             bar.set_description(f'Epoch {epoch + 1} {phase}'.ljust(20))\n#             inputs = inputs['image'].to(device, dtype=torch.float)\n\n#             # zero the parameter gradients\n#             optimizer.zero_grad()\n\n#             # forward\n#             # track history if only in train\n#             with torch.set_grad_enabled(phase == 'train'):\n#                 outputs = model(inputs)\n#                 loss = criterion(outputs, inputs)\n\n#                 # backward + optimize only if in training phase\n#                 if phase == 'train':\n#                     loss.backward()\n#                     optimizer.step()\n\n#             # statistics\n#             running_loss += loss.item() * inputs.size(0)\n#             running_preds += inputs.size(0)\n#             bar.set_postfix(loss=f'{running_loss / running_preds:0.6f}')\n\n#         epoch_loss = running_loss / dataset_sizes[phase]\n#         writer.add_scalar(f'Loss/{phase}', epoch_loss, epoch)\n\n#         # deep copy the model\n# #         if phase == 'val' and epoch_loss < best_loss:\n#         if phase == 'train' and epoch_loss < best_loss:\n#             best_loss = epoch_loss\n#             best_model_wts = copy.deepcopy(model.state_dict())\n#             torch.save(best_model_wts, model_file)\n\n# # load best model weights\n# model.load_state_dict(best_model_wts)\n\n# print(f'Done! Time {timedelta(seconds=time() - t0)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"excluded_patients = ['ID00011637202177653955184', 'ID00052637202186188008618', 'ID00102637202206574119190']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# slc = 0.5\n# sample_id = np.random.randint(len(data))\n# print(f'Inspecting CT Scan {data[sample_id][\"patient_id\"]}')\n\n# fig, axs = plt.subplots(1, 2, figsize=(10, 7))\n\n# sample = data[sample_id]['image'].squeeze(0).numpy()\n# axs[0].imshow(sample[int(40 * slc), :, :], cmap=cm.bone)\n# axs[0].axis('off')\n# imageio.mimsave(\"sample_input.gif\", sample, duration=0.1)\n\n# with torch.no_grad():\n#     img = data[sample_id]['image'].unsqueeze(0).float().to(device)\n#     latent_features = model.encode(img, return_partials=False)\\\n#         .squeeze().cpu().numpy().tolist()\n#     outputs = model(img).squeeze().cpu().numpy()\n\n# axs[1].imshow(outputs[int(40 * slc), :, :], cmap=cm.bone)\n# axs[1].axis('off')\n\n# imageio.mimsave(\"sample_output.gif\", outputs, duration=0.1)\n\n# rmse = ((sample - outputs)**2).mean()\n# plt.show()\n# print(f'Latent features: {latent_features} \\nLoss: {rmse}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"CT Scan preprocessing and feature generation code:"},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class CTScansDataset(Dataset):\n    def __init__(self, root_dir, transform=None, excl_patients=None):\n        self.root_dir = Path(root_dir)\n        self.patients = [p for p in self.root_dir.glob('*') if (p.is_dir() and p.name not in excl_patients)]\n        self.transform = transform\n        self.excl_patients = excl_patients\n\n    def __len__(self):\n        return len(self.patients)\n\n    def __getitem__(self, idx):\n        if torch.is_tensor(idx):\n            idx = idx.tolist()\n\n        image, metadata, msbs = self.load_scan(self.patients[idx])\n        sample = {'image': image, 'metadata': metadata, 'msbs': msbs}\n        if self.transform:\n            sample = self.transform(sample)\n\n        return sample\n\n    def save(self, path):\n        t0 = time()\n        Path(path).mkdir(exist_ok=True, parents=True)\n        print('Saving pre-processed dataset to disk')\n        sleep(1)\n        cum = 0\n        var = 0\n\n        bar = trange(len(self))\n        for i in bar:\n            sample = self[i]\n            image, data = sample['image'], sample['metadata']\n            cum += torch.mean(image).item()\n            var += torch.var(image).item()\n\n            bar.set_description(f'Saving CT scan {data.PatientID}')\n            fname = Path(path) / f'{data.PatientID}.pt'\n            torch.save(image, fname)\n\n        sleep(1)\n        bar.close()\n        print(f'Done! Time {timedelta(seconds=time() - t0)}\\n'\n              f'Mean value: {cum / len(self)}\\n'\n              f'Std value: {(var / len(self)) ** (1/2)}')\n\n    def get_patient(self, patient_id):\n        patient_ids = [str(p.stem) for p in self.patients]\n        return self.__getitem__(patient_ids.index(patient_id))\n\n    @staticmethod\n    def load_scan(path):\n        raw_slices = [pydicom.read_file(p) for p in path.glob('*.dcm')]\n        final_slices = [s for s in raw_slices\n                        if hasattr(s, 'ImagePositionPatient')]\n\n        if len(final_slices) > 0:\n            final_slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n            msbs = 0\n            for i in range(len(final_slices) - 1):\n                a = final_slices[i].ImagePositionPatient[2]\n                b = final_slices[i + 1].ImagePositionPatient[2]\n                msbs += abs(a - b)\n\n            msbs /= len(final_slices)\n            image = np.stack([s.pixel_array.astype(float) for s in final_slices])\n            metadata = final_slices[0]\n\n        else:\n            # Guess based on eyballing other mean spaces between slices:\n            msbs = None\n            if raw_slices[0].PatientID == 'ID00132637202222178761324':\n                msbs = 0.7\n            elif raw_slices[0].PatientID == 'ID00128637202219474716089':\n                msbs = 5.\n\n            warnings.warn(f'Patient {raw_slices[0].PatientID} CT scan does '\n                          f'not have \"ImagePositionPatient\". Assuming '\n                          f'filenames in the right scan order. Also, assuming'\n                          f'mean space between slices of {msbs}')\n\n            image = np.stack([s.pixel_array.astype(float) for s in raw_slices])\n            metadata = raw_slices[0]\n\n        return image, metadata, msbs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset = CTScansDataset(root_dir='../input/osic-pulmonary-fibrosis-progression/train', excl_patients=excluded_patients)\ndataset.patients","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class CropBoundingBox:\n    @staticmethod\n    def bounding_box(img3d: np.array):\n        mid_img = img3d[int(img3d.shape[0] / 2)]\n        same_first_row = (mid_img[0, :] == mid_img[0, 0]).all()\n        same_first_col = (mid_img[:, 0] == mid_img[0, 0]).all()\n        if same_first_col and same_first_row:\n            return True\n        else:\n            return False\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n        if not self.bounding_box(image):\n            return sample\n\n        mid_img = image[int(image.shape[0] / 2)]\n        r_min, r_max = None, None\n        c_min, c_max = None, None\n        for row in range(mid_img.shape[0]):\n            if not (mid_img[row, :] == mid_img[0, 0]).all() and r_min is None:\n                r_min = row\n            if (mid_img[row, :] == mid_img[0, 0]).all() and r_max is None \\\n                    and r_min is not None:\n                r_max = row\n                break\n\n        for col in range(mid_img.shape[1]):\n            if not (mid_img[:, col] == mid_img[0, 0]).all() and c_min is None:\n                c_min = col\n            if (mid_img[:, col] == mid_img[0, 0]).all() and c_max is None \\\n                    and c_min is not None:\n                c_max = col\n                break\n\n        image = image[:, r_min:r_max, c_min:c_max]\n        return {'image': image, 'metadata': data, 'msbs': sample['msbs']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class ConvertToHU:\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n\n        img_type = data.ImageType\n        is_hu = img_type[0] == 'ORIGINAL' and not (img_type[2] == 'LOCALIZER')\n        # if not is_hu:\n        #     warnings.warn(f'Patient {data.PatientID} CT Scan not cannot be'\n        #                   f'converted to Hounsfield Units (HU).')\n\n        intercept = data.RescaleIntercept\n        slope = data.RescaleSlope\n        image = (image * slope + intercept).astype(np.int16)\n        return {'image': image, 'metadata': data, 'msbs': sample['msbs']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class Resample:\n    def __init__(self, new_spacing=(1, 1, 1)):\n        \"\"\"new_spacing means now much every pixel represent in mm, in each\n        dimension. E.g. 2, 2, 2 means every pixel represent 2mm in  every\n        dimension.\n        \"\"\"\n        assert isinstance(new_spacing, tuple)\n        self.new_spacing = new_spacing\n\n    def __call__(self, sample):\n        image, data, msbs = sample['image'], sample['metadata'], sample['msbs']\n\n        spacing = np.array([msbs] + list(data.PixelSpacing), dtype=np.float32)\n        resize_factor = spacing / self.new_spacing\n        new_real_shape = image.shape * resize_factor\n        new_shape = np.round(new_real_shape)\n        real_resize_factor = new_shape / image.shape\n\n        image = zoom(image, real_resize_factor, mode='nearest')\n        return {'image': image, 'metadata': data}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class Clip:\n    def __init__(self, bounds=(-1000, 500)):\n        self.min = min(bounds)\n        self.max = max(bounds)\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n        image[image < self.min] = self.min\n        image[image > self.max] = self.max\n        return {'image': image, 'metadata': data}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class MaskWatershed:\n    def __init__(self, min_hu, iterations, show_tqdm):\n        self.min_hu = min_hu\n        self.iterations = iterations\n        self.show_tqdm = show_tqdm\n\n    def __call__(self, sample):\n        image, data = sample['image'], sample['metadata']\n\n        stack = []\n        if self.show_tqdm:\n            bar = trange(image.shape[0])\n            bar.set_description(f'Masking CT scan {data.PatientID}')\n        else:\n            bar = range(image.shape[0])\n        for slice_idx in bar:\n            sliced = image[slice_idx]\n            stack.append(self.seperate_lungs(sliced, self.min_hu,\n                                             self.iterations))\n\n        return {\n            'image': np.stack(stack),\n            'metadata': sample['metadata']\n        }\n\n    @staticmethod\n    def seperate_lungs(image, min_hu, iterations):\n        h, w = image.shape[0], image.shape[1]\n\n        marker_internal, marker_external, marker_watershed = MaskWatershed.generate_markers(image)\n\n        # Sobel-Gradient\n        sobel_filtered_dx = ndimage.sobel(image, 1)\n        sobel_filtered_dy = ndimage.sobel(image, 0)\n        sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)\n        sobel_gradient *= 255.0 / np.max(sobel_gradient)\n\n        watershed = morphology.watershed(sobel_gradient, marker_watershed)\n\n        outline = ndimage.morphological_gradient(watershed, size=(3,3))\n        outline = outline.astype(bool)\n\n        # Structuring element used for the filter\n        blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],\n                           [0, 1, 1, 1, 1, 1, 0],\n                           [1, 1, 1, 1, 1, 1, 1],\n                           [1, 1, 1, 1, 1, 1, 1],\n                           [1, 1, 1, 1, 1, 1, 1],\n                           [0, 1, 1, 1, 1, 1, 0],\n                           [0, 0, 1, 1, 1, 0, 0]]\n\n        blackhat_struct = ndimage.iterate_structure(blackhat_struct, iterations)\n\n        # Perform Black Top-hat filter\n        outline += ndimage.black_tophat(outline, structure=blackhat_struct)\n\n        lungfilter = np.bitwise_or(marker_internal, outline)\n        lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)\n\n        segmented = np.where(lungfilter == 1, image, min_hu * np.ones((h, w)))\n\n        return segmented\n\n    @staticmethod\n    def generate_markers(image, threshold=-400):\n        h, w = image.shape[0], image.shape[1]\n\n        marker_internal = image < threshold\n        marker_internal = segmentation.clear_border(marker_internal)\n        marker_internal_labels = measure.label(marker_internal)\n\n        areas = [r.area for r in measure.regionprops(marker_internal_labels)]\n        areas.sort()\n\n        if len(areas) > 2:\n            for region in measure.regionprops(marker_internal_labels):\n                if region.area < areas[-2]:\n                    for coordinates in region.coords:\n                        marker_internal_labels[coordinates[0], coordinates[1]] = 0\n\n        marker_internal = marker_internal_labels > 0\n\n        # Creation of the External Marker\n        external_a = ndimage.binary_dilation(marker_internal, iterations=10)\n        external_b = ndimage.binary_dilation(marker_internal, iterations=55)\n        marker_external = external_b ^ external_a\n\n        # Creation of the Watershed Marker\n        marker_watershed = np.zeros((h, w), dtype=np.int)\n        marker_watershed += marker_internal * 255\n        marker_watershed += marker_external * 128\n\n        return marker_internal, marker_external, marker_watershed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class LungVolume:\n    def __init__(self, unit, spacing, min_HU):\n        self.unit = unit\n        self.spacing = spacing\n        self.min_HU = min_HU\n\n    def __call__(self, sample):\n        image = sample['image']\n\n        volume = (image > self.min_HU).sum()\n        volumes_in_perc = volume / np.prod(image.shape)\n        volumes_in_litr = volumes_in_perc * np.prod(image.shape) * \\\n                          np.prod(self.spacing) / 10**6\n\n        if self.unit == '%':\n            sample['volume'] = volumes_in_perc\n        elif self.unit == 'l':\n            sample['volume'] = volumes_in_litr\n\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class ChestCircumference:\n    def __init__(self, spacing, min_HU):\n        self.spacing = spacing\n        self.min_HU = min_HU\n\n    def ellipsis_perimeter(self, img):\n        # trim slice, removing zero rows and cols\n        img = img[~(img == self.min_HU).all(1)].T\n        img = img[~(img == self.min_HU).all(1)]\n        # ellipsis perimeter approximation 2\n        # https://www.mathsisfun.com/geometry/ellipse-perimeter.html\n        a, b = img.shape\n        a /= 2\n        b /= 2\n        p = np.pi * (3 * (a + b) - np.sqrt((3 * a + b) * (a + 3 * b)))\n        return p\n\n    def __call__(self, sample):\n        image = sample['image']\n\n        perimeters = []\n        for slice_idx in range(image.shape[0]):\n            slc = image[slice_idx]\n            if slc.max() != self.min_HU:  # There's something in the slice\n                perimeters.append(self.ellipsis_perimeter(slc))\n            else:\n                perimeters.append(0)\n\n        # Transform in mm\n        cc_in_mm = max(perimeters) * self.spacing[1] * self.spacing[2]\n        sample['chest'] = int(cc_in_mm)\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class LungHeight:\n    def __init__(self, spacing):\n        self.spacing = spacing\n\n    def __call__(self, sample):\n        image = sample['image']\n        sample['height'] = image.shape[0] * self.spacing[0]\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"class HUHistogram:\n    def __init__(self, bounds, bins):\n        self.bounds = bounds\n        self.bins = bins\n\n    def __call__(self, sample):\n        image = sample['image']\n        pixels = image[image != min(self.bounds)]\n\n        if pixels.shape[0] > 0:\n            sample['mean'] = pixels.mean()\n            sample['std'] = pixels.std()\n            sample['skew'] = skew(pixels)\n            sample['kurtosis'] = kurtosis(pixels)\n\n            hist, _ = np.histogram(pixels, bins=self.bins, range=self.bounds)\n            sample['hist'] = hist / pixels.size  # Normalized\n            # Density = True in np.histogram doesn't work as expected\n        else:\n            sample['mean'] = 0\n            sample['std'] = 0\n            sample['skew'] = 0\n            sample['kurtosis'] = 0\n            sample['hist'] = np.zeros(self.bins)\n\n        return sample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def gennerate_row(dataset, i):\n    sample = dataset[i]\n    print(sample['metadata'].PatientID)\n    df = pd.DataFrame(columns=np.arange(8).tolist())\n    df.loc[0] = [sample['metadata'].PatientID, sample[\"volume\"],\n                 sample[\"chest\"], sample[\"height\"], sample[\"mean\"],\n                 sample[\"std\"], sample[\"skew\"], sample[\"kurtosis\"]]\n\n    dfhist = pd.DataFrame()\n    dfhist[0] = sample['hist']\n    dfhist = dfhist.T\n    df = pd.concat([df, dfhist], axis=1, ignore_index=True)\n\n    columns = ['Patient', 'Volume', 'Chest', 'Height', 'Mean',\n               'Std', 'Skew', 'Kurtosis']\n    columns += [f'bin{b + 1}' for b in range(len(dfhist.columns))]\n    df.columns = columns\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"t0 = time()\ndataset = CTScansDataset(\n    root_dir='../input/osic-pulmonary-fibrosis-progression/train',\n    transform=transforms.Compose([\n        CropBoundingBox(),\n        ConvertToHU(),\n        Resample(new_spacing=(1, 1, 1)),\n        Clip(bounds=(-1000, 400)),\n        MaskWatershed(min_hu=-1000, iterations=1, show_tqdm=False),\n        # Features\n        LungVolume(unit='l',\n                   spacing=(1, 1, 1),\n                   min_HU=-1000),\n        ChestCircumference(spacing=(1, 1, 1), min_HU=-1000),\n        LungHeight(spacing=(1, 1, 1)),\n        HUHistogram(bounds=(-1000, 400), bins=14)\n    ]),\n    excl_patients=excluded_patients\n)\n\nf = partial(gennerate_row, dataset)\nwith mp.Pool(processes=mp.cpu_count()) as pool:\n    results = list(tqdm(pool.imap_unordered(f, range(len(dataset))),\n                        total=len(dataset)))\n\nfeatures_train = pd.concat(results, ignore_index=True)\nfeatures_train.to_csv('features_train.csv', index=False)\n\nprint(f'\\nDone! Time: {timedelta(seconds=time() - t0)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"features_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"t0 = time()\ndataset = CTScansDataset(\n    root_dir='../input/osic-pulmonary-fibrosis-progression/test',\n    transform=transforms.Compose([\n        CropBoundingBox(),\n        ConvertToHU(),\n        Resample(new_spacing=(1, 1, 1)),\n        Clip(bounds=(-1000, 400)),\n        MaskWatershed(min_hu=-1000, iterations=1, show_tqdm=False),\n        # Features\n        LungVolume(unit='l',\n                   spacing=(1, 1, 1),\n                   min_HU=-1000),\n        ChestCircumference(spacing=(1, 1, 1), min_HU=-1000),\n        LungHeight(spacing=(1, 1, 1)),\n        HUHistogram(bounds=(-1000, 400), bins=14)\n    ]))\n\nf = partial(gennerate_row, dataset)\nwith mp.Pool(processes=mp.cpu_count()) as pool:\n    results = list(tqdm(pool.imap_unordered(f, range(len(dataset))),\n                        total=len(dataset)))\n\nfeatures_test = pd.concat(results, ignore_index=True)\nfeatures_test.to_csv('features_test.csv', index=False)\n\nprint(f'\\nDone! Time: {timedelta(seconds=time() - t0)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"features_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"features = pd.concat([features_train, features_test])\nfeatures.head(176)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Bayesian Code:"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"scrolled":false},"cell_type":"code","source":"# Preprocessing\n# train is resused/editted. train_raw is left alone\ntrain = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntrain_raw = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\n# train.drop(train[train.Patient == 'ID00197637202246865691526'].index, inplace=True)\n\n# exclude_test_patient_data_from_trainset = False\n\n# if exclude_test_patient_data_from_trainset:\n#     train = train[~train['Patient'].isin(test['Patient'].unique())]\n\ntrain = pd.concat([train, test], axis=0, ignore_index=True)\\\n    .drop_duplicates()\nle_id = LabelEncoder()\ntrain['PatientID'] = le_id.fit_transform(train['Patient'])\n\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def add_baselines(data):\n    aux = data[['Patient', 'Weeks']].groupby('Patient')\\\n        .min().reset_index()\n    aux = pd.merge(aux, data[['Patient', 'Weeks', 'FVC', 'Percent']], how='left', \n                   on=['Patient', 'Weeks'])\n    aux = aux.groupby('Patient').mean().reset_index()\n    aux['Weeks'] = aux['Weeks'].astype(int)\n    aux['FVC'] = aux['FVC'].astype(int)\n    aux['Percent'] = aux['Percent'].astype(float)\n    data = pd.merge(data, aux, how='left', on='Patient', suffixes=('', '_base'))\n    return data\n\ntrain = add_baselines(train)\ntest = add_baselines(test)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def patient_class(row):\n    if row['Sex'] == 'Male':\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 0\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 1\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 2\n    else:\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 3\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 4\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 5\n\ntrain['Class'] = train.apply(patient_class, axis=1)\ntest['Class'] = test.apply(patient_class, axis=1)\ntest.head()\ntrain.loc[train['Patient']=='ID00007637202177411956430']['Class'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"patients = pd.DataFrame()\n# Very simple pre-processing: creating patient indexes\n\n\npatients = train[['Patient', 'PatientID', 'Age', 'Class', 'Weeks_base', 'FVC_base', 'Percent_base']].drop_duplicates()\nfvc_data = train[['Patient', 'PatientID', 'Weeks', 'FVC']]\npatients['FVC_nom_base'] = patients['FVC_base']*(100./patients['Percent_base'])\npatients.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"features.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"train = train.merge(features, on=['Patient'])\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"test = test.merge(features, on=['Patient'])\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"patients = patients.merge(features, on=['Patient'])\npatients.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"patient_ids = patients['Patient']\n#patient_ids\n\nfor patient_id in patient_ids:\n    df_patient = pd.DataFrame(train.loc[train.Patient == patient_id])\n    x = np.array(train.loc[train.Patient == patient_id, 'Weeks'].values).reshape(-1,1)\n    y = train.loc[train.Patient == patient_id, 'FVC'].values\n    regr = linear_model.LinearRegression()\n    regr.fit(x,y)\n    patients.loc[patients.Patient==patient_id,'gradient'] = regr.coef_[0]\n    patients.loc[patients.Patient==patient_id,'intercept'] = regr.intercept_\n    patients.loc[patients.Patient==patient_id,'intercept'] = regr.intercept_\n    y_pred = regr.predict(x)\n    patients.loc[patients.Patient==patient_id,'rmse'] = np.sqrt(mean_squared_error(y, y_pred))\n\npatients.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# for i in range(len(data)):\n#     patient_id = data[i][\"patient_id\"]\n#     with torch.no_grad():\n#         img = data[i]['image'].unsqueeze(0).float().to(device)\n#         latent_features = model.encode(img, return_partials=False).squeeze().cpu().numpy().tolist()\n#         for j in range(len(latent_features)):\n#             patients.loc[patients.Patient==patient_id,'l'+str(j)] = latent_features[j]\n\n# patients.head(176)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# plt.title(\"l4\")\n# mean = patients['l4'].mean()\n# std = patients['l4'].std()\n# sns.distplot(patients['l4'], bins=100);\n# print(f'l4 mean: {mean}')\n# print(f'l4 std: {std}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# plt.title(\"Gradients\")\n# mean = patients['gradient'].mean()\n# std = patients['gradient'].std()\n# sns.distplot(patients['gradient'], bins=100);\n# print(f'Gradient mean: {mean}')\n# print(f'Gradient std: {std}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# plt.title(\"Intercepts\")\n# mean = patients['intercept'].mean()\n# std = patients['intercept'].std()\n# # sns.distplot(patients['intercept'], bins=100);\n# print(f'Intercept mean: {mean}')\n# print(f'Intercept std: {std}')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# plt.title(\"RMSEs\")\n# mean = patients['rmse'].mean()\n# std = patients['rmse'].std()\n# # sns.distplot(patients['rmse'], bins=100);\n# print(f'RMSE mean: {mean}')\n# print(f'RMSE std: {std}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"corr_matrix = patients.corr()\nsns.heatmap(corr_matrix)\ncorr_matrix","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"corr_matrix.sort_values(by=['gradient'])[corr_matrix['gradient'].abs()>0.1]['gradient']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"corr_matrix.sort_values(by=['intercept'])[corr_matrix['intercept'].abs()>0.48]['intercept']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"corr_matrix.sort_values(by=['rmse'])[corr_matrix['rmse'].abs()>0.15]['rmse']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"x = patients['Class']\ny = patients['intercept']\nplt.title(\"intercept\")\nplt.scatter(x,y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# x = patients['FVC_base']\n# y = patients['gradient']\n# fig, ax = plt.subplots(figsize=(10,10));\n# ax.set_title(\"gradient\")\n# sns.scatterplot(x,np.log(y),hue=patients['Class'], size=patients['Age'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# x = patients['Class']\n# y = patients['rmse']\n# plt.title(\"rmse\")\n# plt.scatter(x,y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"num_features = ['Age', 'Class', 'Weeks_base', 'FVC_base',\n       'Percent_base', 'FVC_nom_base', 'gradient', 'intercept']\ntarget = ['rmse']\n\nx = patients[num_features]\ny = patients[target]\nregr = linear_model.LinearRegression()\nregr.fit(x,y)\nplt.scatter(num_features, regr.coef_)\nplt.xticks(rotation=90)\nplt.title(\"RMSE\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# num_features = ['Age', 'Class', 'Weeks_base', 'FVC_base',\n#        'Percent_base', 'FVC_nom_base', 'gradient', 'rmse']\n# target = ['intercept']\n\n# x = patients[num_features]\n# y = patients[target]\n# regr = linear_model.LinearRegression()\n# regr.fit(x,y)\n# plt.scatter(num_features, regr.coef_)\n# plt.xticks(rotation=90)\n# plt.title(\"Intercept\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"# num_features = ['Age', 'Class', 'Weeks_base', 'FVC_base',\n#        'Percent_base', 'FVC_nom_base', 'intercept', 'rmse']\n# target = ['gradient']\n\n# x = patients[num_features]\n# y = patients[target]\n# regr = linear_model.LinearRegression()\n# regr.fit(x,y)\n# plt.scatter(num_features, regr.coef_)\n# plt.xticks(rotation=90)\n# plt.title(\"Gradient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"PatientID = train['Patient'].values\nfvc_b = train.groupby('Patient').first()['FVC_base']\n#fvc_b.values\nage = train.groupby('PatientID').first()['Age']\nage","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def model_fit(data, examine=True):\n    n_patients = data['Patient'].nunique()\n    FVC_obs = data['FVC'].values\n    Weeks = data['Weeks'].values\n    #X = train[['Weeks', 'Male', 'ExSmoker', 'CurrentlySmokes', 'Percent_base']].values\n    PatientID = data['PatientID'].values\n    patient_class = data['Class'].values\n    FVC_b = data.groupby('PatientID').first()['FVC_base']\n    w_b = data.groupby('PatientID').first()['Weeks_base']\n    age = data.groupby('PatientID').first()['Age']\n    # Small denotes 1 per patient, ie. not whole array of class values from all datapoints\n    patient_class_small = data.groupby('PatientID').first()['Class']\n\n    with pm.Model() as model:\n        # create shared variables that can be changed later on\n        FVC_obs_shared = pm.Data(\"FVC_obs_shared\", FVC_obs)\n        Weeks_shared = pm.Data('Weeks_shared', Weeks)\n        PatientID_shared = pm.Data('PatientID_shared', PatientID)\n        patient_class_shared = pm.Data('patient_class_shared', patient_class)\n        FVC_b_shared = pm.Data('FVC_b_shared', FVC_b)\n        w_b_shared= pm.Data('w_b_shared', w_b)\n        age_shared = pm.Data('age_shared', age)\n        patient_class_small_shared = pm.Data('patient_class_small_shared', patient_class_small)\n\n        #mu_a_base = pm.Normal('mu_a_base', mu=0., sigma=100)\n        #mu_a = FVC_b_shared + mu_a_base * w_b_shared\n        mu_a = pm.Normal('mu_a', mu=1700., sigma=400)\n        sigma_a = pm.HalfNormal('sigma_a', 1000.)\n        \n        \n        alpha_s = pm.Normal('alpha_s', -7, sigma=2)\n        beta_cs = pm.Normal('beta_cs', 0, sigma=2.5, shape=6)\n        mu_b = alpha_s + beta_cs[patient_class_small_shared]*age_shared \n        sigma_b = pm.HalfNormal('sigma_b', 5.)\n        \n#         mu_b = pm.Normal('mu_b', mu=-4., sigma=1)\n#         sigma_b = pm.HalfNormal('sigma_b', 5.)\n        \n        #mu_c = pm.Normal('mu_c', mu=0, sigma=100)\n        #sigma_c = pm.HalfNormal('sigma_c', 100)\n\n        a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=n_patients)\n        b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=n_patients)\n        #c = pm.Normal('c', mu=mu_c, sigma=sigma_c, shape=n_patients)\n\n        # Model error\n        #sigma = pm.HalfNormal('sigma', 150.)\n        sigma = pm.HalfNormal('sigma', 150., shape=6)\n        \n        FVC_est = a[PatientID_shared] + b[PatientID_shared] * Weeks_shared\n\n        # Data likelihood\n        FVC_like = pm.Normal('FVC_like', mu=FVC_est,\n                             sigma=sigma[patient_class_shared], observed=FVC_obs_shared)\n\n        # Fitting the model\n        trace = pm.sample(1000, tune=2000, target_accept=.9, init=\"adapt_diag\")\n        \n    if examine:\n        with model:\n            pm.traceplot(trace);\n            \n    return model, trace","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# with model_a:\n#     pm.traceplot(trace_a);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def generate_template(data):\n    pred_template = []\n    for i, patient in enumerate(data['Patient'].unique()):\n        df = pd.DataFrame(columns=['PatientID', 'Weeks', 'Patient', 'Class'])\n        df['Weeks'] = np.arange(-12, 134)\n        df['Patient'] = patient\n        df['Class'] = data.loc[data['Patient']==patient]['Class'].max()\n        pred_template.append(df)\n    pred_template = pd.concat(pred_template, ignore_index=True)\n    pred_template['PatientID'] = le_id.transform(pred_template['Patient'])\n    return pred_template","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"template_train_test = generate_template(test)\ntemplate_train_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"def model_predict(model, trace, template):\n    with model:\n        pm.set_data({\n            \"PatientID_shared\": template['PatientID'].values.astype(int),\n            \"Weeks_shared\": template['Weeks'].values.astype(int),\n            \"FVC_obs_shared\": np.zeros(len(template)).astype(int),\n            \"patient_class_shared\": template['Class'].values.astype(int),\n        })\n        post_pred = pm.sample_posterior_predictive(trace)\n    df = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC_pred', 'sigma'])\n    df['Patient'] = le_id.inverse_transform(template['PatientID'])\n    df['Weeks'] = template['Weeks']\n    df['FVC_pred'] = post_pred['FVC_like'].T.mean(axis=1)\n    df['sigma'] = (post_pred['FVC_like'].T.std(axis=1))\n    df['FVC_inf'] = df['FVC_pred'] - df['sigma']\n    df['FVC_sup'] = df['FVC_pred'] + df['sigma']\n    df = pd.merge(df, train[['Patient', 'Weeks', 'FVC']], how='left', on=['Patient', 'Weeks'])\n    #if exclude_test_patient_data_from_trainset:\n    #df=df.append(train_raw[train_raw['Patient'].isin(test['Patient'].unique())][['Patient', 'Weeks', 'FVC']])\n    df = df.rename(columns={'FVC': 'FVC_true'})\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def examine_predictions(data):\n    n = (data['Patient'].nunique())+1\n    f, axes = plt.subplots((n//3)+1, 3, figsize=(15, 5*((n//3)+1)))\n    #for i, patient in enumerate(np.random.choice(data['Patient'].unique(), size=3, replace=False)):\n    for i, patient in enumerate(data['Patient'].unique()):\n        ax = axes[i//3, i%3]\n        df = data[data['Patient'] == patient]\n        x = df['Weeks']\n        ax.set_title(patient)\n        ax.plot(x, df['FVC_true'], 'o')\n        ax.plot(x, df['FVC_pred'])\n        ax = sns.regplot(x, df['FVC_true'], ax=ax, ci=None, line_kws={'color':'red'})\n        ax.fill_between(x, df[\"FVC_inf\"], df[\"FVC_sup\"],alpha=0.5, color='#ffcd3c')\n        ax.set_ylabel('FVC')\n    axes[n//3,n%3].plot()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def evaluate_predictions(df, use_only_last_3_measures=False, examine=False):\n    if use_only_last_3_measures:\n        y = df.dropna().groupby('Patient').tail(3)\n    else:\n        y = df.dropna()\n\n    rmse = ((y['FVC_pred'] - y['FVC_true']) ** 2).mean() ** (1/2)\n    mae = (y['FVC_pred'] - y['FVC_true']).abs()\n    mae_mean = (np.sqrt((y['FVC_pred'] - y['FVC_true']) ** 2)).mean()\n    mae_sd = (np.sqrt((y['FVC_pred'] - y['FVC_true']) ** 2)).mean()\n    mae_max = (np.sqrt((y['FVC_pred'] - y['FVC_true']) ** 2)).mean()\n    #print(f'RMSE (mean): {rmse:.1f} ml')\n    #print(f'MAE (mean): {rmse:.1f} ml')\n    sigma_c = y['sigma'].values\n    sigma_c[sigma_c < 70] = 70\n    delta = (y['FVC_pred'] - y['FVC_true']).abs()\n    delta[delta > 1000] = 1000\n    lll = - np.sqrt(2) * delta / sigma_c - np.log(np.sqrt(2) * sigma_c)\n\n    #print(f'Laplace Log Likelihood: {lll.mean():.4f}')\n    y['sigma_c'] = y['sigma']\n    y['sigma_c'].values[y['sigma_c'].values < 70] = 70\n    y['delta_c'] = (y['FVC_pred'] - y['FVC_true']).abs()\n    y['delta_c'].values[y['delta_c'].values > 1000] = 1000\n    y['main_loss'] = y['delta_c']/y['sigma_c']\n    y['lll'] = - np.sqrt(2) * y['delta_c'] / y['sigma_c'] - np.log(np.sqrt(2) * y['sigma_c'])\n    patient_llls = y.groupby('Patient')['lll'].mean()\n    if examine:\n        #plt.hist(y['main_loss'], bins=100)\n        plt.hist(patient_llls, bins=100);\n    #bad_patients = y[y['main_loss'] > 2.5]\n    #bad_patients\n    #bad_patients['Patient'].unique()\n\n    #plt.scatter(delta, sigma_c);\n    #plt.xlabel(\"delta\");\n    #plt.ylabel(\"sigma_c\");\n    return lll.mean(), patient_llls","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def evaluation_cycle(train, valid, examine_trace=True, examine_preds=True):\n    # Fit model\n    print(\"Fit model ...\")\n    model,trace = model_fit(train, examine=examine_trace)\n\n    # Model predictions with training data\n    print(\"Examine true vs predictions for training data ...\")\n    template_train = generate_template(train)\n    template_train.head()\n    pred_train = model_predict(model,trace,template_train)\n    if examine_preds:\n        examine_predictions(pred_train)\n    lll_train, train_patient_llls = evaluate_predictions(pred_train)\n\n    # Model predictions with validation data (unseen)\n    if valid is not None:\n        print(\"Examine true vs predictions for validation data ...\")\n        template_valid = generate_template(valid)\n        pred_valid = model_predict(model,trace,template_valid)\n        if examine_preds:\n            examine_predictions(pred_valid)\n        lll_valid, valid_patient_llls = evaluate_predictions(pred_valid)\n    return pred_train, pred_valid, lll_train, lll_valid, train_patient_llls, valid_patient_llls","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"#Evaluate model\n\nlll_trains = []\nlll_valids = []\n\nfold = 0\nkfold = KFold(3, shuffle=True, random_state=1)\nall_patients = train['Patient'].unique()\n# Fold\n# for train_index, validation_index in kfold.split(all_patients):\n#     validation_patients = all_patients[validation_index]\n#     fold += 1\n\n# Random\nfor fold in range(1):\n    validation_patients = np.random.choice(all_patients, size=30, replace=False)\n    \n    examine_lll = True\n    \n    df_valid = train[train['Patient'].isin(validation_patients)]\n    df_train = train[~train['Patient'].isin(validation_patients)]\n\n    df_valid_first_readings = df_valid.groupby('Patient').head(1)\n    df_train = pd.concat([df_train, df_valid_first_readings], axis=0, ignore_index=True)\n\n    print(f'Fold: {fold}')\n    pred_train, pred_valid, lll_train, lll_valid, train_patient_llls, valid_patient_llls = evaluation_cycle(df_train, df_valid, examine_trace=True, examine_preds=False)\n\n    print(f'Laplace Log Likelihoods for fold: {fold}')\n    print(f'Training:     {lll_train:.4f}')\n    print(f'Validation:   {lll_valid:.4f}')\n    print(\"\")\n        \n    train_worst = train_patient_llls.nsmallest(10).index\n    valid_worst = valid_patient_llls.nsmallest(10).index\n    \n    print(f'Training - worst:   {*train_worst,}')\n    print(f'Validation - worst: {*valid_worst,}')\n    print(\"\")\n    pred_train_worst = pred_train[pred_train['Patient'].isin(train_worst)]\n    pred_valid_worst = pred_valid[pred_valid['Patient'].isin(valid_worst)]\n    \n    examine_predictions(pred_train_worst)\n    examine_predictions(pred_valid_worst)\n    \n    lll_trains.append(lll_train)\n    lll_valids.append(lll_valid)\n    \n    evaluate_predictions(pred_train, use_only_last_3_measures=True, examine=examine_lll)\n    evaluate_predictions(pred_valid, use_only_last_3_measures=True, examine=examine_lll)\n\nprint(f'Training LLLs:   {*lll_trains,}')\nprint(f'Validation LLLs: {*lll_valids,}')\nprint(\"\")\n\nlll_trains_array = np.array(lll_trains)\nlll_valids_array = np.array(lll_valids)\n\nprint(f'Training LLLs Mean:   {lll_trains_array.mean():.4f}')\nprint(f'Training LLLs STD:    {lll_trains_array.std():.4f}')\nprint(\"\")\nprint(f'Validation LLLs Mean: {lll_valids_array.mean():.4f}')\nprint(f'Validation LLLs STD:  {lll_valids_array.std():.4f}')\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# Fit model using all of training data\nprint(\"Fit model ...\")\nmodel,trace = model_fit(train, examine=False)\nprint(\"\")\n\n# Make predictions for test data\nprint(\"Make predictions for test data ...\")\ntemplate_test = generate_template(test)\ntemplate_test.head()\npred_test = model_predict(model,trace,template_test)\n#examine_predictions(pred_test)\n\n# Prepare final submission\nfinal = pd.DataFrame(columns=['Patient_Week', 'FVC', 'Confidence'])\nfinal['Patient_Week'] = pred_test['Patient'] + '_' + pred_test['Weeks'].astype(str)\nfinal['FVC'] = pred_test['FVC_pred']\nfinal['Confidence'] = pred_test['sigma']\nfinal.head()\nfinal.to_csv('submission.csv', index=False)\nprint(final.shape)\nfinal.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# pred_template = []\n# for p in test['Patient'].unique():\n#     df = pd.DataFrame(columns=['PatientID', 'Weeks'])\n#     df['Weeks'] = np.arange(-12, 134)\n#     df['Patient'] = p\n#     pred_template.append(df)\n# pred_template = pd.concat(pred_template, ignore_index=True)\n# pred_template['PatientID'] = le_id.transform(pred_template['Patient'])\n\n# with model_a:\n#     pm.set_data({\n#         \"PatientID_shared\": pred_template['PatientID'].values.astype(int),\n#         \"Weeks_shared\": pred_template['Weeks'].values.astype(int),\n#         \"FVC_obs_shared\": np.zeros(len(pred_template)).astype(int),\n#     })\n#     post_pred = pm.sample_posterior_predictive(trace_a)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"# df = pd.DataFrame(columns=['Patient', 'Weeks', 'Patient_Week', 'FVC', 'Confidence'])\n# df['Patient'] = pred_template['Patient']\n# df['Weeks'] = pred_template['Weeks']\n# df['Patient_Week'] = df['Patient'] + '_' + df['Weeks'].astype(str)\n# df['FVC'] = post_pred['FVC_like'].T.mean(axis=1)\n# df['Confidence'] = post_pred['FVC_like'].T.std(axis=1)\n# final = df[['Patient_Week', 'FVC', 'Confidence']]\n# final.to_csv('submission.csv', index=False)\n# print(final.shape)\n# final.head()","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}