{"cells":[{"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\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\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\n\nimport os\nfor 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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport random\nimport math\nimport gc\nimport glob\n\nfrom sklearn.model_selection import KFold, GroupKFold, train_test_split\nfrom sklearn.cluster import KMeans\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_squared_error\n\n\nimport operator\nimport typing as tp\nfrom logging import getLogger, INFO, StreamHandler, FileHandler, Formatter\nfrom functools import partial\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom tqdm.notebook import tqdm\n\nimport pydicom\n\nprint(os.listdir('/kaggle/input/osic-pulmonary-fibrosis-progression/'))\n\nfrom time import time, strftime, gmtime\n\nstart = time()\n#print(start)\n\nimport datetime\nprint(str(datetime.datetime.now()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seed = 2019\nseed_everything(seed)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_path = '/kaggle/input/osic-pulmonary-fibrosis-progression/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(input_path + 'train.csv')\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(input_path + 'test.csv')\ntest","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample = pd.read_csv(input_path + 'sample_submission.csv')\nsample","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Number of patients in the train set: ', train['Patient'].nunique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\n\nplt.title('Age')\nsns.distplot(train['Age'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\n\nplt.title('Weeks')\nsns.distplot(train['Weeks'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\n\nplt.title('FVC')\nsns.distplot(train['FVC'], bins = int(1 + math.log2(train['Patient'].nunique())))\n#use Sturgess's formula to find the appropriate number of classes in the histogram","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Max. FVC: ', train['FVC'].max())\nprint('Min. FVC: ', train['FVC'].min())\nprint('Mean FVC: ', train['FVC'].mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['SmokingStatus'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize = (15, 15))\n#plt.figure(figsize = (10, 10))\nplt.suptitle('Smoking Status')\nsns.countplot(train['SmokingStatus'], ax = ax[0])\n\nlbls, freqs = np.unique(train['SmokingStatus'].values, return_counts = True)\n\nax[1].pie(freqs, labels = lbls, autopct = '%1.1f%%', shadow = False, startangle = 90)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize = (15, 15))\n\nplt.suptitle('Sex')\nsns.countplot(train['Sex'], ax = ax[0])\n\nlbls, freqs = np.unique(train['Sex'].values, return_counts = True)\n\nax[1].pie(freqs, labels = lbls, autopct = '%1.1f%%', shadow = False, startangle = 90)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\nplt.title('Smoking Status by Sex')\nsns.countplot(train['SmokingStatus'], hue = train['Sex'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\nplt.title('Smoking Status by Sex')\nsns.countplot(train['Age'], hue = train['SmokingStatus'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\nsns.heatmap(train.corr(), annot = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Correlation coeff between Age and FVC is: ', train.corr()['Age']['FVC'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (8, 8))\nsns.scatterplot(data = train, x = 'Age', y = 'FVC')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Corr for smokers\ntrain_cs = train.loc[train['SmokingStatus'] == 'Currently smokes']\n\nplt.figure(figsize = (8, 8))\nsns.scatterplot(data = train_cs, x = 'Age', y = 'FVC')\n\nprint('Correlation coeff between Age and FVC (Current Smokers) is: ', train_cs.corr()['Age']['FVC'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dcm = glob.glob(input_path + 'train/*/*')\ntest_dcm = glob.glob(input_path + 'test/*/*')\n\nprint('Num of train dicom: ', len(train_dcm))\nprint('Num of test dicom: ', len(test_dcm))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_imgs_pid = [len(os.listdir(input_path + 'train/' + path)) for path in os.listdir(input_path + 'train/')]\nplt.figure(figsize = (8, 8))\nplt.hist(num_imgs_pid)\nplt.ylabel('Number of patients')\nplt.xlabel('DICOM files')\nplt.title('DICOM Images per patient')\nplt.show()\nprint('Max. number of dicom images per patient: ', max(num_imgs_pid))\nprint('Min. number of dicom images per patient: ', min(num_imgs_pid))\nprint('Mean. number of dicom images per patient: ', np.mean(num_imgs_pid))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pydicom.dcmread(train_dcm[10])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"random_dcm_train = np.random.choice(train_dcm, 6)\nfig, ax = plt.subplots(2, 3, figsize = (15, 10))\n\nax = ax.ravel()\n\nfor i, file in enumerate(random_dcm_train):\n    img = pydicom.dcmread(file)\n    img = img.pixel_array\n    # Since the scanning equipment is cylindrical in nature and image output is square,\n    # we set the out-of-scan pixels to 0\n    img[img == -2000] = 0\n    ax[i].imshow(img, cmap = plt.cm.bone)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2, 3, figsize = (15, 10))\n\nax = ax.ravel()\n\nfor i, file in enumerate(random_dcm_train):\n    img = pydicom.dcmread(file)\n    img = img.pixel_array\n    img[img == -2000] = 0\n    ax[i].imshow(img, cmap = plt.cm.Reds)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_dicom_images(pid = None, df = None, feature = None, dcm_path = None):\n    fig, ax = plt.subplots(2, 3, figsize = (15, 10))\n    ax = ax.ravel()\n    for i in range(len(ax)):\n        dcm = pydicom.dcmread(input_path + '/train/' + pid + '/' + str(i + 1) + '.dcm')\n        img = dcm.pixel_array\n        img[img == -2000] = 0\n        ax[i].imshow(img, cmap = plt.cm.bone)\n    return None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_dicom_images(pid = np.random.choice(train['Patient'].values, 1)[0], df = train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['Patient_Week'] = train['Patient'].astype(str) + '_' + train['Weeks'].astype(str)\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_out = pd.DataFrame()\n\ntrain_grp = train.groupby('Patient')\n\nfor _, df_out in tqdm(train_grp):\n    df_pid = pd.DataFrame()\n    for wk, temp in df_out.groupby('Weeks'):\n        rename_cols = {'Weeks': 'base_Week', 'FVC': 'base_FVC', 'Percent': 'base_Percent', 'Age': 'base_Age'}\n        temp = temp.drop(columns = 'Patient_Week').rename(columns = rename_cols)\n        drop_cols = ['Age', 'Sex', 'SmokingStatus', 'Percent']\n        _df_pid = df_out.drop(columns = drop_cols).rename(columns = {'Weeks': 'predict_Week'}).merge(temp, on = 'Patient')\n        _df_pid['Week_passed'] = _df_pid['predict_Week'] - _df_pid['base_Week']\n        df_pid = pd.concat([df_pid, _df_pid])\n    train_out = pd.concat([train_out, df_pid])\n\ntrain = train_out[train_out['Week_passed'] != 0].reset_index(drop = True)\ntrain","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = test.rename(columns = {'Weeks': 'base_Week', 'FVC': 'base_FVC', 'Percent': 'base_Percent', 'Age': 'base_Age'})\n\nsample['Patient'] = sample['Patient_Week'].apply(lambda x: x.split('_')[0])\nsample['predict_Week'] = sample['Patient_Week'].apply(lambda x: x.split('_')[1]).astype(int)\ntest = sample.drop(columns = ['FVC', 'Confidence']).merge(test, on = 'Patient')\ntest['Week_passed'] = test['predict_Week'] - test['base_Week']\ntest","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"folds = train[['Patient_Week', 'Patient', 'FVC']].copy()\ngkf = GroupKFold(n_splits = 5)\ngroups = folds['Patient'].values\nfor i, (trn_idx, val_idx) in enumerate(gkf.split(folds, folds['FVC'], groups)):\n    folds.loc[val_idx, 'fold'] = int(i)\nfolds['fold'] = folds['fold'].astype(int)\nfolds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgbm\n\nparams = {\n    'num_class': 2,\n    'metric': 'None',\n    'boosting_type': 'gbdt',\n    'learning_rate': 5e-02,\n    'seed': seed,\n    \"subsample\": 0.4,\n    \"subsample_freq\": 1,\n    'max_depth': 1,\n    'verbosity': -1,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class OSICLossForLGBM:\n    \"\"\"\n    Custom Loss for LightGBM.\n    \n    * Objective: return grad & hess of NLL of gaussian\n    * Evaluation: return competition metric\n    \"\"\"\n    \n    def __init__(self, epsilon: float=1) -> None:\n        \"\"\"Initialize.\"\"\"\n        self.name = \"osic_loss\"\n        self.n_class = 2  # FVC & Confidence\n        self.epsilon = epsilon\n    \n    def __call__(self, preds: np.ndarray, labels: np.ndarray, weight: tp.Optional[np.ndarray]=None) -> float:\n        \"\"\"Calc loss.\"\"\"\n        sigma_clip = np.maximum(preds[:, 1], 70)\n        Delta = np.minimum(np.abs(preds[:, 0] - labels), 1000)\n        loss_by_sample = - np.sqrt(2) * Delta / sigma_clip - np.log(np.sqrt(2) * sigma_clip)\n        loss = np.average(loss_by_sample, weight)\n        \n        return loss\n    \n    def _calc_grad_and_hess(\n        self, preds: np.ndarray, labels: np.ndarray, weight: tp.Optional[np.ndarray]=None\n    ) -> tp.Tuple[np.ndarray]:\n        \"\"\"Calc Grad and Hess\"\"\"\n        mu = preds[:, 0]\n        sigma = preds[:, 1]\n        \n        sigma_t = np.log(1 + np.exp(sigma))\n        grad_sigma_t = 1 / (1 + np.exp(- sigma))\n        hess_sigma_t = grad_sigma_t * (1 - grad_sigma_t)\n        \n        grad = np.zeros_like(preds)\n        hess = np.zeros_like(preds)\n        grad[:, 0] = - (labels - mu) / sigma_t ** 2\n        hess[:, 0] = 1 / sigma_t ** 2\n        \n        tmp = ((labels - mu) / sigma_t) ** 2\n        grad[:, 1] = 1 / sigma_t * (1 - tmp) * grad_sigma_t\n        hess[:, 1] = (\n            - 1 / sigma_t ** 2 * (1 - 3 * tmp) * grad_sigma_t ** 2\n            + 1 / sigma_t * (1 - tmp) * hess_sigma_t\n        )\n        if weight is not None:\n            grad = grad * weight[:, None]\n            hess = hess * weight[:, None]\n        return grad, hess\n    \n    def return_loss(self, preds: np.ndarray, data: lgbm.Dataset) -> tp.Tuple[str, float, bool]:\n        \"\"\"Return Loss for lightgbm\"\"\"\n        labels = data.get_label()\n        weight = data.get_weight()\n        n_example = len(labels)\n        \n        # # reshape preds: (n_class * n_example,) => (n_class, n_example) =>  (n_example, n_class)\n        preds = preds.reshape(self.n_class, n_example).T\n        # # calc loss\n        loss = self(preds, labels, weight)\n        \n        return self.name, loss, True\n    \n    def return_grad_and_hess(self, preds: np.ndarray, data: lgbm.Dataset) -> tp.Tuple[np.ndarray]:\n        \"\"\"Return Grad and Hess for lightgbm\"\"\"\n        labels = data.get_label()\n        weight = data.get_weight()\n        n_example = len(labels)\n        \n        # # reshape preds: (n_class * n_example,) => (n_class, n_example) =>  (n_example, n_class)\n        preds = preds.reshape(self.n_class, n_example).T\n        # # calc grad and hess.\n        grad, hess =  self._calc_grad_and_hess(preds, labels, weight)\n\n        # # reshape grad, hess: (n_example, n_class) => (n_class, n_example) => (n_class * n_example,) \n        grad = grad.T.reshape(n_example * self.n_class)\n        hess = hess.T.reshape(n_example * self.n_class)\n        \n        return grad, hess","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_features = ['Sex', 'SmokingStatus']\nnum_features = [col for col in train.columns if (train[col].dtype != 'object') & (col not in cat_features)]\nprint(cat_features, num_features)\nfeatures = cat_features + num_features\nfeatures = [col for col in features if col not in ['Patient_Week', 'FVC', 'predict_Week', 'base_Week']]\nfeatures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import category_encoders as catenc\n\ntest['FVC'] = np.nan\n\nordenc = catenc.OrdinalEncoder(cols = cat_features, handle_unknown = 'impute')\nordenc.fit(train)\ntrain = ordenc.transform(train)\ntest = ordenc.transform(test)\nprint('Categorical features encoded..')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nb_splits = 5\noof = np.zeros((len(train), 2))\npredictions = np.zeros((len(test), 2))\nfeature_importance_df = pd.DataFrame()\nosic_loss = OSICLossForLGBM()\n\nfor n_folds in range(nb_splits):\n    print()\n    print('Fold No: ', n_folds + 1)\n    trn_idx = folds[folds['fold'] != n_folds].index\n    val_idx = folds[folds['fold'] == n_folds].index\n    #print(trn_idx, val_idx)\n    ltrain = lgbm.Dataset(train.iloc[trn_idx][features], label = train.iloc[trn_idx]['FVC'])\n    lvalid = lgbm.Dataset(train.iloc[val_idx][features], label = train.iloc[val_idx]['FVC'])\n    \n    clf = lgbm.train(params, ltrain, \n                    num_boost_round = 10000, \n                    verbose_eval = 100, \n                    early_stopping_rounds = 400, \n                    valid_sets = [ltrain, lvalid], \n                    fobj = osic_loss.return_grad_and_hess,\n                    feval = osic_loss.return_loss\n                    )\n    oof[val_idx] = clf.predict(train.iloc[val_idx][features], num_iteration = clf.best_iteration)\n    \n    # RMSE\n    print(\"CV RMSE score: {:<8.5f}\".format(np.sqrt(mean_squared_error(train['FVC'], oof[:, 0]))))\n    # Metric\n    print(\"CV Metric: {:<8.5f}\".format(osic_loss(oof, train['FVC'])))\n    \n    fold_imp_df = pd.DataFrame()\n    fold_imp_df['feature'] = train[features].columns\n    fold_imp_df['importance'] = clf.feature_importance()\n    fold_imp_df['fold'] = n_folds + 1\n    feature_importance_df = pd.concat([feature_importance_df, fold_imp_df], axis = 0)\n    \n    predictions += clf.predict(test[features], num_iteration = clf.best_iteration) / nb_splits","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = (feature_importance_df[['feature', 'importance']]\n        .groupby('feature')\n        .mean()\n        .sort_values(by = 'importance', ascending = False).index)\nbest_features = feature_importance_df.loc[feature_importance_df.feature.isin(cols)]\n\nplt.figure(figsize = (10,14))\nsns.barplot(x = 'importance', y = 'feature', data = best_features.sort_values(by = 'importance', ascending = False))\nplt.title('LightGBM Features')\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions, oof","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[\"FVC_pred\"] = oof[:, 0]\ntrain[\"Confidence\"] = oof[:, 1]\ntest[\"FVC_pred\"] = predictions[:, 0]\ntest[\"Confidence\"] = predictions[:, 1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(input_path + 'sample_submission.csv')\nsub","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = sub.drop(columns = ['FVC', 'Confidence']).merge(test[['Patient_Week', 'FVC_pred', 'Confidence']], \n                                                           on = 'Patient_Week')\nsubmission.columns = sub.columns\nsubmission.to_csv('./submission.csv', index = False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(submission['FVC'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(submission['Confidence'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"finish = time()\nprint(strftime(\"%H:%M:%S\", gmtime(finish - start)))","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}