{"cells":[{"metadata":{},"cell_type":"markdown","source":"# IMPORT"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import copy\nimport random\nimport os\n\nfrom datetime import datetime, timedelta\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nfrom tqdm import tqdm\nfrom tqdm.notebook import trange\nfrom time import time\n\nimport torch\nfrom torch.utils.data import DataLoader, Subset\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset\nfrom torch.autograd import Variable\n\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.model_selection import GroupKFold\n\nLR_SCHEDULE = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# GLOBAL VARIABLES"},{"metadata":{"trusted":true},"cell_type":"code","source":"class CFG:\n    root_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression')\n    model_dir = Path('/kaggle/working')\n    num_kfolds = 5\n    cpu_workers = 4\n    batch_size = 1\n    learning_rate = 1e-2\n    num_epochs = 100\n    quantiles = [0.2, 0.5, 0.8]\n    # LSTM parameters\n    seq_length = 146\n    n_features = 8\n    n_layers = 2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# UTILS"},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=42):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = True\n    torch.backends.cudnn.enabled = True\n    \n    \n# Helper generator that yields kfold PyTorch datasets\ndef group_kfold(dataset, groups, n_splits):\n    gkf = GroupKFold(n_splits=n_splits)\n    for train_idx, val_idx in gkf.split(dataset, dataset, groups):\n        train = Subset(dataset, train_idx)\n        val = Subset(dataset, val_idx)\n        yield train, val\n        \n            \n# Helper function with competition metric\ndef metric(p0, p1, p2, targets):\n    sigma = p2 - p0\n    sigma[sigma < 70] = 70\n    delta = np.absolute(p1 - targets)\n    delta[delta > 1000] = 1000\n    return np.mean(-np.sqrt(2) * delta / sigma - np.log(np.sqrt(2) * sigma), 1)\n\n# Loss\ndef pinball_loss(preds, targets, q):\n    assert not targets.requires_grad\n    assert preds.size(0) == targets.size(0)\n    e = targets - preds\n    loss = torch.max((q - 1) * e, q * e)\n    loss = torch.mean(torch.sum(loss, dim=1))\n    return loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Configure tabular data"},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler_fvc = MinMaxScaler()\nscaler_percent = MinMaxScaler()\nscaler_age = MinMaxScaler()\n\n# Read train csv\ntrain_df = pd.read_csv(Path(CFG.root_dir)/\"train.csv\")\ntrain_df.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\n\n# Normalize features\ntrain_df[['FVC']] = scaler_fvc.fit_transform(train_df[['FVC']])\ntrain_df[['Percent']] = scaler_percent.fit_transform(train_df[['Percent']])\ntrain_df[['Age']] = scaler_age.fit_transform(train_df[['Age']])\n\n# Create sequences \ndf = pd.merge(train_df.groupby('Patient')['Weeks'].apply(list).to_frame().reset_index(), \\\n              train_df.groupby('Patient')['FVC'].apply(list).to_frame().reset_index(), on=\"Patient\")\ndf = pd.merge(df,  train_df.groupby('Patient')['Percent'].apply(list).to_frame().reset_index(), on=\"Patient\")\ndf = pd.merge(df, train_df.groupby('Patient')['Age'].first(), on=\"Patient\")\ndf = pd.merge(df, train_df.groupby('Patient')['Sex'].first(), on=\"Patient\")\ndf = pd.merge(df, train_df.groupby('Patient')['SmokingStatus'].first(), on=\"Patient\")\n\n# Convert sex and smoking status to one-hot encoding\nCOLS = ['Sex', 'SmokingStatus']\nfor col in COLS:\n    for mod in df[col].unique():\n        df[mod] = (df[col] == mod).astype(int)\n\n# Remove useless columns\ndf.drop(columns=['Sex', 'SmokingStatus'], inplace=True)\n        \ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"t = df['FVC'].iloc[90]\nplt.plot(t)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# DEFINE DATASET"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"class ClinicalDataset(Dataset):\n    def __init__(self, df):\n        self.df = df\n        self.N = CFG.seq_length\n        self.feat = CFG.n_features\n\n    def __len__(self):\n        return len(self.df.Patient.values)\n\n    def __getitem__(self, idx):\n        if torch.is_tensor(idx):\n            idx = idx.tolist()\n\n        # Extract row from dataframe\n        row = self.df.iloc[idx].values\n        \n        # Get features \n        patiend_id = row[0]\n        weeks = np.array(row[1]) + 12\n        fvc = np.array(row[2])\n        percent = np.array(row[3])\n        age =  np.array(row[4])\n        male = np.array(row[5])\n        female = np.array(row[6])\n        ex_smok = np.array(row[7])\n        n_smok = np.array(row[8])\n        smok = np.array(row[9])\n        \n        # Create input FVC sequence\n        seq = np.zeros((self.N, self.feat))\n        seq[0, 0] = fvc[0]\n        seq[0, 1] = percent[0]\n        seq[:, 2].fill(age)\n        seq[:, 3].fill(male)\n        seq[:, 4].fill(female)\n        seq[:, 5].fill(ex_smok)\n        seq[:, 6].fill(n_smok)\n        seq[:, 7].fill(smok)\n        \n        # Create expected FVCs (starting just after the initial FVC)\n        gt = np.zeros(self.N)\n        out_ids = weeks[1:] - weeks[0] - 1\n        gt[out_ids] = fvc[1:]\n\n        return torch.tensor(seq), torch.tensor(gt), out_ids.tolist()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Check dataset**"},{"metadata":{"trusted":true},"cell_type":"code","source":"df.iloc[0].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset = ClinicalDataset(df)\n\nprint('Input sequence \\n', dataset[0][0])\nprint('\\nOutput sequence \\n', dataset[0][1])\nprint('\\nOutput ids \\n', dataset[0][2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# DEFINE NETWORK AND LOSS"},{"metadata":{"trusted":true},"cell_type":"code","source":"class LSTM(nn.Module):\n\n    def __init__(self, device, input_size, hidden_size, num_layers):\n        super(LSTM, self).__init__()\n        \n        self.device = device\n        self.num_layers = num_layers\n        self.input_size = input_size\n        self.hidden_size = hidden_size\n        \n        self.lstm = nn.LSTM(input_size=input_size, hidden_size=hidden_size,\n                            num_layers=num_layers, batch_first=True, dropout=0.5)\n\n    def forward(self, x):\n        h_0 = Variable(torch.zeros(\n            self.num_layers, x.size(0), self.hidden_size).to(self.device))\n        \n        c_0 = Variable(torch.zeros(\n            self.num_layers, x.size(0), self.hidden_size).to(self.device))\n        \n        # Propagate input through LSTM\n        ula, (hn, _) = self.lstm(x, (h_0, c_0))\n        \n        y = hn.view(-1, self.hidden_size)\n        h_out = hn.view(self.num_layers, x.size(0), self.hidden_size)[-1]\n        \n        h_out = h_out.view(-1, self.hidden_size)\n       \n        return h_out","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# TRAIN"},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"# Get current device\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\nseed_everything()\n\nmodels = []\nscores = []\n\n# KFold \ngkf = GroupKFold(n_splits=CFG.num_kfolds)\n\nfor i, (train_idx, val_idx) in enumerate(gkf.split(df.values, df.values, df['Patient'].values)): \n    df.loc[val_idx, 'fold'] = i\n\ndf['fold'] = df['fold'].astype(int)\n\nt0 = time()\n\n# Loop through folds\nfor fold in range(CFG.num_kfolds):\n    \n    # Get fold ids\n    trn_idx = df[df['fold'] != fold].index\n    val_idx = df[df['fold'] == fold].index\n\n    df_train = df.iloc[trn_idx].reset_index(drop=True)\n    df_valid = df.iloc[val_idx].reset_index(drop=True)\n    \n    # Create dataset\n    train_dataset = ClinicalDataset(df_train)\n    valid_dataset = ClinicalDataset(df_valid)\n    \n    # Create dataloaders\n    dataset_sizes = {'train': len(train_dataset), \n                     'val': len(valid_dataset)}\n    \n    dataloaders = {\n        'train': DataLoader(train_dataset, batch_size=CFG.batch_size,\n                            shuffle=True, num_workers=CFG.cpu_workers),\n        'val': DataLoader(valid_dataset, batch_size=CFG.batch_size,\n                          shuffle=False, num_workers=CFG.cpu_workers)\n    }\n    # Create the models\n    modelq0 = LSTM(device=device, input_size=CFG.n_features, hidden_size=CFG.seq_length, num_layers=CFG.n_layers)\n    modelq1 = LSTM(device=device, input_size=CFG.n_features, hidden_size=CFG.seq_length, num_layers=CFG.n_layers)\n    modelq2 = LSTM(device=device, input_size=CFG.n_features, hidden_size=CFG.seq_length, num_layers=CFG.n_layers)\n    \n    modelq0.to(device)\n    modelq1.to(device)\n    modelq2.to(device)\n    \n    # Optimizer\n    optimizerq0 = torch.optim.Adam(modelq0.parameters(), lr=CFG.learning_rate)\n    optimizerq1 = torch.optim.Adam(modelq1.parameters(), lr=CFG.learning_rate)\n    optimizerq2 = torch.optim.Adam(modelq2.parameters(), lr=CFG.learning_rate)\n    \n    # Schedulers\n    if LR_SCHEDULE == True:\n        schedulerq0 = torch.optim.lr_scheduler.CosineAnnealingLR(optimizerq0, CFG.num_epochs)\n        schedulerq1 = torch.optim.lr_scheduler.CosineAnnealingLR(optimizerq1, CFG.num_epochs)\n        schedulerq2 = torch.optim.lr_scheduler.CosineAnnealingLR(optimizerq2, CFG.num_epochs)\n    \n    epoch_loss_train_0 = np.zeros(CFG.num_epochs)\n    epoch_loss_train_1 = np.zeros(CFG.num_epochs)\n    epoch_loss_train_2 = np.zeros(CFG.num_epochs)\n    \n    epoch_loss_val_0 = np.zeros(CFG.num_epochs)\n    epoch_loss_val_1 = np.zeros(CFG.num_epochs)\n    epoch_loss_val_2 = np.zeros(CFG.num_epochs)\n    m = np.zeros(CFG.num_epochs)\n    \n    # Loop through epochs\n    bar = trange(CFG.num_epochs, desc=f'Training fold {fold + 1}')\n    for epoch in bar:\n        # Train\n        modelq0.train()\n        modelq1.train()\n        modelq2.train()\n        \n        for in_seq, out_seq, mask in dataloaders['train']:\n            \n            # Get input and target sequences\n            inputs = in_seq.float().to(device) # [bs, N, n_feat]\n            targets = out_seq.to(device) # [bs, N]\n\n            # First quantile model\n            optimizerq0.zero_grad()\n            preds = modelq0(inputs) # [bs, N]\n            loss0 = pinball_loss(preds[:, mask], targets[:, mask], CFG.quantiles[0])\n            loss0.backward()\n            optimizerq0.step()\n            epoch_loss_train_0[epoch] += loss0.item()\n            \n            # Second quantile model\n            optimizerq1.zero_grad()\n            preds = modelq1(inputs) # [bs, N]\n            loss1 = pinball_loss(preds[:, mask], targets[:, mask], CFG.quantiles[1])\n            loss1.backward()\n            optimizerq1.step()\n            epoch_loss_train_1[epoch] += loss1.item()\n            \n            # Third quantile model\n            optimizerq2.zero_grad()\n            preds = modelq2(inputs) # [bs, N]\n            loss2 = pinball_loss(preds[:, mask], targets[:, mask], CFG.quantiles[2])\n            loss2.backward()\n            optimizerq2.step()\n            epoch_loss_train_2[epoch] += loss2.item()\n\n        # Epoch losses\n        epoch_loss_train_0[epoch] = epoch_loss_train_0[epoch] / dataset_sizes['train']\n        epoch_loss_train_1[epoch] = epoch_loss_train_1[epoch] / dataset_sizes['train']\n        epoch_loss_train_2[epoch] = epoch_loss_train_2[epoch] / dataset_sizes['train']\n        \n        # Validate\n        modelq0.eval()\n        modelq1.eval()\n        modelq2.eval()\n        \n        for in_seq, out_seq, mask in dataloaders['val']:\n            \n            # Get input and target sequences\n            inputs = in_seq.float().to(device) # [bs, N, n_feat]\n            targets = out_seq.to(device) # [bs, N]\n            \n            # Inference\n            preds0 = modelq0(inputs) # [bs, N]\n            preds1 = modelq1(inputs) # [bs, N]\n            preds2 = modelq2(inputs) # [bs, N]\n            \n            # Losses\n            epoch_loss_val_0[epoch] += pinball_loss(preds0[:, mask], targets[:, mask], CFG.quantiles[0]).item()\n            epoch_loss_val_1[epoch] += pinball_loss(preds1[:, mask], targets[:, mask], CFG.quantiles[1]).item()\n            epoch_loss_val_2[epoch] += pinball_loss(preds2[:, mask], targets[:, mask], CFG.quantiles[2]).item()\n            \n            # Metric\n            p0 = scaler_fvc.inverse_transform(preds0[:, mask].cpu().detach().numpy())\n            p1 = scaler_fvc.inverse_transform(preds1[:, mask].cpu().detach().numpy())\n            p2 = scaler_fvc.inverse_transform(preds2[:, mask].cpu().detach().numpy())\n            gt = scaler_fvc.inverse_transform(targets[:, mask].cpu().detach().numpy())\n            m[epoch] += (metric(p0[:, :-3], p1[:, :-3], p2[:, :-3], gt[:, :-3]).sum()) # evaluate on last 3 FVC \n            \n        if LR_SCHEDULE == True:\n            schedulerq0.step()\n            schedulerq1.step()\n            schedulerq2.step()\n            \n        # Epoch losses\n        epoch_loss_val_0[epoch] = epoch_loss_val_0[epoch] / dataset_sizes['val']\n        epoch_loss_val_1[epoch] = epoch_loss_val_1[epoch] / dataset_sizes['val']\n        epoch_loss_val_2[epoch] = epoch_loss_val_2[epoch] / dataset_sizes['val']\n        m[epoch] = m[epoch] / dataset_sizes['val']\n            \n        # Update progress bar\n        bar.set_postfix(q0_loss_train=f'{epoch_loss_train_0[epoch]:0.4f}', q1_loss_train=f'{epoch_loss_train_1[epoch]:0.4f}', q2_loss_train=f'{epoch_loss_train_2[epoch]:0.4f}', \\\n                       q0_loss_val=f'{epoch_loss_val_0[epoch]:0.4f}', q1_loss_val=f'{epoch_loss_val_1[epoch]:0.4f}', q2_loss_val=f'{epoch_loss_val_2[epoch]:0.4f}', \\\n                       metric=f'{m[epoch]:0.4}')    \n    # Plot losses \n    fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20,5))\n    \n    ax1.plot(epoch_loss_train_0, color='blue', label='train loss')\n    ax1.plot(epoch_loss_val_0, color='red', label='val loss')\n    handles, labels = ax1.get_legend_handles_labels()\n    ax1.legend(handles, labels, loc='upper left')\n    \n    ax2.plot(epoch_loss_train_1, color='blue', label='train loss')\n    ax2.plot(epoch_loss_val_1, color='red', label='val loss')\n    handles, labels = ax2.get_legend_handles_labels()\n    ax2.legend(handles, labels, loc='upper left')\n    \n    ax3.plot(epoch_loss_train_2, color='blue', label='train loss')\n    ax3.plot(epoch_loss_val_2, color='red', label='train val')\n    handles, labels = ax3.get_legend_handles_labels()\n    ax3.legend(handles, labels, loc='upper left')\n    \n    ax4.plot(m, color='green', label='metric val')\n    handles, labels = ax4.get_legend_handles_labels()\n    ax4.legend(handles, labels, loc='upper left')\n    \n    plt.show()","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}