{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# EDA","metadata":{}},{"cell_type":"code","source":"# imports\nimport numpy as np\nimport pandas as pd\nimport glob\nfrom pydicom.data import get_testdata_file\nfrom pydicom import dcmread\nimport matplotlib.pyplot as plt\nfrom tqdm import trange","metadata":{"execution":{"iopub.status.busy":"2021-08-21T15:24:31.768102Z","iopub.execute_input":"2021-08-21T15:24:31.768517Z","iopub.status.idle":"2021-08-21T15:24:32.097821Z","shell.execute_reply.started":"2021-08-21T15:24:31.768423Z","shell.execute_reply":"2021-08-21T15:24:32.096838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# read in data\ntrain_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv')\ntest_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:34.007472Z","iopub.execute_input":"2021-08-16T23:51:34.007909Z","iopub.status.idle":"2021-08-16T23:51:34.038785Z","shell.execute_reply.started":"2021-08-16T23:51:34.007868Z","shell.execute_reply":"2021-08-16T23:51:34.037939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find total number of unique patients\ntrain_df['Patient'].nunique()","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:34.676685Z","iopub.execute_input":"2021-08-16T23:51:34.677245Z","iopub.status.idle":"2021-08-16T23:51:34.696014Z","shell.execute_reply.started":"2021-08-16T23:51:34.677211Z","shell.execute_reply":"2021-08-16T23:51:34.695096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# finding number of patients that are male and female\nprint(test_df['Sex'].unique())\nprint(train_df[train_df['Sex'] == 'Female'].shape)\nprint(train_df[train_df['Sex'] == 'Male'].shape)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:35.179195Z","iopub.execute_input":"2021-08-16T23:51:35.17956Z","iopub.status.idle":"2021-08-16T23:51:35.188157Z","shell.execute_reply.started":"2021-08-16T23:51:35.179527Z","shell.execute_reply":"2021-08-16T23:51:35.186844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# getting list of paths to all subjects folders\nsubject_paths = glob.glob('/kaggle/input/osic-pulmonary-fibrosis-progression/train/*')","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:36.39003Z","iopub.execute_input":"2021-08-16T23:51:36.390673Z","iopub.status.idle":"2021-08-16T23:51:36.438243Z","shell.execute_reply.started":"2021-08-16T23:51:36.390619Z","shell.execute_reply":"2021-08-16T23:51:36.437305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# getting list of all dicom paths\ndicom_paths = []\nfor subject in subject_paths:\n    subject_dicoms = glob.glob(subject + '/*')\n    dicom_paths.append(subject_dicoms)","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:36.911906Z","iopub.execute_input":"2021-08-16T23:51:36.912427Z","iopub.status.idle":"2021-08-16T23:51:46.642333Z","shell.execute_reply.started":"2021-08-16T23:51:36.912389Z","shell.execute_reply":"2021-08-16T23:51:46.641093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# reading in a particular dicom path\nds = dcmread(dicom_paths[100][13])","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:46.64488Z","iopub.execute_input":"2021-08-16T23:51:46.64522Z","iopub.status.idle":"2021-08-16T23:51:46.696983Z","shell.execute_reply.started":"2021-08-16T23:51:46.645188Z","shell.execute_reply":"2021-08-16T23:51:46.695876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# displaying the pixel array in the dicom\nplt.figure(figsize = (7, 7))\nplt.imshow(ds.pixel_array, cmap=\"plasma\")\nplt.axis('off');","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:46.699474Z","iopub.execute_input":"2021-08-16T23:51:46.699988Z","iopub.status.idle":"2021-08-16T23:51:46.935544Z","shell.execute_reply.started":"2021-08-16T23:51:46.699923Z","shell.execute_reply":"2021-08-16T23:51:46.934048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find max value in the pixel array\nds.pixel_array.max()","metadata":{"execution":{"iopub.status.busy":"2021-08-16T23:51:51.238057Z","iopub.execute_input":"2021-08-16T23:51:51.23845Z","iopub.status.idle":"2021-08-16T23:51:51.245433Z","shell.execute_reply.started":"2021-08-16T23:51:51.238416Z","shell.execute_reply":"2021-08-16T23:51:51.244578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#  Model Training","metadata":{}},{"cell_type":"code","source":"# checking which GPU was allocated\n!nvidia-smi","metadata":{"execution":{"iopub.status.busy":"2021-08-21T15:24:39.650819Z","iopub.execute_input":"2021-08-21T15:24:39.651193Z","iopub.status.idle":"2021-08-21T15:24:40.417973Z","shell.execute_reply.started":"2021-08-21T15:24:39.651163Z","shell.execute_reply":"2021-08-21T15:24:40.416812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Libraries Utilized\n* numpy - used for arrays and fast random selection\n* pandas - used for reading in csv file and accessing row data for generating model inputs and targets\n* glob - used for finding all dicom files belonging to a particular patient (finds dicom files in patient folder)\n* matplotlib - used for visualizing training and validation loss curves\n* pydicom - used as an interface to read in dicom information and pixel values\n* cv2 - used to resize images to a smaller and consistent shape\n* scikit-learn - used for ```train_test_split``` and for r^2 and RMSE calculation\n* rich - used for a progress bar to visualize model training runtime\n* pytorch - used for model creation, training, and data preparation\n* efficientnet_pytorch - used for pretrained EfficientNet models for use in transfer learning","metadata":{}},{"cell_type":"code","source":"!pip install rich efficientnet_pytorch","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:21.524054Z","iopub.execute_input":"2021-08-29T21:28:21.524403Z","iopub.status.idle":"2021-08-29T21:28:23.576475Z","shell.execute_reply.started":"2021-08-29T21:28:21.524370Z","shell.execute_reply":"2021-08-29T21:28:23.575513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# imports\nimport math\nimport numpy as np\nimport pandas as pd\nimport glob\nimport matplotlib.pyplot as plt\nfrom pydicom import dcmread\nimport cv2\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom rich.progress import track","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:23.580070Z","iopub.execute_input":"2021-08-29T21:28:23.580362Z","iopub.status.idle":"2021-08-29T21:28:24.783085Z","shell.execute_reply.started":"2021-08-29T21:28:23.580333Z","shell.execute_reply":"2021-08-29T21:28:24.781168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Torch-specific imports\nimport torch\nfrom torch import nn\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision import transforms\nfrom efficientnet_pytorch import EfficientNet","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:24.783837Z","iopub.status.idle":"2021-08-29T21:28:24.784385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to read in dicom pixel values from file, and to reshape and resize to smaller size - 224x224\ndef get_img(path):\n    # read in file with pydicom\n    d = dcmread(path)\n    # rescale pixel intensities to match actual values\n    pxls = (d.pixel_array - d.RescaleIntercept) / (d.RescaleSlope * 1000)\n    # resize pixel intensities to a 224 by 224 array\n    return cv2.resize(pxls, (224, 224))","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:24.788823Z","iopub.status.idle":"2021-08-29T21:28:24.789217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PyTorch Dataset Class\n### ```__init__``` function:\n* Variables:\n    * ```self.mode``` to tell function where to find dicom images for particular subject (train or test folder)\n    * ```self.df``` specifies the dataframe containing train, validation, or test patients/FVC measurements\n    * ```self.transforms``` specifies particular transforms applied on image - here, the image is only converted to a pytorch tensor\n\n### ```__len__``` function:\n * Function that tells PyTorch the size of the dataset being used\n * set to ```len(self.df)``` because the dataframe contains the total number of pids and each image/target is selected once per patient\n\n### ```__getitem__``` function:\n* particular data point at specified idx selected\n* ```pid``` stores Patient ID for the particular datapoint\n* ```label``` stores FVC measurement for a particular Patient ID\n* ```img``` contains a random image chosen from all patient dicoms after being transformed\n* Tensors of type double are returned (model requires double data type)","metadata":{}},{"cell_type":"code","source":"class OsicDataset(Dataset):\n    def __init__(self, df, mode = 'train'):\n        self.mode = mode\n        self.df = df\n        self.transforms = transforms.Compose([\n            transforms.ToTensor(),\n#             transforms.Normalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225)),\n        ])\n    \n    def __len__(self):\n        return len(self.df)\n    \n    def __getitem__(self, idx):\n        # particular row of dataset\n        row = self.df.iloc[idx,:]\n        # patient id for certain row\n        pid = row['Patient']\n        # fvc measurement for particular patient\n        label = row['FVC']\n        # list of all dicom file paths within patient directory is obtained with ```glob.glob```\n        img_ids = glob.glob('/kaggle/input/osic-pulmonary-fibrosis-progression/' + self.mode + '/' + pid + '/*')\n        # random dicom file path is selected and transformed through the ```get_img``` function\n        img = get_img(np.random.choice(img_ids))\n        # Image is converted to PyTorch tensor\n        img = self.transforms(img)\n        \n        # Tensors are changed to ```torch.double``` datatype\n        # returns dictionary of image and target\n        return {\n            'img': img.type(torch.DoubleTensor),\n            'target': torch.tensor(label, dtype=torch.double),\n        }","metadata":{"execution":{"iopub.status.busy":"2021-08-17T00:28:36.003531Z","iopub.execute_input":"2021-08-17T00:28:36.003859Z","iopub.status.idle":"2021-08-17T00:28:36.010802Z","shell.execute_reply.started":"2021-08-17T00:28:36.003831Z","shell.execute_reply":"2021-08-17T00:28:36.009756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model Specification\nTwo models are created:\n* CustomModel\n    * custom model created using prior code\n* EffModel\n    * model that makes use of pretrained neural networks: created with help from https://www.kaggle.com/noelmat/pytorch-efficientnet-with-better-dataloaders\n    * ```eff_name``` specifies which model (from b0-b7) is to be used as the pretrained model - here, ```b1``` is used","metadata":{}},{"cell_type":"code","source":"class CustomModel(nn.Module):\n    def __init__(self):\n        super().__init__()\n        self.conv1 = nn.Conv2d(1,32,kernel_size=(4,4),padding=1,stride=2)\n        self.max_pool = nn.MaxPool2d(kernel_size=(2, 2))\n        self.conv2 = nn.Conv2d(32, 64, kernel_size=(4, 4), padding=1, stride=2)\n        self.bn1 = nn.BatchNorm2d(64)\n        self.dropout = nn.Dropout(0.15)\n        self.conv3 = nn.Conv2d(64, 4, kernel_size=(4, 4))\n        self.bn2 = nn.BatchNorm2d(4)\n        self.ada_pool = nn.AdaptiveMaxPool2d((28, 28))\n        self.fc1 = nn.Linear(28 * 28 * 4, 16)\n        self.fc2 = nn.Linear(16, 8)\n        self.fc3 = nn.Linear(8, 1)\n        \n        self.relu = nn.ReLU()\n        \n    def forward(self, x):\n        # first convolutional layer with convs of size 4 and 32 channels\n        x = self.conv1(x)\n        # thresholded with relu\n        x = self.relu(x)\n        # max pooling to reduce image size\n        x = self.max_pool(x)\n        # second convolutional layer with convs of size 4 and 64 channels\n        x = self.conv2(x)\n        x = self.relu(x)\n        # batch normalization used to normalize outputs so far\n        x = self.bn1(x)\n        # max pooling to reduce image size\n        x = self.max_pool(x)\n        # dropout to reduce model overfitting\n        x = self.dropout(x)\n        # third convolutional layer with convs of size 4 and 4 channels\n        x = self.conv3(x)\n        x = self.relu(x)\n        x = self.bn2(x)\n        # adaptive pooling to resize images to 228 by 228\n        x = self.ada_pool(x)\n        \n        # flatten image to be passed through fully connected layers\n        x = x.view(x.size(0), -1)\n        x = self.fc1(x)\n        x = self.relu(x)\n        x = self.fc2(x)\n        x = self.relu(x)\n        x = self.fc3(x)\n        \n        return x","metadata":{"execution":{"iopub.status.busy":"2021-08-17T01:50:07.321965Z","iopub.execute_input":"2021-08-17T01:50:07.3223Z","iopub.status.idle":"2021-08-17T01:50:07.333693Z","shell.execute_reply.started":"2021-08-17T01:50:07.32227Z","shell.execute_reply":"2021-08-17T01:50:07.332723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class EffModel(nn.Module):\n    def __init__(self, eff_name='b0'):\n        super().__init__()\n        self.conv = nn.Conv2d(1, 3, kernel_size=3, padding=1, stride=2)\n        self.bn = nn.BatchNorm2d(3)\n        self.model = EfficientNet.from_pretrained(f'efficientnet-{eff_name}')\n        self.fc1 = nn.Linear(1000, 128)\n        self.fc2 = nn.Linear(128, 1)\n        self.relu = nn.ReLU()\n    \n    def forward(self, x):\n        # one initial convolution layer with convs of size 3 and 3 channels\n        x = self.conv(x)\n        # batch normalization\n        x = self.bn(x)\n        x = self.relu(x)\n        \n        # pass through pretrained model\n        x = self.model(x)\n        # image passed through two additional fully connected layers\n        x = self.fc1(x)\n        x = self.relu(x)\n        x = self.fc2(x)\n        \n        return x","metadata":{"execution":{"iopub.status.busy":"2021-08-17T01:52:54.135447Z","iopub.execute_input":"2021-08-17T01:52:54.135815Z","iopub.status.idle":"2021-08-17T01:52:54.143272Z","shell.execute_reply.started":"2021-08-17T01:52:54.135785Z","shell.execute_reply":"2021-08-17T01:52:54.142378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# b1 model is being used\nnet = EffModel(eff_name='b1')\n# model set to type torch.double\nnet.double();","metadata":{"execution":{"iopub.status.busy":"2021-08-17T00:42:13.629025Z","iopub.execute_input":"2021-08-17T00:42:13.629468Z","iopub.status.idle":"2021-08-17T00:42:15.181212Z","shell.execute_reply.started":"2021-08-17T00:42:13.629416Z","shell.execute_reply":"2021-08-17T00:42:15.180236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train Function\n* Created with use of previous mimic_bert code\n* Allows GPU training\n* Mean Squared Error used as loss function\n* Adam used as optimizer function\n* MSE loss, RMSE, and R^2 calculated during training and validation\n* RMSE and R^2 calculated during train validation, and test\n* Loss curve of training and validation loss is plotted","metadata":{}},{"cell_type":"code","source":"def train(model, train_dataloader, valid_dataloader, test_dataloader, epochs):\n    # have torch recognize GPU if it exists, otherwise use CPU\n    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n    # load model onto the device\n    model.to(device)\n    \n    # set loss to MSELoss and optimizer to Adam\n    optimizer = torch.optim.Adam(model.parameters(), lr=0.001)\n    criterion = nn.MSELoss()\n    \n    print(\"start of training loop\")\n    \n    # creating lists to store training and validation losses\n    train_losses = []\n    val_losses = []\n    \n    # iterating model training for each epoch\n    for epoch in track(range(epochs), description='Training...'):\n        model.train()\n        train_loss = 0\n        train_rmse = 0\n        train_rsq = 0\n        count = 0\n        \n        for batch in train_dataloader:\n            # set to zero so gradients are not accumulated\n            optimizer.zero_grad()\n            \n            # load images and targets onto device\n            imgs = batch['img'].to(device)\n            targets = batch['target'].to(device)\n            \n            # pass images through model\n            out = model(imgs)\n            \n            # calculate loss and backpropagate through network\n            loss = criterion(out, targets.view(-1, 1))\n            loss.backward()\n            optimizer.step()\n            \n            # add train_loss for particular training step to epoch training loss\n            train_loss += loss.item() * targets.size(0)\n            \n            # get predicted and actual values for targets\n            predicted = out.cpu().detach().numpy()\n            actual = targets.cpu().detach().numpy()\n            \n            # use predicted and actual values for target to calculate rmse and r^2 score\n            train_rmse += mean_squared_error(actual, predicted, squared=False)\n            train_rsq  += r2_score(actual, predicted)\n            count += 1\n        \n        # divide train_loss by number of training steps in order to get epoch train_loss\n        train_loss /= len(train_dataloader.sampler)\n        train_losses.append(train_loss)\n        \n        # divide rmse and r^2 to get epoch rmse and r^2 score\n        train_rmse /= count\n        train_rsq /= count\n        \n        # Validation\n        model.eval()\n        val_loss = 0\n        val_rmse = 0\n        val_rsq = 0\n        count = 0\n        \n        for batch in val_dataloader:\n            optimizer.zero_grad()\n            \n            imgs = batch['img'].to(device)\n            targets = batch['target'].to(device)\n            \n            out = model(imgs)\n            loss = criterion(out, targets.view(-1, 1))\n            \n            val_loss += loss.item() * targets.size(0)\n            \n            predicted = out.cpu().detach().numpy()\n            actual = targets.cpu().detach().numpy()\n            \n            val_rmse += mean_squared_error(actual, predicted, squared=False)\n            val_rsq  += r2_score(actual, predicted)\n            count += 1\n        \n        val_loss /= len(val_dataloader.sampler)\n        val_losses.append(val_loss)\n        \n        val_rmse /= count\n        val_rsq /= count\n        \n        # print metrics\n        print(\n            \"\\n\",\n            \"\\n\",\n            f\"Epoch {epoch+1}/{epochs}:\\n\",\n            f\"Train loss: {train_loss:.3f}...\\n\",\n            f\"Valid loss: {val_loss:.3f}...\\n\",\n            \"\\n\",\n            f\"Train RMSE: {train_rmse:.3f}...\\n\",\n            f\"Valid RMSE: {val_rmse:.3f}...\\n\",\n            \"\\n\",\n            f\"Train R^2: {train_rsq}...\\n\",\n            f\"Valid R^2: {val_rsq}...\\n\",\n        )\n    \n    \n    # Test\n    model.eval()\n    \n    test_rmse = 0\n    test_rsq = 0\n    count = 0\n\n    for batch in test_dataloader:\n\n        imgs = batch['img'].to(device)\n        targets = batch['target'].to(device)\n\n        out = model(imgs)\n    \n        predicted = out.cpu().detach().numpy()\n        actual = targets.cpu().detach().numpy()\n        \n        test_rmse += mean_squared_error(actual, predicted, squared=False)\n        test_rsq  += r2_score(actual, predicted)\n        count += 1\n\n    test_rmse /= count\n    test_rsq /= count\n\n    print(\n        \"\\n\",\n        \"\\n\",\n        \"\\n\",\n        \"Training Finished!\\n\",\n        f\"Test RMSE: {test_rmse:.3f}...\\n\",\n        f\"Test R^2: {test_rsq:.3f}...\\n\\n\",\n        f\"Predictions: {predicted}\\n\",\n        f\"Actual: {actual} \\n\\n\",\n    )\n    \n    # display loss curves\n    print('Loss Curves: ')\n    plt.xlabel('Epochs')\n    plt.ylabel('Loss')\n    # plot train losses per epoch\n    plt.plot(list(range(epochs)), train_losses, label='train')\n    # plot validation losses per epoch\n    plt.plot(list(range(epochs)), val_losses, label='valid')\n    plt.legend()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-17T00:42:15.88321Z","iopub.execute_input":"2021-08-17T00:42:15.883552Z","iopub.status.idle":"2021-08-17T00:42:15.903085Z","shell.execute_reply.started":"2021-08-17T00:42:15.883514Z","shell.execute_reply":"2021-08-17T00:42:15.90204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocessing\n\n* provided csv files are read into ```train_df``` and ```test_df``` respectively\n* subjects with dicom files that cause errors due to being stored in a different compressed format were not included\n* Provided train dataframe was partitioned into train and validation with in a ratio of 0.7:0.3","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv') \ntest_df  = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv')","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:17:44.767323Z","iopub.execute_input":"2021-08-20T18:17:44.767793Z","iopub.status.idle":"2021-08-20T18:17:44.787767Z","shell.execute_reply.started":"2021-08-20T18:17:44.76775Z","shell.execute_reply":"2021-08-20T18:17:44.786562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = train_df.drop(np.nonzero(np.array(train_df['Patient'] == 'ID00011637202177653955184',dtype=float))[0], axis=0).reset_index(drop=True)\ntrain_df = train_df.drop(np.nonzero(np.array(train_df['Patient'] == 'ID00052637202186188008618',dtype=float))[0], axis=0).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:17:45.524552Z","iopub.execute_input":"2021-08-20T18:17:45.524984Z","iopub.status.idle":"2021-08-20T18:17:45.538245Z","shell.execute_reply.started":"2021-08-20T18:17:45.524949Z","shell.execute_reply":"2021-08-20T18:17:45.536764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df, val_df = train_test_split(train_df, test_size=0.3)","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:17:46.174148Z","iopub.execute_input":"2021-08-20T18:17:46.174526Z","iopub.status.idle":"2021-08-20T18:17:46.182471Z","shell.execute_reply.started":"2021-08-20T18:17:46.174494Z","shell.execute_reply":"2021-08-20T18:17:46.181072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Loading\n* ```train_df```,```val_df```, and ```test_df``` are passed into train, validation, and test instances of the ```OsicDataset``` class\n* PyTorch DataLoaders are created that prepare data in batches of 32 and randomly shuffle data","metadata":{}},{"cell_type":"code","source":"train_dataset = OsicDataset(train_df)\nval_dataset = OsicDataset(val_df)\ntest_dataset = OsicDataset(test_df, mode='test')","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:17:47.375677Z","iopub.execute_input":"2021-08-20T18:17:47.376076Z","iopub.status.idle":"2021-08-20T18:17:47.452651Z","shell.execute_reply.started":"2021-08-20T18:17:47.376043Z","shell.execute_reply":"2021-08-20T18:17:47.451011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dataloader = DataLoader(train_dataset, batch_size=32, shuffle=True)\nval_dataloader = DataLoader(val_dataset, batch_size=32, shuffle=True)\ntest_dataloader = DataLoader(test_dataset, batch_size=32, shuffle=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-20T18:17:47.759892Z","iopub.execute_input":"2021-08-20T18:17:47.760283Z","iopub.status.idle":"2021-08-20T18:17:47.785781Z","shell.execute_reply.started":"2021-08-20T18:17:47.760247Z","shell.execute_reply":"2021-08-20T18:17:47.783944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Loop\n* ```net``` (convolutional neural network chosen) is passed into train function\n* train, validation, and test dataloaders are passed into train function\n* 10 epochs of training are specified","metadata":{}},{"cell_type":"code","source":"train(\n    net,\n    train_dataloader,\n    val_dataloader,\n    test_dataloader,\n    epochs=10\n)","metadata":{"execution":{"iopub.status.busy":"2021-08-17T00:42:31.595025Z","iopub.execute_input":"2021-08-17T00:42:31.595359Z","iopub.status.idle":"2021-08-17T00:46:08.668349Z","shell.execute_reply.started":"2021-08-17T00:42:31.595331Z","shell.execute_reply":"2021-08-17T00:46:08.667375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Cleanup\n* gc package used for garbage collection to free RAM memory\n* ```torch.cuda.empty_cache()``` used to free GPU VRAM memory","metadata":{}},{"cell_type":"code","source":"import gc\ngc.collect()\ntorch.cuda.empty_cache()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Combining both a CNN and RNN to include Tabular Data","metadata":{}},{"cell_type":"markdown","source":"## Preprocessing","metadata":{}},{"cell_type":"code","source":"# read in data\ntrain_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv')\ntest_df = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv')","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:30.989502Z","iopub.execute_input":"2021-08-29T21:28:30.989836Z","iopub.status.idle":"2021-08-29T21:28:31.018423Z","shell.execute_reply.started":"2021-08-29T21:28:30.989807Z","shell.execute_reply":"2021-08-29T21:28:31.017608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# dropping patients who have dicoms that raise errors when trying to open with pydicom\ntrain_df = train_df.drop(np.nonzero(np.array(train_df['Patient'] == 'ID00011637202177653955184',dtype=float))[0], axis=0).reset_index(drop=True)\ntrain_df = train_df.drop(np.nonzero(np.array(train_df['Patient'] == 'ID00052637202186188008618',dtype=float))[0], axis=0).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:32.436322Z","iopub.execute_input":"2021-08-29T21:28:32.436644Z","iopub.status.idle":"2021-08-29T21:28:32.452392Z","shell.execute_reply.started":"2021-08-29T21:28:32.436616Z","shell.execute_reply":"2021-08-29T21:28:32.451467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_df.columns)\nprint(train_df['SmokingStatus'].unique())","metadata":{"execution":{"iopub.status.busy":"2021-08-29T21:28:33.125917Z","iopub.execute_input":"2021-08-29T21:28:33.126238Z","iopub.status.idle":"2021-08-29T21:28:33.133739Z","shell.execute_reply.started":"2021-08-29T21:28:33.126208Z","shell.execute_reply":"2021-08-29T21:28:33.132872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# dropping Percent and rearranging columns\ncols = ['Patient', 'Age', 'Sex', 'SmokingStatus', 'Weeks', 'FVC']\ntrain_df = train_df[cols]\ntest_df = test_df[cols]","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.034149Z","iopub.status.idle":"2021-08-29T17:41:41.034836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.035966Z","iopub.status.idle":"2021-08-29T17:41:41.036529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# label encoding Sex and SmokingStatus\ntrain_df['Sex'] = train_df['Sex'].astype('category').cat.codes\ntrain_df['SmokingStatus'] = train_df['SmokingStatus'].astype('category').cat.codes\n\ntest_df['Sex'] = test_df['Sex'].astype('category').cat.codes\ntest_df['SmokingStatus'] = test_df['SmokingStatus'].astype('category').cat.codes","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.037597Z","iopub.status.idle":"2021-08-29T17:41:41.038315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train_df.info())\nprint(test_df.info())","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.039365Z","iopub.status.idle":"2021-08-29T17:41:41.040085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# taken from https://machinelearningmastery.com/multivariate-time-series-forecasting-lstms-keras/\ndef series_to_supervised(data, n_in=1, n_out=1, dropnan=True):\n    # get number of variables to shift\n    n_vars = 1 if type(data) is list else data.shape[1]\n    df = pd.DataFrame(data)\n    cols, names = list(), list()\n    \n    # shift feature data backward in order to use past n_in time points to predict next fvc values\n    for i in range(n_in, 0, -1):\n        cols.append(df.shift(i))\n        names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]\n    \n    # create target column(s) - n_out time points to be predicted\n    for i in range(0, n_out):\n        cols.append(df.shift(-i))\n        if i == 0:\n            names += [('var%d(t)' % (j+1)) for j in range(n_vars)]\n        else:\n            names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]\n    \n    # concatenating all shift columns together\n    agg = pd.concat(cols, axis=1)\n    agg.columns = names\n    \n    # drop nan values\n    if dropnan:\n        agg.dropna(inplace=True)\n    return agg","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.041153Z","iopub.status.idle":"2021-08-29T17:41:41.041762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get dfs for each patient\ntrain_split = [subj_data for _, subj_data in train_df.groupby('Patient', as_index=False)]\n# get values of each df in train_split\ntrain_vals = [df.values for df in train_split]\n# apply series to supervised function on the values\ntrain_time_series = [series_to_supervised(val) for val in train_vals]\n# sort values in time series by week just in case measurements in the original dataset are not in order\ntrain_time_series = [ts.sort_values(by='var5(t-1)') for ts in train_time_series]","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.043529Z","iopub.execute_input":"2021-08-29T17:41:41.043891Z","iopub.status.idle":"2021-08-29T17:41:41.062826Z","shell.execute_reply.started":"2021-08-29T17:41:41.04385Z","shell.execute_reply":"2021-08-29T17:41:41.061417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# concatenate list of time series for each patient\ntrain_ts = pd.concat(train_time_series)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.064064Z","iopub.status.idle":"2021-08-29T17:41:41.064634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# only predicting next fvc measurement, not other variables\ndrop_idxs = list(range(6, train_ts.shape[1] - 1))\ntrain_ts = train_ts.drop(train_ts.columns[drop_idxs], axis=1)\ntest_ts = test_df","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.065906Z","iopub.status.idle":"2021-08-29T17:41:41.066468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ts","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.103439Z","iopub.execute_input":"2021-08-29T17:41:41.103759Z","iopub.status.idle":"2021-08-29T17:41:41.114657Z","shell.execute_reply.started":"2021-08-29T17:41:41.103726Z","shell.execute_reply":"2021-08-29T17:41:41.11312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# label encoding the patient column\ntrain_ts['Patient'] = train_ts['var1(t-1)'].astype('category').cat.codes","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.116023Z","iopub.status.idle":"2021-08-29T17:41:41.116807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# rearrange columns to make the df look nicer\ntrain_ts = train_ts[['var1(t-1)', 'Patient', 'var2(t-1)', 'var3(t-1)', 'var4(t-1)', 'var5(t-1)', 'var6(t-1)', 'var6(t)']]","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.118132Z","iopub.status.idle":"2021-08-29T17:41:41.1189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ts","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.120227Z","iopub.status.idle":"2021-08-29T17:41:41.120993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_ts","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.122074Z","iopub.status.idle":"2021-08-29T17:41:41.122786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\nsubmission = pd.DataFrame(submission.Patient_Week.str.split('_',1).tolist(),\n                                 columns = ['Patient','Weeks'])","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.124214Z","iopub.execute_input":"2021-08-29T17:41:41.124698Z","iopub.status.idle":"2021-08-29T17:41:41.140359Z","shell.execute_reply.started":"2021-08-29T17:41:41.124661Z","shell.execute_reply":"2021-08-29T17:41:41.139085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.141444Z","iopub.status.idle":"2021-08-29T17:41:41.142146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CombinedDataset Class\n* Relatively similar to the ```OsicDataset``` but includes other feature data to be used in time series prediction (smoking status, gender, etc.)","metadata":{}},{"cell_type":"code","source":"class CombinedDataset(Dataset):\n    def __init__(self, df):\n        self.df = df\n        self.transforms = transforms.Compose([\n            transforms.ToTensor(),\n#             transforms.Normalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225)),\n        ])\n    \n    def __len__(self):\n        return len(self.df)\n    \n    def __getitem__(self, idx):\n        # particular row of dataset\n        row = self.df.iloc[idx,:]\n\n        # patient id for certain row\n        pid = row['var1(t-1)']\n        # fvc measurement for particular patient\n        label = row['var6(t)']\n        \n        ## Image Processing\n        \n        # list of all dicom file paths within patient directory is obtained with ```glob.glob```\n        img_ids = glob.glob('/kaggle/input/osic-pulmonary-fibrosis-progression/train/' + pid + '/*')\n        # random dicom file path is selected and transformed through the ```get_img``` function\n        img = get_img(np.random.choice(img_ids))\n        # Image is converted to PyTorch tensor\n        img = self.transforms(img)\n        \n        ## Time Series Processing\n        # don't include pid char string or label at the end\n        ts_data = np.array(list(row[1:-1].values))\n        ts_data = ts_data.reshape(1, 6)\n        # Tensors are changed to ```torch.double``` datatype\n        # returns dictionary of image and target\n        return {\n            'img': img.type(torch.DoubleTensor),\n            'ts': torch.tensor(ts_data, dtype=torch.double),\n            'target': torch.tensor(label, dtype=torch.double),\n        }","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.143216Z","iopub.status.idle":"2021-08-29T17:41:41.143965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CombinedModel Class\n* similar to ```EffModel``` but includes the time series data through the use of an LSTM\n* hidden_size of the LSTM can be changed\n* output tensors from RNN and CNN are concatenated and passed through a last fully connected layer to create output","metadata":{}},{"cell_type":"code","source":"class CombinedModel(nn.Module):\n    def __init__(self, hidden_size, eff_name='b1'):\n        super().__init__()\n        \n        ## img layers\n        self.conv = nn.Conv2d(1, 3, kernel_size=3, padding=1, stride=2)\n        self.bn = nn.BatchNorm2d(3)\n        self.model = EfficientNet.from_pretrained(f'efficientnet-{eff_name}')\n        self.fc1 = nn.Linear(1000, 500)\n        \n        ## ts layer\n        self.hidden_size = hidden_size\n        self.rnn = nn.LSTM(input_size=6,hidden_size=self.hidden_size, num_layers=2, dropout=0.4, batch_first=True)\n        \n        self.output = nn.Linear(500 + hidden_size, 1)\n        self.relu = nn.ReLU()\n        \n    def forward(self, img, ts):\n        # apply conv to increase channels to 3\n        img = self.conv(img)\n        # batch norm to \n        img = self.bn(img)\n        img = self.relu(img)\n        # pass through efficientnet model\n        img = self.model(img)\n        # fc layer to condense outputs down to 500\n        img = self.fc1(img)\n        img = self.relu(img)\n        \n        # run rnn on time series data\n        ts, _ = self.rnn(ts)\n        \n        # concatenate image features with time series features\n        x = torch.cat([img, ts.view(-1,self.hidden_size)], dim=1)\n        # pass through last fc layer\n        return self.output(x)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.145074Z","iopub.status.idle":"2021-08-29T17:41:41.145776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_loop(model, train_dataloader, val_dataloader, epochs):\n    \n    # have torch recognize GPU if it exists, otherwise use CPU\n    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n    # load model onto the device\n    model.to(device)\n    \n    # set loss to MSELoss and optimizer to Adam\n    optimizer = torch.optim.Adam(model.parameters(), lr=0.001)\n    criterion = nn.L1Loss()\n    \n    print(\"start of training loop\")\n    \n    # creating lists to store training and validation losses\n    train_losses = []\n    val_losses = []\n    \n    # iterating model training for each epoch\n    for epoch in track(range(epochs), description='Training...'):\n        model.train()\n        train_loss = 0\n        train_rmse = 0\n        train_rsq = 0\n        count = 0\n        \n        for batch in train_dataloader:\n            # set to zero so gradients are not accumulated\n            optimizer.zero_grad()\n            \n            # load images, time series data, and targets onto device\n            imgs = batch['img'].to(device)\n            ts = batch['ts'].to(device)\n            targets = batch['target'].to(device)\n            \n            # pass images and time series data through model\n            out = model(imgs, ts)\n            \n            # calculate loss and backpropagate through network\n            loss = criterion(out, targets.view(-1, 1))\n            loss.backward()\n            optimizer.step()\n            \n            # add train_loss for particular training step to epoch training loss\n            train_loss += loss.item() * targets.size(0)\n            \n            # get predicted and actual values for targets\n            predicted = out.cpu().detach().numpy()\n            actual = targets.cpu().detach().numpy()\n            \n            # use predicted and actual values for target to calculate rmse and r^2 score\n            train_rmse += mean_squared_error(actual, predicted, squared=False)\n            train_rsq  += r2_score(actual, predicted)\n            count += 1\n        \n        # divide train_loss by number of training steps in order to get epoch train_loss\n        train_loss /= len(train_dataloader.sampler)\n        train_losses.append(train_loss)\n        \n        # divide rmse and r^2 to get epoch rmse and r^2 score\n        train_rmse /= count\n        train_rsq /= count\n        \n        # Validation\n        model.eval()\n        val_loss = 0\n        val_rmse = 0\n        val_rsq = 0\n        count = 0\n        \n        for batch in val_dataloader:\n            optimizer.zero_grad()\n            \n            imgs = batch['img'].to(device)\n            ts = batch['ts'].to(device)\n            targets = batch['target'].to(device)\n            \n            out = model(imgs, ts)\n            loss = criterion(out, targets.view(-1, 1))\n            \n            val_loss += loss.item() * targets.size(0)\n            \n            predicted = out.cpu().detach().numpy()\n            actual = targets.cpu().detach().numpy()\n            \n            val_rmse += mean_squared_error(actual, predicted, squared=False)\n            val_rsq  += r2_score(actual, predicted)\n            count += 1\n        \n        val_loss /= len(val_dataloader.sampler)\n        val_losses.append(val_loss)\n        \n        val_rmse /= count\n        val_rsq /= count\n        \n        # print metrics\n        print(\n            \"\\n\",\n            \"\\n\",\n            f\"Epoch {epoch+1}/{epochs}:\\n\",\n            f\"Train loss: {train_loss:.3f}...\\n\",\n            f\"Valid loss: {val_loss:.3f}...\\n\",\n            \"\\n\",\n            f\"Train RMSE: {train_rmse:.3f}...\\n\",\n            f\"Valid RMSE: {val_rmse:.3f}...\\n\",\n            \"\\n\",\n            f\"Train R^2: {train_rsq}...\\n\",\n            f\"Valid R^2: {val_rsq}...\\n\",\n        )\n    \n    \n    # Test\n    model.eval()\n    \n    test_rmse = 0\n    test_rsq = 0\n    count = 0\n\n#     for batch in test_dataloader:\n\n#         imgs = batch['img'].to(device)\n#         targets = batch['target'].to(device)\n\n#         out = model(imgs, targets)\n    \n#         predicted = out.cpu().detach().numpy()\n#         actual = targets.cpu().detach().numpy()\n        \n#         test_rmse += mean_squared_error(actual, predicted, squared=False)\n#         test_rsq  += r2_score(actual, predicted)\n#         count += 1\n\n#     test_rmse /= count\n#     test_rsq /= count\n\n    print(\n        \"\\n\",\n        \"\\n\",\n        \"\\n\",\n        \"Training Finished!\\n\",\n#         f\"Predictions vs Actual: {submission}\",\n    )\n    \n    # display loss curves\n    print('Loss Curves: ')\n    plt.xlabel('Epochs')\n    plt.ylabel('Loss')\n    # plot train losses per epoch\n    plt.plot(list(range(epochs)), train_losses, label='train')\n    # plot validation losses per epoch\n    plt.plot(list(range(epochs)), val_losses, label='valid')\n    plt.legend()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.150941Z","iopub.execute_input":"2021-08-29T17:41:41.15149Z","iopub.status.idle":"2021-08-29T17:41:41.243234Z","shell.execute_reply.started":"2021-08-29T17:41:41.151455Z","shell.execute_reply":"2021-08-29T17:41:41.242341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# randomly sample certain patients to be in the train set and others to be in validation\nrng = np.random.default_rng()\nnum_patients = train_ts['Patient'].nunique()\nval_idxs = rng.choice(num_patients, size=math.floor(num_patients * 0.3), replace=False)\ntrain = train_ts.loc[~train_ts['Patient'].isin(val_idxs)].copy()\nvalid = train_ts.loc[train_ts['Patient'].isin(val_idxs)].copy()","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.2449Z","iopub.execute_input":"2021-08-29T17:41:41.245442Z","iopub.status.idle":"2021-08-29T17:41:41.266553Z","shell.execute_reply.started":"2021-08-29T17:41:41.245406Z","shell.execute_reply":"2021-08-29T17:41:41.264884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train.shape)\nprint(valid.shape)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.267732Z","iopub.status.idle":"2021-08-29T17:41:41.268567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# instantiate datasets and dataloaders with a batch size of 64\ntrain_dataset = CombinedDataset(train)\nvalid_dataset = CombinedDataset(valid)\n\ntrain_dataloader = DataLoader(train_dataset, batch_size=64)\nvalid_dataloader = DataLoader(train_dataset, batch_size=64)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.269769Z","iopub.status.idle":"2021-08-29T17:41:41.270354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create new efficientnet b3 model and convert model to work with double tensors\nnet = CombinedModel(32, eff_name='b3')\nnet.double();","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.271656Z","iopub.status.idle":"2021-08-29T17:41:41.272202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# run training loop\ntrain_loop(net, train_dataloader, valid_dataloader, epochs=30)","metadata":{"execution":{"iopub.status.busy":"2021-08-29T17:41:41.273552Z","iopub.status.idle":"2021-08-29T17:41:41.274094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cleanup","metadata":{}},{"cell_type":"code","source":"import gc\ngc.collect()\ntorch.cuda.empty_cache()","metadata":{},"execution_count":null,"outputs":[]}]}