{"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":"## Imports & Setup","metadata":{}},{"cell_type":"code","source":"# installing necessary libraries\n!pip install rich efficientnet_pytorch\n!pip install --upgrade wandb","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-12-13T00:53:49.984512Z","iopub.execute_input":"2021-12-13T00:53:49.984983Z","iopub.status.idle":"2021-12-13T00:54:12.409145Z","shell.execute_reply.started":"2021-12-13T00:53:49.984881Z","shell.execute_reply":"2021-12-13T00:54:12.408212Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport glob\nimport random\nfrom pydicom import dcmread\nimport matplotlib.pyplot as plt\nfrom rich.progress import track\nfrom sklearn.metrics import mean_squared_error, r2_score\n\nimport torch\nfrom torch import nn\nfrom torch.utils.data import Dataset, DataLoader","metadata":{"execution":{"iopub.status.busy":"2021-12-13T00:54:35.865181Z","iopub.execute_input":"2021-12-13T00:54:35.865495Z","iopub.status.idle":"2021-12-13T00:54:38.193424Z","shell.execute_reply.started":"2021-12-13T00:54:35.865460Z","shell.execute_reply":"2021-12-13T00:54:38.192514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import wandb\n# wandb.login()","metadata":{"execution":{"iopub.status.busy":"2021-12-13T00:54:39.920880Z","iopub.execute_input":"2021-12-13T00:54:39.921175Z","iopub.status.idle":"2021-12-13T00:54:51.446778Z","shell.execute_reply.started":"2021-12-13T00:54:39.921141Z","shell.execute_reply":"2021-12-13T00:54:51.445825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# set seeds\nrandom.seed(42)\nnp.random.seed(42)\ntorch.manual_seed(42)","metadata":{"execution":{"iopub.status.busy":"2021-12-13T00:55:58.284531Z","iopub.execute_input":"2021-12-13T00:55:58.284894Z","iopub.status.idle":"2021-12-13T00:55:58.294235Z","shell.execute_reply.started":"2021-12-13T00:55:58.284862Z","shell.execute_reply":"2021-12-13T00:55:58.293635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Preparation","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')\n# find num of patients in train_df\ntrain_df['Patient'].nunique()","metadata":{"execution":{"iopub.status.busy":"2021-12-13T00:55:58.873650Z","iopub.execute_input":"2021-12-13T00:55:58.874145Z","iopub.status.idle":"2021-12-13T00:55:58.926339Z","shell.execute_reply.started":"2021-12-13T00:55:58.874113Z","shell.execute_reply":"2021-12-13T00:55:58.925663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# encoding patient, sex, and smoking status as nums\ntrain_df['Patient'] = train_df['Patient'].astype('category').cat.codes\ntrain_df['Sex'] = train_df['Sex'].astype('category').cat.codes\ntrain_df['SmokingStatus'] = train_df['SmokingStatus'].astype('category').cat.codes\ntrain_df","metadata":{"execution":{"iopub.status.busy":"2021-12-13T00:55:59.159187Z","iopub.execute_input":"2021-12-13T00:55:59.159482Z","iopub.status.idle":"2021-12-13T00:55:59.195047Z","shell.execute_reply.started":"2021-12-13T00:55:59.159447Z","shell.execute_reply":"2021-12-13T00:55:59.194308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get shifted future values for Weeks and FVC -> goal is to predict FVC at next time point: FVC(t+1)\ndef shift(x: pd.DataFrame):\n    x['Weeks(t+1)'] = x['Weeks'].shift(-1)\n    x['FVC(t+1)'] = x['FVC'].shift(-1)\n    x = x.iloc[:-1]\n    return x\n\n# get dataframe for each patient and map shift function to each dataframe\nsplit = [y for _, y in train_df.groupby('Patient', as_index=False)]\nsplit = list(map(shift, split))\n\n\n# split data into train, validation, and test\nsplit_len = len(split)\nrand_idx = list(range(split_len))\nrandom.shuffle(rand_idx)\n\n# train: 70%, validation: 15%, test 15%\ntrain_list = [split[x] for x in rand_idx[:int(split_len * 0.7)]]\nval_list = [split[x] for x in rand_idx[int(split_len * 0.7): int(split_len * 0.85)]]\ntest_list = [split[x] for x in rand_idx[int(split_len * 0.85):]]\n\ntrain_df = pd.concat(train_list)\nval_df = pd.concat(val_list)\ntest_df = pd.concat(test_list)\ntrain_df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save dataframes\ntrain_df.to_csv('/kaggle/working/train.csv')\nval_df.to_csv('/kaggle/working/val.csv')\ntest_df.to_csv('/kaggle/working/test.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# normalize FVC values and store FVC mean and standard deviation\ntarget_mean = train_df['FVC(t+1)'].mean()\ntarget_stdev = train_df['FVC(t+1)'].std()\nmean_list = []\nstdev_list = []\nfor c in train_df.columns:\n    if c in ['FVC', 'FVC(t+1)']:\n        mean = train_df[c].mean()\n        stdev = train_df[c].std()\n        mean_list.append(mean)\n        stdev_list.append(stdev)\n        train_df[c] = (train_df[c] - mean) / stdev\n        val_df[c] = (val_df[c] - mean) / stdev\n        test_df[c] = (test_df[c] - mean) / stdev\n\n        \n# remove any potential nans before start of training\ntrain_df.dropna(inplace=True)\ntrain_df.reset_index(drop=True, inplace=True)\nval_df.dropna(inplace=True)\nval_df.reset_index(drop=True, inplace=True)\ntest_df.dropna(inplace=True)\ntest_df.reset_index(drop=True, inplace=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## PyTorch Dataset and Model Classes","metadata":{}},{"cell_type":"code","source":"class OASISDataset(Dataset):\n    def __init__(self, dataframe, sequence_length=10):\n        self.sequence_length = sequence_length\n        \n        # get column names for target and features\n        self.target = dataframe.columns[-1]\n        self.features = dataframe.columns[:-1]\n        \n        # save df values as torch tensors\n        self.y = torch.tensor(dataframe[self.target].values).float()\n        self.X = torch.tensor(dataframe[self.features].values).float()\n\n    def __len__(self):\n        return self.X.shape[0]\n\n    def __getitem__(self, i):\n        if i >= self.sequence_length - 1:\n            # get a slice of dataframe for `sequence_length` values\n            i_start = i - self.sequence_length + 1\n            x = self.X[i_start:(i+1), :]\n        else:\n            # add padding of first value if not enough examples before it\n            padding = self.X[0].repeat(self.sequence_length - i - 1, 1)\n            x = self.X[0:(i+1), :]\n            x = torch.cat((padding, x), 0)\n            \n        # return feature and target as dict\n        return {\n            'X': x.type(torch.float),\n            'y': self.y[i].type(torch.float)\n        }","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class OASISModel(nn.Module):\n    def __init__(self, input_size, hidden_units, num_layers):\n        super(OASISModel, self).__init__()\n        self.input_size = input_size\n        self.hidden_units = hidden_units\n        self.num_layers = num_layers\n\n        self.lstm = nn.LSTM(\n            input_size=self.input_size,\n            hidden_size=self.hidden_units,\n            batch_first=True,\n            num_layers=self.num_layers\n        )\n\n        self.linear = nn.Linear(hidden_units, 1)\n\n    def forward(self, x):\n        # simple lstm with zero initialization, output to single linear layer\n        batch_size = x.shape[0]\n        h0 = torch.zeros(self.num_layers, batch_size, self.hidden_units).requires_grad_()\n        c0 = torch.zeros(self.num_layers, batch_size, self.hidden_units).requires_grad_()\n\n        _, (hn, _) = self.lstm(x, (h0, c0))\n        out = self.linear(hn[0])\n        return out","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# model hyperparameters\nbatch_size = 16\nsequence_length = 10\nlearning_rate = 5e-4\nnum_hidden_units = 10\n\n# create dataset class with specified sequence length\ntrain_dataset = OASISDataset(train_df, sequence_length=sequence_length)\nvalid_dataset = OASISDataset(val_df, sequence_length=sequence_length)\ntest_dataset = OASISDataset(test_df, sequence_length=sequence_length)\n\n# create dataloader that does not shuffle data\ntrain_dataloader = DataLoader(train_dataset, batch_size=batch_size, shuffle=False)\nvalid_dataloader = DataLoader(valid_dataset, batch_size=batch_size, shuffle=False)\ntest_dataloader = DataLoader(test_dataset, batch_size=len(test_dataset), shuffle=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# instantiate model\nmodel = OASISModel(input_size=train_dataset.X.shape[1], hidden_units=num_hidden_units, num_layers=3)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_reg(model, train_oader, valid_loader, test_loader, 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=learning_rate)\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#    wandb.init(project='osic_lstm', entity='ashcher51')\n#    wandb.watch(model)\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            # get batch and load on to device\n            X, y = batch['X'].to(device), batch['y'].to(device)\n            \n            # get model output and calculate loss\n            out = model(X)\n            loss = criterion(out, y.view(-1, 1))\n\n            # update params\n            loss.backward()\n            optimizer.step()\n\n            # add train_loss at particular training step to epoch training loss\n            train_loss += loss.item() * y.size(0)\n\n            # get predicted and actual values for targets\n            predicted = out.cpu().detach().numpy()\n            actual = y.cpu().detach().numpy()\n\n            # use predicted and actual 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\n            count += 1\n\n        # divide train_loss by number of training steps 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#        wandb.log({'train_loss':train_loss,'train_rmse':train_rmse,'train_rsq':train_rsq,})\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 valid_dataloader:\n            optimizer.zero_grad()\n\n            X, y = batch['X'].to(device), batch['y'].to(device)\n\n            out = model(X)\n            loss = criterion(out, y.view(-1, 1))\n\n            val_loss += loss.item() * y.size(0)\n\n            predicted = out.cpu().detach().numpy()\n            actual = y.cpu().detach().numpy()\n\n            val_rmse += mean_squared_error(actual, predicted, squared=False)\n            val_rsq += r2_score(actual, predicted)\n        \n\n            count += 1\n\n        val_loss /= len(valid_dataloader.sampler)\n        val_losses.append(val_loss)\n\n        val_rmse /= count\n        val_rsq /= count\n            \n#        wandb.log({'val_loss':val_loss,'val_rmse':val_rmse,'val_rsq':val_rsq,})\n\n\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    # 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        X, y = batch['X'].to(device), batch['y'].to(device)\n\n        out = model(X)\n        loss = criterion(out, y.view(-1, 1))\n\n        test_predicted = out.cpu().detach().numpy()\n        test_actual = y.cpu().detach().numpy()\n\n        test_rmse += mean_squared_error(test_actual, test_predicted, squared=False)\n        test_rsq += r2_score(test_actual, test_predicted)\n\n\n        count += 1\n\n    test_rmse /= count\n    test_rsq /= count\n\n#    wandb.log({'test_rmse':test_rmse,'test_rsq':test_rsq,})\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: {test_predicted}\\n\",\n        f\"Actual: {test_actual} \\n\\n\",\n    )\n    \n    print(f'\\nlen of predicted: {len(test_predicted)}')\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()\n\n    # display scatterplot\n    plt.scatter(test_actual, test_predicted)\n    # plt.scatter(test_actual * target_stdev + target_mean, test_predicted * target_stdev + target_mean) # <- scaled version of scatterplot\n    plt.xlabel('Actual')\n    plt.ylabel('Predicted')\n    plt.show()\n    \n    # get sample patient progression\n    sub_pred = [np.array(test_df.iloc[9, 2] * stdev_list[0] + mean_list[0])] + list(predicted[9:16] * stdev_list[1] + mean_list[1])\n    sub_actual = [np.array(test_df.iloc[9, 2] * stdev_list[0] + mean_list[0])] + list(actual[9:16] * stdev_list[1] + mean_list[1])\n    sub_time = [np.array(test_df.iloc[9, 1])] + list(test_df.iloc[9:16, -2])\n    \n    print(f'actual: {sub_actual}')\n    print(f'pred: {sub_pred}')\n    \n    print(np.array(sub_time).reshape(1, 8))\n    print(np.array(sub_pred).reshape(1, 8))\n    print(np.array(sub_actual).reshape(1, 8))\n    \n    plt.plot(np.array(sub_time).reshape(8, 1), np.array(sub_pred).reshape(8, 1), label='predicted')\n    plt.plot(np.array(sub_time).reshape(8, 1), np.array(sub_actual).reshape(8, 1), label='actual')\n    plt.legend()\n    plt.show()\n\n    # get results and save test_df, scaled FVC, actual values, and predicted values to file\n    res_df = pd.concat([test_df, pd.Series(test_df['FVC'] * stdev_list[0] + mean_list[0], name = 'FVC scaled'), pd.Series(np.array(test_actual) * stdev_list[1] + mean_list[1], name='actual'), pd.Series(np.array([predict[0] for predict in test_predicted]) * stdev_list[1] + mean_list[1], name='predicted')], axis=1)\n    res_df.to_csv('/kaggle/working/results.csv', index=False)\n    \n#    wandb.finish()\n    # simulate progression\n#     sub = 0\n#     weeks = [-1.442765, -1.123454, -0.963798, -0.857361, -0.804142, -0.484831, 0.153792, 0.792414, 2.687218]\n#     FVC = 2.865436\n#     percent = 2.188235\n#     age = 0\n#     sex = 0.513665\n#     smoking = -0.397809\n    \n#     predicted_fvc = []\n    \n#     for i, week in enumerate(weeks[:-1]):\n#         model_in = torch.tensor([sub, week, FVC, percent, age, sex, smoking, weeks[i+1]], dtype=torch.float)\n#         model_out = model(model_in)\n#         predicted_fvc.append(model_out * target_stdev + target_mean)\n#         FVC = model_out\n    \n#     plt.plot(weeks, predicted_fvc)\n#     plt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_reg(model, train_dataloader, valid_dataloader, test_dataloader, 500)","metadata":{},"execution_count":null,"outputs":[]}]}