{"cells":[{"metadata":{},"cell_type":"markdown","source":"Borrowed pieces from the excellent kernel https://www.kaggle.com/carlossouza/quantile-regression-pytorch-tabular-data-only. Some people had followed my other similar work with the bayesian neural network, tabular data estimating mu and sigma and dropout as bayesian estimation so I've done the quantile regression in approximately the same framework.\n\nI'm working on a neat kernel that combines all of these into one watch this space.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n#Load the dependancies\nfrom tqdm.notebook import tqdm\nimport copy\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport torch\nimport torch.nn as nn\nfrom torch import optim\nimport random\nfrom torchvision import models\nimport torch.multiprocessing as mp\nfrom torch.nn import functional as F\nimport matplotlib.pyplot as plt\nplt.style.use('ggplot')\nimport pydicom.pixel_data_handlers.gdcm_handler as gdcm_handler \nimport cv2\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nimport pydicom\nimport os\nfrom torch.utils.data import DataLoader, Dataset\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n#    for filename in filenames:\n#        print(os.path.join(dirname, filename))\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_all(seed=42):\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = False\n    \nseed_all()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_metric_loss(pred_fvc,sigma,true_fvc):\n    true_fvc=np.reshape(true_fvc,pred_fvc.shape)\n    sigma[sigma<70]=70\n    delta=np.abs(pred_fvc-true_fvc)\n    delta[delta>1000]=1000\n    metric=-(np.sqrt(2)*delta/sigma)-np.log(np.sqrt(2)*sigma)\n    return -metric","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class main_model(nn.Module):\n    def __init__(self, n_additional_features, n_outputs):\n        super(main_model, self).__init__()\n        self.fc1 = nn.Linear(n_additional_features, 100)\n        self.fc2 = nn.Linear(100, 100)\n        self.fc3 = nn.Linear(100, n_outputs)\n\n    def forward(self, additional_features):\n        out=additional_features\n        out = F.relu(self.fc1(out))\n        out = F.relu(self.fc2(out))\n        out = self.fc3(out)\n        return out\n    \n    def metric_loss(self,pred_fvc,true_fvc):\n        sigma = pred_fvc[:, 2] - pred_fvc[:, 0]\n        true_fvc=torch.reshape(true_fvc,pred_fvc[:,1].shape)\n        sigma_clipped=torch.clamp(sigma,min=70)\n        delta=torch.clamp(torch.abs(pred_fvc[:,1]-true_fvc),max=1000)\n        metric=torch.div(-torch.sqrt(torch.tensor([2.0]).to(device))*delta,sigma_clipped)-torch.log(torch.sqrt(torch.tensor([2.0]).to(device))*sigma_clipped)\n        return -metric\n    \n    def fvc_loss(self,pred_fvc,true_fvc):\n        true_fvc=torch.reshape(true_fvc,pred_fvc[:,1].shape)\n        fvc_err=torch.abs(pred_fvc-true_fvc)\n        return fvc_err\n    \n    def quantile_loss(self,preds, target, quantiles):\n        assert not target.requires_grad\n        assert preds.size(0) == target.size(0)\n        losses = []\n        for i, q in enumerate(quantiles):\n            errors = target - preds[:, i]\n            losses.append(torch.max((q - 1) * errors, q * errors).unsqueeze(1))\n\n        loss = torch.mean(torch.sum(torch.cat(losses, dim=1), dim=1))\n        return loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_training_loss(train, val,title='loss'):\n    plt.figure()\n    plt.plot(train, label='Train')\n    plt.plot(val, label='Val')\n    if title=='loss':\n        plt.title('Model Training Loss')\n    else:\n        plt.title('Model Metric Loss')\n    plt.xlabel('Epoch')\n    plt.ylabel('Loss')\n    plt.legend()\n    plt.savefig('training_loss')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Dataframes","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nsubmission=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prepare Training Data (Tabular)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train['base_Weeks']=train.groupby(['Patient'])['Weeks'].transform('min')\nbase=train[train.Weeks==train.base_Weeks]\nbase = base.rename(columns={'FVC': 'base_FVC','Percent': 'base_Percent'})\nbase.drop_duplicates(subset=['Patient', 'Weeks'], keep='first',inplace=True)\ntrain=train.merge(base[['Patient','base_FVC','base_Percent']],on='Patient',how='left')\ntrain['Week_passed'] = train['Weeks'] - train['base_Weeks']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prepare Test Data (tabular)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test = test.rename(columns={'Weeks': 'base_Weeks', 'FVC': 'base_FVC','Percent': 'base_Percent'})\n\n# Adding Sample Submission\nsubmission = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/sample_submission.csv\")\n\n# In submisison file, format: ID_'week', using lambda to split the ID\nsubmission['Patient'] = submission['Patient_Week'].apply(lambda x:x.split('_')[0])\n\n# In submisison file, format: ID_'week', using lambda to split the Week\nsubmission['Weeks'] = submission['Patient_Week'].apply(lambda x:x.split('_')[1]).astype(int)\n\ntest = submission.drop(columns = [\"FVC\", \"Confidence\"]).merge(test, on = 'Patient')\n\ntest['Week_passed'] = test['Weeks'] - test['base_Weeks']\n\ntest=test[train.columns.drop(['FVC','Percent'])]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# OH Encode Sex and Smoking\nWith thanks to https://www.kaggle.com/ulrich07/osic-keras-starter-with-custom-metrics","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"COLS = ['Sex','SmokingStatus']\nfor col in COLS:\n    for mod in train[col].unique():\n        train[mod] = (train[col] == mod).astype(int)\n        \n        test[mod] = (test[col] == mod).astype(int)\n    train.drop(col,axis=1,inplace=True)\n    test.drop(col,axis=1,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Rescale based on train data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\nrobust_scaler = preprocessing.RobustScaler()\ntrain[train.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])]=robust_scaler.fit_transform(train[train.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# I've used this to avoid the files dcmread can't load (future version with CNN)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pickle\nwith open('../input/file-dictionary/good_files.pickle', 'rb') as handle:\n    good_file_dict = pickle.load(handle)\nwith open('../input/file-dictionary/bad_files.pickle', 'rb') as handle:\n    bad_file_dict = pickle.load(handle)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class OSIC(Dataset):\n    def __init__(self,patient_ids,df,file_dict,train=True, transform=None,nims=10):\n        self.df=df[df.Patient.isin(patient_ids)]\n        self.train=train\n        if self.train:\n            self.fvc=self.df['FVC'].values\n        else:\n            self.df[self.df.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])]=robust_scaler.transform(self.df[self.df.columns.difference(['Patient','FVC','Percent','Weeks','base_Weeks'])])\n        self.data=self.df[self.df.columns.difference(['FVC','Patient','Percent'])].values\n        self.patients=self.df['Patient'].values\n        self.file_dict=file_dict\n        self.nims=nims\n        \n    def __len__(self):\n        return len(self.df)\n\n    def __getitem__(self, idx):\n        if self.train:\n            data = {'fvc': self.fvc[idx],\n                   'data': self.data[idx]}\n        else:\n            data = {'data': self.data[idx]}\n        return data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training Parameters","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"epochs=400\nbatch_size=32\nnum_workers=3\nquantiles = (0.159, 0.5, 0.841)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Split training data into train and val by patient (80:20)\n'ID00011637202177653955184' has no images we can load. Therefore I'm going to drop.\n\nWe shuffle the train data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ids=train.Patient.unique()\nindex = np.argwhere(ids=='ID00011637202177653955184')\nids = list(np.delete(ids, index))\nrandom.shuffle(ids)\nids=np.array(ids)\n\ntrain_ids,val_ids=np.split(ids, [int(round(0.9 * len(ids), 0))])\n\ntrain_dataset = OSIC(train_ids,train,good_file_dict)  \ntrain_dataloader = DataLoader(train_dataset, batch_size=batch_size,shuffle=True,num_workers=num_workers)\n\nval_dataset = OSIC(val_ids,train,good_file_dict)  \nval_dataloader = DataLoader(val_dataset, batch_size=batch_size,shuffle=True,num_workers=num_workers)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = main_model(train_dataset.data.shape[1], 3).to(device)\nprint('Number of parameters:')\nprint(sum(p.numel() for p in model.parameters() if p.requires_grad))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Pre-train with mse of FVC prediction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"optimizer = optim.Adam(model.parameters())\nepoch_train_metric=[]\nepoch_val_metric=[]\nepoch_train_loss=[]\nepoch_val_loss=[]\nepoch=0\nmin_val_loss = 1e+100\nearly_stop = False\npatience=10\n#Start by training for fvc\nwhile epoch<epochs and not early_stop:\n    epoch+=1\n    train_loss=0\n    train_metric=0\n    model.train()\n    for batch_idx, data in enumerate(train_dataloader):\n        optimizer.zero_grad()\n        model_output = model(data['data'].float().to(device))\n        quantile_loss = model.quantile_loss(model_output, data['fvc'].to(device), quantiles)\n        metric_loss = model.metric_loss(model_output,data['fvc'].to(device)).mean()\n        loss=quantile_loss\n        loss.backward()\n        train_loss += loss.item()\n        train_metric += metric_loss.item()\n        optimizer.step()\n    print('====> Epoch: {} Average train quantile loss: {:.4f}'.format(\n                        epoch, train_loss / len(train_dataloader)))\n    print('====> Epoch: {} Average train metric: {:.4f}'.format(\n                        epoch, train_metric / len(train_dataloader)))\n    \n    val_loss=0\n    val_metric=0\n    model.eval()\n    with torch.no_grad():\n        for batch_idx, data in enumerate(val_dataloader):\n            model_output = model(data['data'].float().to(device))\n            quantile_loss = model.quantile_loss(model_output, data['fvc'].to(device), quantiles)\n            metric_loss = model.metric_loss(model_output,data['fvc'].to(device)).mean()\n            loss=quantile_loss\n            val_loss += loss.item()\n            val_metric += metric_loss.item()\n        print('====> Epoch: {} Average val quantile loss: {:.4f}'.format(\n                        epoch, val_loss / len(val_dataloader)))\n        print('====> Epoch: {} Average val metric: {:.4f}'.format(\n                        epoch, val_metric / len(val_dataloader)))\n    \n    epoch_train_loss.append(train_loss/ len(train_dataloader))\n    epoch_val_loss.append(val_loss / len(val_dataloader))\n    epoch_train_metric.append(train_metric/ len(train_dataloader))\n    epoch_val_metric.append(val_metric / len(val_dataloader))\n    \n    if val_loss < min_val_loss:\n        min_val_loss = val_loss\n        best_model = copy.deepcopy(model.state_dict())\n        print('Min loss %0.2f' % min_val_loss)\n        epochs_no_improve = 0\n\n    else:\n        epochs_no_improve += 1\n        # Check early stopping condition\n        if epochs_no_improve == patience:\n            print('Early stopping!')\n            early_stop = True\n            model.load_state_dict(best_model)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Plot training curves","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_training_loss(epoch_train_loss, epoch_val_loss)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"possibly a learning rate problem/something to do with the different scaling of percent and fvc","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_training_loss(epoch_train_metric, epoch_val_metric,title='metric')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission=pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_ids=test.Patient.unique()\ntest_dataset = OSIC(test_ids,test,good_file_dict,train=False)  \ntest_dataloader = DataLoader(test_dataset, batch_size=batch_size,shuffle=False,num_workers=num_workers)\nmodel.eval()\nfvc_pred = []\nsigma_pred = []\nwith torch.no_grad():\n    for batch_idx, data in enumerate(test_dataloader):\n        model_output = model(data['data'].float().to(device))\n        fvc_pred.append(model_output[:,1])\n        sigma_pred.append(model_output[:,2]-model_output[:,0])\nfvc_pred=torch.cat(fvc_pred, dim=0)\nsigma_pred=torch.cat(sigma_pred, dim=0)\ntest['FVC']=fvc_pred.cpu().numpy()\ntest['Confidence']=sigma_pred.cpu().numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test['Patient_Week']=test[\"Patient\"] + '_' + test['Weeks'].apply(str)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission=submission[['Patient_Week']].merge(test[['Patient_Week','FVC','Confidence']],on='Patient_Week')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False, float_format='%.1f')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Post-Match Analysis","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Test Predictions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(submission['FVC'],submission['Confidence'])\nplt.title('Test')\nplt.xlabel('FVC')\nplt.ylabel('Confidence')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fvc_pred = []\nsigma_pred = []\ntrain_dataloader = DataLoader(train_dataset, batch_size=batch_size,shuffle=False,num_workers=num_workers)\nwith torch.no_grad():\n    for batch_idx, data in enumerate(train_dataloader):\n        model_output = model(data['data'].float().to(device))\n        fvc_pred.append(model_output[:,1])\n        sigma_pred.append(model_output[:,2]-model_output[:,0])\nfvc_pred_train=torch.cat(fvc_pred, dim=0)\nsigma_pred_train=torch.cat(sigma_pred, dim=0)\n\nprint('train metric', calc_metric_loss(fvc_pred_train.cpu().numpy(),sigma_pred_train.cpu().numpy(),train_dataset.fvc).mean())\n\nplt.scatter(fvc_pred_train.cpu().numpy(),sigma_pred_train.cpu().numpy())\nplt.title('Train')\nplt.xlabel('FVC')\nplt.ylabel('Confidence')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(train_dataset.fvc,fvc_pred_train.cpu().numpy())\nplt.title('Train: predicted FVC vs true FVC')\nplt.xlabel('True FVC')\nplt.ylabel('Predicted FVC')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Val","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fvc_pred = []\nsigma_pred = []\nval_dataloader = DataLoader(val_dataset, batch_size=batch_size,shuffle=False,num_workers=num_workers)\nwith torch.no_grad():\n    for batch_idx, data in enumerate(val_dataloader):\n        model_output = model(data['data'].float().to(device))\n        fvc_pred.append(model_output[:,1])\n        sigma_pred.append(model_output[:,2]-model_output[:,0])\nfvc_pred_val=torch.cat(fvc_pred, dim=0)\nsigma_pred_val=torch.cat(sigma_pred, dim=0)\n\nprint('val metric', calc_metric_loss(fvc_pred_val.cpu().numpy(),sigma_pred_val.cpu().numpy(),val_dataset.fvc).mean())\n\nplt.scatter(fvc_pred_val.cpu().numpy(),sigma_pred_val.cpu().numpy())\nplt.title('Val')\nplt.xlabel('FVC')\nplt.ylabel('Confidence')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(val_dataset.fvc,fvc_pred_val.cpu().numpy())\nplt.title('Val: predicted FVC vs true FVC')\nplt.xlabel('True FVC')\nplt.ylabel('Predicted FVC')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## All","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(submission['FVC'], alpha=0.5,label='test')\nplt.hist(fvc_pred_train.cpu().numpy(), alpha=0.5,label='train')\nplt.hist(fvc_pred_val.cpu().numpy(), alpha=0.5,label='val')\nplt.legend()\nplt.title('Histogram of FVC predictions')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(submission['Confidence'], alpha=0.5,label='test')\nplt.hist(sigma_pred_train.cpu().numpy(), alpha=0.5,label='train')\nplt.hist(sigma_pred_val.cpu().numpy(), alpha=0.5,label='val')\nplt.legend()\nplt.title('Histogram of Confidence predictions')","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}