{"cells":[{"metadata":{},"cell_type":"markdown","source":"The purpose of this notebook is to train a simple Pytorch CNN, which predicts the slope of a linear model based on dcm images. Based on the predicted slopes you can create a very simple linear model, that then predicts for each patient the FVC and the confidence interval. The trained  model does not use any tabular data."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import pydicom\nimport pandas as pd\nimport numpy as np\nfrom pathlib import Path\nimport torch\nimport random\nimport os\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nfrom tqdm.notebook import tqdm\nimport torch.nn.functional as F\nimport torchvision\nimport torchvision.transforms as transforms\nimport cv2\nimport glob, os\nimport re\nimport matplotlib.pyplot as plt\n\nimport gc","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"pd.options.mode.chained_assignment = None\npd.options.display.max_columns = 50","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Path"},{"metadata":{"trusted":true},"cell_type":"code","source":"path = Path('/kaggle/input/osic-pulmonary-fibrosis-progression/')\nassert path.exists()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_TYPES={\"Patient\": \"category\", \n         \"Weeks\": \"int16\", \"FVC\": \"int32\", 'Percent': 'float32', \"Age\": \"uint8\",\n        \"Sex\": \"category\", \"SmokingStatus\": \"category\" }\nSUBMISSION_TYPES={\"Patient_Week\": \"category\", \"FVC\": \"int32\", \"Confidence\": \"int16\"}\n\n\ndef read_data(path):\n    train_df = pd.read_csv(path/'train.csv', dtype = TRAIN_TYPES)\n    test_df = pd.read_csv(path/'test.csv', dtype = TRAIN_TYPES)\n    submission_df = pd.read_csv(path/'sample_submission.csv', dtype = SUBMISSION_TYPES)\n    train_df.drop_duplicates(keep='first', inplace=True, subset=['Patient','Weeks'])\n    return train_df, test_df, submission_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df, test_df, submission_df = read_data(path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def prepare_submission(df, test_df):\n    df['Patient'] = df['Patient_Week'].apply(lambda x:x.split('_')[0])\n    df['Weeks'] = df['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\n    df = df[['Patient','Weeks','Confidence','Patient_Week']]\n    df = df.merge(test_df.drop('Weeks', axis=1).copy(), on=['Patient'])\n    return df\n\nsubmission_df = prepare_submission(submission_df, test_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Feature Engineering"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['WHERE'] = 'train'\ntest_df['WHERE'] = 'val'\nsubmission_df['WHERE'] = 'test'\ndata = train_df.append([test_df, submission_df])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient','FVC', 'Percent']].copy()\nbase.columns = ['Patient','min_FVC', 'min_Percent']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"###### Calculate the slope for the training patients"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import datasets, linear_model\n\nslope_map = {}\nintercept_map = {}\nfor i, p in tqdm(enumerate(train_df.Patient.unique())):\n    sub = train_df.loc[train_df.Patient == p, :] \n    fvc = sub.FVC.values\n    weeks = sub.Weeks.values\n    c = np.vstack([weeks, np.ones(len(weeks))]).T\n    lin_model = linear_model.Ridge()\n    a, b = np.linalg.lstsq(c, fvc, rcond=None)[0]\n    lin_model.fit(X = weeks.reshape([-1, 1]), y = fvc)\n    slope_map[p] = lin_model.coef_[0]\n    intercept_map[p] = np.log(lin_model.intercept_) * -1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create datasets"},{"metadata":{"trusted":true},"cell_type":"code","source":"BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_img(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize(d.pixel_array / 2**11, IMG_DIM)\ntry:\n    get_img(f'{path}/train/ID00219637202258203123958/9.dcm').shape\nexcept:\n    print('Image not found')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_image_slope(ds_type = 'train'):\n    img_folders = [t[0] for t in os.walk(path/f\"{ds_type}\") if re.match(r'.+ID\\d+$', t[0])]\n    training_data = []\n    folder_count = {}\n    for f in img_folders:\n        patient = re.sub(r'.+/(.+)', r'\\1', f)\n        if patient not in BAD_ID:\n            slope = slope_map[patient]\n            ldir = glob.glob(str(f'{f}/*.dcm'))\n            folder_count[patient] = len(ldir)\n            for img_file in ldir:\n                try:\n                    file_name = re.sub(r'.+/(.+)', r'\\1', img_file)\n                    if re.match(r'\\d+\\..*', file_name):\n                        if int(file_name[:-4]) / len(ldir) < 0.8 and int(file_name[:-4]) / len(ldir) > 0.15:\n                            training_data.append((img_file, slope))\n                except:\n                    print(f'Failed on f{img_file}')\n    return training_data, folder_count","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_slope_data, folder_count = read_image_slope()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create training and validation datasets"},{"metadata":{"trusted":true},"cell_type":"code","source":"IMG_DIM = (256, 256)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split \n\ntrain_image_slope, val_image_slope = train_test_split(img_slope_data, shuffle=True, train_size= 0.9) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"###### Getting the mean and standard deviation on the images"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_mean_std(ds):\n    # VAR[X] = E[X**2] - E[X]**2\n    channels_sum, channels_squared_sum, num_batches = 0, 0, 0\n    \n    for data, _ in ds:\n        channels_sum += torch.mean(data, dim=[-2, -1])\n        channels_squared_sum += torch.mean(data**2, dim=[-2, -1])\n        num_batches += 1\n    \n    mean = channels_sum / num_batches\n    std = (channels_squared_sum / num_batches - mean ** 2) ** 0.5\n    return mean, std","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# mean,std = get_mean_std(CombinedImageDataset(train_image_slope))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mean = torch.tensor([-0.0149])\nstd = torch.tensor([0.5000])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class CombinedImageDataset(Dataset):\n    def __init__(self, img_slope_data, transforms=transforms.Compose([transforms.ToTensor()])):\n        self.data = img_slope_data\n        self.transforms = transforms\n    \n    def __len__(self):\n        return len(self.data)\n    \n    def __getitem__(self, i):\n        path, slope_intercept = self.data[i]\n        slope= slope_intercept\n        image = get_img(path)\n        if self.transforms is not None:\n            image = self.transforms(image)\n        return image.float(), slope\n    \n    def get_patient(self, i):\n        path = self.data[i][0]\n        return re.sub(r'.+/(ID.+?)/.+', r'\\1', path)\n        \n    def __repr__(self):\n        return  f'patients: {len(self.data)}, image_path: {self.image_path}, transforms: {self.transforms}'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize(mean,std)])\ntest_transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize(mean,std)])\nvalid_transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize(mean,std)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_ds = CombinedImageDataset(train_image_slope, transforms=train_transform)\nlen(train_ds)\nval_ds = CombinedImageDataset(val_image_slope, transforms=test_transform)\nlen(train_ds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 64\nNUM_WORKERS = 4\n\ntrain_dl = DataLoader(train_ds, BATCH_SIZE, shuffle=False, num_workers=NUM_WORKERS)\nval_dl = DataLoader(val_ds, BATCH_SIZE, shuffle=False, num_workers=NUM_WORKERS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def move_to_dev(items):\n    return [i.to(device) for i in items]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def freeze(model, layers=7, requires_grad=False):\n    ct = 0\n    for name, child in model.named_children():\n        ct += 1\n        if ct < layers:\n            for _, params in child.named_parameters():\n                params.requires_grad = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def unfreeze_all(model):\n    # Unfreeze model weights\n    for param in model.parameters():\n        param.requires_grad = True","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_model():\n    model = torchvision.models.resnet.ResNet(torchvision.models.resnet.Bottleneck, [3, 4, 6, 3])\n    model.load_state_dict(torch.load('/kaggle/input/resnet50/resnet50.pth'))\n    num_ftrs = model.fc.in_features\n    model.conv1 = nn.Conv2d(1, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)\n    model.fc = nn.Linear(num_ftrs, 1)\n    freeze(model)\n    model = model.to(device);\n    return model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##### Test the model"},{"metadata":{"trusted":true},"cell_type":"code","source":"model1 = create_model()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_img, sample_slope = move_to_dev(next(iter(train_dl)))\nsample_img.shape, sample_slope","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_out = model1(sample_img)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slope_mse = torch.mean((sample_out[:,0] - sample_slope) ** 2)\nslope_mse","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Training"},{"metadata":{"trusted":true},"cell_type":"code","source":"LR=1e-4\ncriterion = nn.MSELoss()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_lr(optimizer):\n    for param_group in optimizer.param_groups:\n        return param_group['lr']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def eval_loop(valid_dl, model):\n    with torch.no_grad():\n        model.eval()\n        total_eval_loss = 0\n        total_eval_score = 0\n        for val_vals in valid_dl:\n            x, y_slope = move_to_dev(val_vals)\n            output = model(x)\n            loss = criterion(y_slope.unsqueeze(-1), output)\n            total_eval_loss += loss.item()\n            total_eval_score += loss.item()\n\n        avg_val_loss = total_eval_loss / len(valid_dl)\n        avg_val_score = total_eval_score / len(valid_dl)\n        return {\n            'avg_val_loss': avg_val_loss,\n            'avg_val_score': avg_val_score\n        }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_loop(epochs, train_dl, valid_dl, model, lr = 1e-3, print_score=False, model_name='test'):\n    steps = len(train_dl) * epochs\n    optimizer = torch.optim.Adam(model.parameters(), lr=lr)\n    scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=lr, steps_per_epoch=len(train_dl), epochs=epochs)\n    avg_train_losses = []\n    avg_val_losses = []\n    avg_val_scores = []\n    lr = []\n    best_avg_val_score = 1000\n    for epoch in tqdm(range(epochs), total=epochs):\n        model.train()\n        total_train_loss = 0.0\n        t = tqdm(enumerate(train_dl), total=len(train_dl))\n        for i, train_vals in t:\n            x, y_slope = move_to_dev(train_vals)\n            model.zero_grad()\n            output = model(x)\n            loss = criterion(y_slope.unsqueeze(-1), output)\n            t.set_postfix({'loss': loss.item()})\n            total_train_loss += loss.item()\n            \n            # Backward Pass and Optimization\n            loss.backward()\n            optimizer.step()\n            scheduler.step()\n            lr.append(get_lr(optimizer))\n        \n        avg_train_loss = total_train_loss / len(train_dl)\n        avg_train_losses.append(avg_train_loss)\n        eval_res = eval_loop(valid_dl, model)\n        avg_val_loss = eval_res['avg_val_loss']\n        avg_val_score = eval_res['avg_val_score']\n        avg_val_losses.append(avg_val_loss)\n        avg_val_scores.append(avg_val_score)\n        if best_avg_val_score > avg_val_score:\n            best_avg_val_score = avg_val_score\n            # save best model\n            print(f'Best model: {best_avg_val_score}')\n#             torch.save(model.state_dict(), model_path/f'best_model_images_{model_name}.pt')\n        if print_score:\n            print(f'{epoch}: avg_val_score: {avg_val_score}')\n    return pd.DataFrame({'avg_train_losses': avg_train_losses, 'avg_val_losses': avg_val_losses, 'avg_val_scores': avg_val_scores}), pd.DataFrame({'lr': lr})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NUM_EPOCHS = 2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"res_df, lr_df = train_loop(NUM_EPOCHS, train_dl, val_dl, model1, lr=LR, print_score=True)\nres_df[['avg_train_losses', 'avg_val_losses']].plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"UNFROZEN_EPOCHS=2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"unfreeze_all(model1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"res_df, lr_df = train_loop(UNFROZEN_EPOCHS, train_dl, val_dl, model1, lr=LR / 10, print_score=True)\nres_df[['avg_train_losses', 'avg_val_losses']].plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_orig_df = pd.read_csv(path/'test.csv', dtype = TRAIN_TYPES)\nsub_df = pd.read_csv(path/'sample_submission.csv', dtype = SUBMISSION_TYPES)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"current_folder = 'test'\n\ndef load_data(current_df, current_folder, ldir, p, x):\n    for i in ldir:\n        if int(i[:-4]) / len(ldir) < 0.8 and int(i[:-4]) / len(ldir) > 0.15:\n            x.append(get_img(f'/kaggle/input/osic-pulmonary-fibrosis-progression/{current_folder}/{p}/{i}'))\n\n\nquantiles = [0.2, 0.5, 0.8]\nunique_patients = test_orig_df.Patient.unique()\n\nwith torch.no_grad():\n    model1.eval()\n    A_test, B_test, P_test, WEEK = {}, {}, {},{}\n    for p in unique_patients:\n        x = []\n        ldir = os.listdir(f'/kaggle/input/osic-pulmonary-fibrosis-progression/{current_folder}/{p}/')\n        load_data(test_orig_df, current_folder, ldir, p, x)\n        if len(x) <= 1:\n            continue\n\n        x = np.expand_dims(x, axis=-1)\n        _a = model1(torch.tensor(x).squeeze().unsqueeze(1).float().to(device))\n\n        # A = slopes, B = intercepts\n        A_test[p] = [np.quantile(_a.cpu(), q) for q in quantiles]\n        B_test[p] = test_orig_df.FVC.values[test_orig_df.Patient == p] - A_test[p] * test_orig_df.Weeks.values[test_orig_df.Patient == p]\n        P_test[p] = test_orig_df.Percent.values[test_orig_df.Patient == p]\n        WEEK[p] = test_orig_df.Weeks.values[test_orig_df.Patient == p]\n\ndef fvc_calc(a, x, b): return a * x + b\n        \nfor k in sub_df.Patient_Week.values:\n    p, w = k.split('_')\n    w = int(w)\n\n    fvc = fvc_calc(A_test[p][1], w, B_test[p][1])\n    sub_df.loc[sub_df.Patient_Week == k, 'FVC'] = fvc\n    conf_1 = np.abs(P_test[p] - A_test[p][1] * abs(WEEK[p] - w))\n    conf_2 = np.abs(fvc_calc(A_test[p][2], w, B_test[p][2]) - fvc_calc(A_test[p][0], w, B_test[p][0]))\n    sub_df.loc[sub_df.Patient_Week == k, 'Confidence'] = np.clip(np.average([conf_1, conf_2], axis=0), 100, 1000)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_df[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_final_df = pd.read_csv(\"submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for p in test_df['Patient'].unique():\n    submission_final_df[submission_final_df['Patient_Week'].str.find(p) == 0]['FVC'].plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}