{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pydicom\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold\nimport cv2\nimport time \nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=1000000007):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\nprint(time.time())\nseed_everything(int(time.time()))\nROOT = \"../input/osic-pulmonary-fibrosis-progression\"\nBATCH_SIZE=128\ntr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nchunk = pd.read_csv(f\"{ROOT}/test.csv\")\n\nsub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")\n\ntr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\nAll_data = tr.append([chunk, sub])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"All_data['min_week'] = All_data['Weeks']\nAll_data.loc[All_data.WHERE=='test','min_week'] = np.nan\nAll_data['min_week'] = All_data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = All_data.loc[All_data.Weeks == All_data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"All_data = All_data.merge(base, on='Patient', how='left')\nAll_data['base_week'] = All_data['Weeks'] - All_data['min_week']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"COLS = ['Sex','SmokingStatus'] #,'Age'\nFE = []\nfor col in COLS:\n    for mod in All_data[col].unique():\n        FE.append(mod)\n        All_data[mod] = (All_data[col] == mod).astype(int)\n#=================","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def csv_normlize(data,attr1,attr):\n    data[attr1] = (data[attr] - data[attr].min() ) / ( data[attr].max() - data[attr].min() )\ncsv_normlize(All_data,'age','Age')\ncsv_normlize(All_data,'BASE','min_FVC')\ncsv_normlize(All_data,'week','base_week')\ncsv_normlize(All_data,'percent','Percent')\n\nFE += ['age','percent','week','BASE']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import copy\nfrom datetime import timedelta, datetime\nimport imageio\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport multiprocessing\nimport numpy as np\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\nimport warnings\nimport scipy\n\nroot_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/train'\ntest_dir = '/kaggle/input/osic-pulmonary-fibrosis-progression/test'\ndesk_dir = 'save_output'\nmodel_file = '/kaggle/working/diophantus.pt'\nresize_dims = (40, 256, 256)\nclip_bounds = (-1000, 200)\nwatershed_iterations = 1\npre_calculated_mean = 0.02865046213070556\nlatent_features = 3\nbatch_size = 16\nlearning_rate = 1e-4\nnum_epochs = 3\nval_size = 0.2\ntensorboard_dir = '/kaggle/working/runs'\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":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            try:\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            except:\n                print(i)\n                pass\n        sleep(1)\n        bar.close()\n        print(f'Done! Time {timedelta(seconds=time() - t0)}\\n'\n              f'Tot value: {cum }')\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]\nclass 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}\nclass 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\nclass 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'''\nclass 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        z_spa = data.SliceThickness\n        pix_spa = [float(i) for i in list(data.PixelSpacing)]\n        spacing = np.array([z_spa]+pix_spa, dtype=np.float32)\n        \n        resize_factor = spacing / self.output_size\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        image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n        return {'image': image, 'metadata': data}\n'''\n\nclass 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}\nclass 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\nclass 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\nclass 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    \nclass 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":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()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''\ntest = 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\nlist_imgs = [test[i]['image'] for i in range(len(test))]\nshow(list_imgs)\n'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = 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(\n            min_hu=min(clip_bounds),\n            iterations=watershed_iterations,\n            show_tqdm=False),\n        Normalize(bounds=clip_bounds),\n        ToTensor()\n    ]))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":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},"cell_type":"code","source":"train = CTTensorsDataset(\n    root_dir='../input/osic-cached-dataset',\n    transform=ZeroCenter(pre_calculated_mean=pre_calculated_mean)\n)\ncum = 0\nfor i in range(len(train)):\n    sample = train[i]['image']\n    cum += torch.mean(sample).item()\n\nassert cum / len(train) == pytest.approx(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":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},"cell_type":"code","source":"t0 = time()\n\n# Load the data\ndata = CTTensorsDataset(\n    root_dir='../input/osic-cached-dataset',\n    transform=ZeroCenter(pre_calculated_mean=pre_calculated_mean)\n)\ntrain_set, val_set = data.random_split(val_size)\ndatasets = {'train': train_set, 'val': val_set}\ndataloaders = {\n    x: DataLoader(\n        datasets[x],\n        batch_size=batch_size,\n        shuffle=(x == 'train'),\n        num_workers=2\n    ) for x in ['train', 'val']}\n\ndataset_sizes = {x: len(datasets[x]) for x in ['train', 'val']}\n\n# Prepare for training\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\nmodel = AutoEncoder(latent_features=latent_features).to(device)\ncriterion = torch.nn.MSELoss()\noptimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\nbest_model_wts = None\nbest_loss = np.inf\n\ndate_time = datetime.now().strftime(\"%Y%m%d-%H%M\")\nlog_dir = Path(tensorboard_dir) / f'{date_time}'\nwriter = SummaryWriter(log_dir)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfor epoch in range(num_epochs):\n\n    # Each epoch has a training and validation phase\n    for phase in ['train', 'val']:\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            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\nmodel.load_state_dict(best_model_wts)\n\nprint(f'Done! Time {timedelta(seconds=time() - t0)}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slc = 0.5\nsample_id = np.random.randint(len(data))\nprint(f'Inspecting CT Scan {data[sample_id][\"patient_id\"]}')\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 7))\n\nsample = data[sample_id]['image'].squeeze(0).numpy()\naxs[0].imshow(sample[int(40 * slc), :, :], cmap=cm.bone)\naxs[0].axis('off')\nimageio.mimsave(\"sample_input.gif\", sample, duration=0.0001)\n\nwith 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\naxs[1].imshow(outputs[int(40 * slc), :, :], cmap=cm.bone)\naxs[1].axis('off')\n\nimageio.mimsave(\"sample_output.gif\", outputs, duration=0.0001)\n\nrmse = ((sample - outputs)**2).mean()\nplt.show()\nprint(f'Latent features: {latent_features} \\nLoss: {rmse}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"latent_Dicom={}\nfor i in range(len(test)):\n    Pid = test[i]['metadata'].PatientID\n   \n    print(f'Inspecting CT Scan {Pid}')\n    sample = test[i]['image'].squeeze(0).numpy()\n    with torch.no_grad():\n        img = test[i]['image'].unsqueeze(0).float().to(device)\n        latent_features = model.encode(img, return_partials=False)\\\n            .squeeze().cpu().numpy().tolist()\n        latent_Dicom[Pid]=latent_features\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(len(data)):\n    Pid=data[i][\"patient_id\"]\n    \n    sample = data[i]['image'].squeeze(0).numpy()\n\n    with torch.no_grad():\n        img = data[i]['image'].unsqueeze(0).float().to(device)\n        latent_features = model.encode(img, return_partials=False)\\\n            .squeeze().cpu().numpy().tolist()\n        latent_Dicom[Pid]=latent_features\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"latent_features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"latent_features_tot=3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pl = All_data['Patient'].tolist()\npl2=[latent_Dicom[i] if i in latent_Dicom else np.nan for i in pl]\npl3=np.array(pl2)\nfor i in range(latent_features_tot):\n    stri = 'latent_features'+str(i)\n    pl4=pl3[:,i]\n    All_data[stri]=pl4\n    FE.append(stri)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = All_data.loc[All_data.WHERE=='train']\nchunk = All_data.loc[All_data.WHERE=='val']\nsub = All_data.loc[All_data.WHERE=='test']\ntr.shape, chunk.shape, sub.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return K.mean(metric)\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    qs = [0.2, 0.50, 0.8]\n    q = tf.constant(np.array([qs]), dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return K.mean(v)\n#=============================#\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss\n#=================\ndef make_model(nh):\n    z = L.Input((nh,), name=\"Patient\")\n    x = L.Dense(100, activation=\"relu\", name=\"d1\")(z)\n    x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n    p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n    \n    model = M.Model(z, preds, name=\"CNN\")\n    model.compile(loss=mloss(0.8), optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = tr['FVC'].values\nz = tr[FE].values\nze = sub[FE].values\nnh = z.shape[1]\npe = np.zeros((ze.shape[0], 3))\npred = np.zeros((z.shape[0], 3))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(len(tr))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = make_model(nh)\nprint(net.summary())\nprint(net.count_params())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD = 5\nkf = KFold(n_splits=NFOLD)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def score2(y_true, y_pred):\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    \n    #sigma_clip = sigma + C1\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = np.abs(y_true - fvc_pred)\n    delta[delta>1000]=1000\n    sq2 = np.sqrt(2)\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return np.mean(metric)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ncnt = 0\nEPOCHS = 800\nfor tr_idx, val_idx in kf.split(z):\n    cnt += 1\n    print(f\"FOLD {cnt}\")\n    net = make_model(nh)\n    net.fit(z[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, \n            validation_data=(z[val_idx], y[val_idx]), verbose=0) #\n    print(\"train\", net.evaluate(z[tr_idx], y[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"val\", net.evaluate(z[val_idx], y[val_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"predict val...\")\n    pred[val_idx] = net.predict(z[val_idx], batch_size=BATCH_SIZE, verbose=0)\n    print(\"predict test...\")\n    pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD\n#==============","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nprint(score2(y,pred))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(unc.min(), unc.mean(), unc.max(), (unc>=0).mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(unc)\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pe","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub2=sub.copy(deep=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub2['FVC'] = pe[:, 1]\nsub2['Confidence'] = pe[:, 2] - pe[:, 0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"otest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nfor i in range(len(otest)):\n    sub2.loc[sub2['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'FVC'] = otest.FVC[i]\n    sub2.loc[sub2['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'Confidence'] = 0.1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub2[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","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}