{"cells":[{"metadata":{},"cell_type":"markdown","source":"Main ideas\n\n* easly create subset of data to train and validate on (using  `train_mask` and `valid_mask`) \n* transform the target variable with `transform` and `inverse_transform`","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os, re\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\n\nimport xgboost\nimport lightgbm\nimport catboost\n\nfrom sklearn import ensemble, metrics, model_selection, linear_model, preprocessing, utils\n\nimport tensorflow as tf\nfrom tensorflow import keras","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_likelihood(y, p):\n    m, s = p\n    diff = np.minimum(1000, np.abs(y - m))\n    s = np.maximum(70, s)\n    lik = -np.sqrt(2) * diff / s - np.log(np.sqrt(2) * s)\n    return np.mean(lik)\n\n\ndef laplace_likelihood_bound(y, p):\n    m = p\n    diff = np.minimum(1000, np.abs(y - m))\n    s = np.maximum(70, np.sqrt(2) * diff)\n    lik = -np.sqrt(2) * diff / s - np.log(np.sqrt(2) * s)\n    return np.mean(lik)\n\n\ndef laplace_likelihood_avg(y, p):\n    m = p\n    diff = np.minimum(1000, np.abs(y - m))\n    s = np.sqrt(2) * metrics.mean_absolute_error(y, m)\n    lik = -np.sqrt(2) * diff / s - np.log(np.sqrt(2) * s)\n    return np.mean(lik)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_folds(X, y, group=None, k=5, shuffle=False, train_mask=None, valid_mask=None):\n    if isinstance(X, pd.DataFrame): X = X.values\n    if isinstance(y, pd.DataFrame): y = y.values\n    if isinstance(group, pd.DataFrame): group = group.values\n    if isinstance(train_mask, pd.DataFrame): train_mask = train_mask.values\n    if isinstance(valid_mask, pd.DataFrame): valid_mask = valid_mask.values\n    if group is None:\n        folds = list(model_selection.KFold(k, shuffle=True).split(X, y))\n    else:\n        idx = utils.shuffle(np.arange(X.shape[0]))\n        if shuffle:\n            Xs, ys, groups = X.copy()[idx], y.copy()[idx], group.copy()[idx]\n            folds = list(model_selection.GroupKFold(k).split(np.array(Xs), np.array(ys), np.array(groups)))\n        else:\n            folds = list(model_selection.GroupKFold(k).split(X, y, group))\n    if train_mask is not None:\n        for i in range(k):\n            folds[i] = (np.array([j for j in folds[i][0] if train_mask[j]]), folds[i][1])\n    if valid_mask is not None:\n        for i in range(k):\n            folds[i] = (folds[i][0], np.array([j for j in folds[i][1] if valid_mask[j]]))\n    return folds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class GroupMeanTransformer:\n    \n    def __init__(self, grouper, cols):\n        self.grouper = grouper\n        self.cols = cols\n        self.generated_features = []\n        self.means = dict()\n        pass\n    \n    def fit_transform(self, x, y=None):\n        self.fit(x, y)\n        return self.transform(x)\n        \n    def fit(self, x, y=None):\n        for gr in self.grouper:\n            group_name = '_'.join(gr) if isinstance(gr, list) else gr\n            for col in self.cols:\n                self.means[(group_name, col)] = x.groupby(gr)[col].mean()\n                self.means[(group_name, col)].name = '%s_mean_groupby_%s' % (col, group_name)\n                self.generated_features += ['%s_mean_groupby_%s' % (col, group_name)]\n                self.generated_features += ['delta_%s_mean_groupby_%s' % (col, group_name)]\n        pass\n    \n    def transform(self, x):\n        x = x.copy()\n        for k in self.means.keys():\n            col = k[1]\n            m = self.means[k]\n            x[m.name] = x[m.index.names].merge(m, on=m.index.names)[m.name]\n            x['delta_' + m.name] = x[col] - x[m.index.names].merge(m, on=m.index.names)[m.name]\n        return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def _transform(x, mode: str=None, df=None):\n    if mode is None: return x\n    if mode.upper() == \"DIFF\": return x - df.FVC_base\n    if mode.upper() == \"DIFF_FVC\": return np.where(df.n_weeks == 0, 0.0, x - df.FVC_base)\n    if mode.upper() == \"PERC\": return x / df.FVC_base\n    if mode.upper() == \"PERC_FVC\": return np.where(df.n_weeks == 0, 1.0, x / df.FVC_base)\n    if mode.upper() == \"WEEKLY\": return x / np.maximum(1.0, np.abs(df.n_weeks))\n    if mode.upper() == \"SYMLOG_WEEKLY\": return x / np.log(1 + np.maximum(1.0, np.abs(df.n_weeks)))\n    if re.match(\"^LOG\", mode.upper()): \n        gamma = re.sub(\"^LOG[:]*\", \"\", mode.upper())\n        gamma = float(gamma) if len(gamma) > 0 else 0.0\n        return np.log(gamma + x)\n    if re.match(\"^BOXCOX\", mode.upper()):\n        gamma = re.sub(\"^BOXCOX[:]*\", \"\", mode.upper())\n        gamma = float(gamma) if len(gamma) > 0 else 1.0\n        return (x ** gamma - 1.0) / gamma\n    return None\n\n\ndef _inverse_transform(x, mode=None, df=None):\n    if mode is None: return x\n    if mode.upper() == \"DIFF\": return x + df.FVC_base\n    if mode.upper() == \"DIFF_FVC\": return np.where(df.n_weeks == 0, df.FVC_base, x + df.FVC_base)\n    if mode.upper() == \"PERC\": return x * df.FVC_base\n    if mode.upper() == \"PERC_FVC\": return np.where(df.n_weeks == 0, df.FVC_base, x * df.FVC_base)\n    if mode.upper() == \"WEEKLY\": return x * np.maximum(1.0, np.abs(df.n_weeks))\n    if mode.upper() == \"SYMLOG_WEEKLY\": return x * np.log(1 + np.maximum(1.0, np.abs(df.n_weeks)))\n    if re.match(\"^LOG\", mode.upper()): \n        gamma = re.sub(\"^LOG[:]*\", \"\", mode.upper())\n        gamma = float(gamma) if len(gamma) > 0 else 0.0\n        return np.exp(x) - gamma\n    if re.match(\"^BOXCOX\", mode.upper()):\n        gamma = re.sub(\"^BOXCOX[:]*\", \"\", mode.upper())\n        gamma = float(gamma) if len(gamma) > 0 else 1.0\n        return (1.0 + x * gamma) ** (1 / gamma)\n    return None\n    \n\ndef transform(x, mode=None, df=None):\n    out = x.copy()\n    if isinstance(mode, list):\n        for m in mode:\n            out = _transform(out, mode=m, df=df)\n    if isinstance(mode, str):\n        out = _transform(out, mode=mode, df=df)\n    return out\n\n\ndef inverse_transform(x, mode=None, df=None):\n    out = x.copy()\n    if isinstance(mode, list):\n        for m in mode[::-1]:\n            out = _inverse_transform(out, mode=m, df=df)\n    if isinstance(mode, str):\n        out = _inverse_transform(out, mode=mode, df=df)\n    return out","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def feature_eng(df):\n    df = df.copy()\n    df['n_weeks'] = df['Weeks_target'] - df['Weeks_base']\n    df['symlog_n_weeks'] = np.sign(df['n_weeks']) * np.log(1 + np.abs(df['n_weeks']))\n    df['symlog_n_weeks2'] = np.sign(df['n_weeks']) * np.log(1 + np.abs(df['n_weeks']) ** 2)\n    df['expdecay_n_weeks'] = np.exp(-np.abs(df['n_weeks']))\n    df['Sex_female'] = (df['Sex'] == 'Female').astype('float')\n    df['Smoking_ex'] = (df['SmokingStatus'] == 'Ex-smoker').astype(int)\n    df['Smoking_currently'] = (df['SmokingStatus'] == 'Currently smokes').astype(int)\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_folder = \"../input/osic-pulmonary-fibrosis-progression\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(os.path.join(data_folder, \"train.csv\")).drop_duplicates(keep=False, subset=['Patient', 'Weeks'])\ndf_train['n_obs'] = df_train.groupby('Patient').Weeks.cumcount()\ndf_train['till_last'] = df_train.groupby('Patient').Weeks.transform('count') - 1 - df_train['n_obs']\ndf_train = df_train.merge(df_train.drop(['Age', 'Sex', 'Percent', 'SmokingStatus'], 1), on='Patient', suffixes=['_base', '_target'])\ncols_before_fe = df_train.columns\ndf_train = feature_eng(df_train)\nFEATURES = ['FVC_base', 'Percent', 'Age'] + [c for c in df_train.columns if c not in cols_before_fe]\ndf_train[FEATURES].head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mask = (df_train.n_weeks > 0) & df_train.n_obs_base.between(0, 5) & df_train.till_last_target.between(0, 3)\nvalid_mask = (df_train.n_weeks > 0) & df_train.n_obs_base.between(0, 0) & df_train.till_last_target.between(0, 2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test = pd.read_csv(os.path.join(data_folder, \"test.csv\")).rename(columns={'Weeks': 'Weeks_base', 'FVC': 'FVC_base'})\ndf_test = df_test.assign(k=0).merge(pd.DataFrame({'Weeks_target': np.arange(-12, 133 + 1), 'k': 0})).drop('k', 1)\ndf_test = feature_eng(df_test)\ndf_test.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(50 * \"#\")\nprint(\"##  Train:  %4d rows with %3d unique patients  ##\" % (df_train[train_mask].shape[0], df_train[train_mask].Patient.nunique()))\nprint(\"##  Valid:  %4d rows with %3d unique patients  ##\" % (df_train[valid_mask].shape[0], df_train[valid_mask].Patient.nunique()))\nprint(\"##  Test:   %4d rows with %3d unique patients  ##\" % (df_test.shape[0], df_test.Patient.nunique()))\nprint(50 * \"#\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gr = GroupMeanTransformer(['Sex', 'SmokingStatus'], ['FVC_base'])\ngr.fit(df_train)\ndf_train = gr.transform(df_train)\ndf_test = gr.transform(df_test)\nFEATURES += [f for f in gr.generated_features if re.match(\"^delta_\", f)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"MODE = [\"PERC_FVC\", \"LOG\"]\nSIGMA_MODE = [\"PERC\"]\n\nX = df_train[FEATURES].copy().values\ny = df_train['FVC_target'].copy().values\ngroup = df_train['Patient'].values\nX_test = df_test[FEATURES].copy().values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train[mean_target == mean_target.max()]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N_FOLDS = 10\nfolds = build_folds(X, y, group, k=N_FOLDS, train_mask=train_mask, valid_mask=valid_mask)\nvalid_folds = pd.DataFrame({'idx': np.concatenate([np.array(f[1]) for fi, f in enumerate(folds)]),\n                            'fold': np.concatenate([np.repeat(np.array(fi).reshape(1), len(f[1])) for fi, f in enumerate(folds)])}).set_index('idx')\n\nprep = preprocessing.Normalizer()\nZ = prep.fit_transform(X)\nZ_test = prep.transform(X_test)\n\nmean_target = transform(y, MODE, df_train)\nplt.hist(mean_target, 40); ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_oof = np.nan * np.zeros((N_FOLDS, X.shape[0]))\npred_test = np.nan * np.zeros((N_FOLDS, X_test.shape[0]))\nfor i, (idxT, idxV) in enumerate(tqdm(folds)):\n    model = linear_model.LassoCV()\n    model.fit(Z[idxT], mean_target[idxT])\n    pred_oof[i] = model.predict(Z)\n    pred_oof[i, idxT] = np.nan\n    pred_test[i] = model.predict(Z_test)\npred_mean = inverse_transform(np.nanmean(pred_oof, 0), MODE, df_train)\ntest_mean = inverse_transform(np.nanmean(pred_test, 0), MODE, df_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"opt_sigma = np.sqrt(2) * np.clip(np.abs(y - pred_mean), 0, 1000)\nsigma_target = transform(opt_sigma, SIGMA_MODE, df_train)\nplt.hist(sigma_target[train_mask | valid_mask], 50);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_oof = np.nan * np.zeros((N_FOLDS, X.shape[0]))\npred_test = np.nan * np.zeros((N_FOLDS, X_test.shape[0]))\nfor i, (idxT, idxV) in enumerate(tqdm(folds)):\n    model = linear_model.LinearRegression()\n    model.fit(Z[idxT], sigma_target[idxT])\n    pred_oof[i] = model.predict(Z)\n    pred_oof[i, idxT] = np.nan\n    pred_test[i] = model.predict(Z_test)\npred_sigma = inverse_transform(np.nanmean(pred_oof, 0), SIGMA_MODE, df_train)\ntest_sigma = inverse_transform(np.nanmean(pred_test, 0), SIGMA_MODE, df_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Monitoring","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"lll_overall = laplace_likelihood(y[valid_mask], [pred_mean[valid_mask], pred_sigma[valid_mask]])\nlllb_overall = laplace_likelihood_bound(y[valid_mask], pred_mean[valid_mask])\ny_rmse_overall = metrics.mean_squared_error(y[valid_mask], pred_mean[valid_mask]) ** 0.5\ny_mae_overall = metrics.mean_absolute_error(y[valid_mask], pred_mean[valid_mask])\ns_rmse_overall = metrics.mean_squared_error(opt_sigma[valid_mask], pred_sigma[valid_mask]) ** 0.5\ns_mae_overall = metrics.mean_absolute_error(opt_sigma[valid_mask], pred_sigma[valid_mask])\n\n\n_df_base = pd.DataFrame({'y': y, 'so': opt_sigma, 'p': pred_mean, 's': pred_sigma}).merge(valid_folds, left_index=True, right_index=True)\nlllb   = _df_base.groupby('fold').apply(lambda x: laplace_likelihood_bound(x['y'], x['p']))\nlll    = _df_base.groupby('fold').apply(lambda x: laplace_likelihood(x['y'], [x['p'], x['s']]))\ny_rmse = _df_base.groupby('fold').apply(lambda x: metrics.mean_squared_error(x['y'], x['p']) ** 0.5)\ny_mae  = _df_base.groupby('fold').apply(lambda x: metrics.mean_absolute_error(x['y'], x['p']))\ns_rmse = _df_base.groupby('fold').apply(lambda x: metrics.mean_squared_error(x['so'], x['s']) ** 0.5)\ns_mae  = _df_base.groupby('fold').apply(lambda x: metrics.mean_absolute_error(x['so'], x['s']))\n\n#\nprint(\"        OVERALL  ||  FOLD-WISE\")\nprint(\"Bound:  %7.4f  ||  %7.4f  (+- %7.4f)\" % (lllb_overall, lllb.mean(), 2 * lllb.std()))\nprint(\"Score:  %7.4f  ||  %7.4f  (+- %7.4f)\" % (lll_overall, lll.mean(), 2 * lll.std()))\nprint(\"yRMSE:  %7.2f  ||  %7.2f  (+- %7.2f)\" % (y_rmse_overall, y_rmse.mean(), 2 * y_rmse.std()))\nprint(\"y-MAE:  %7.2f  ||  %7.2f  (+- %7.2f)\" % (y_mae_overall, y_mae.mean(), 2 * y_mae.std()))\nprint(\"sRMSE:  %7.2f  ||  %7.2f  (+- %7.2f)\" % (s_rmse_overall, s_rmse.mean(), 2 * s_rmse.std()))\nprint(\"s-MAE:  %7.2f  ||  %7.2f  (+- %7.2f)\" % (s_mae_overall, s_mae.mean(), 2 * s_mae.std()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame({'LLL': lll, 'Bound': lllb}).boxplot(figsize=(12, 4));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 6))\nplt.subplot(1, 2, 1)\nplt.scatter(y[valid_mask], pred_mean[valid_mask], alpha=0.2);\nplt.xlim(min(plt.xlim()[0], plt.ylim()[0]), max(plt.xlim()[1], plt.ylim()[1]))\nplt.ylim(min(plt.xlim()[0], plt.ylim()[0]), max(plt.xlim()[1], plt.ylim()[1]))\nplt.plot(plt.xlim(), plt.ylim(), color='tab:red', alpha=0.5, linestyle='--')\nplt.subplot(1, 2, 2)\nplt.scatter(opt_sigma[valid_mask], pred_sigma[valid_mask], alpha=0.2)\nplt.xlim(min(plt.xlim()[0], plt.ylim()[0]), max(plt.xlim()[1], plt.ylim()[1]))\nplt.ylim(min(plt.xlim()[0], plt.ylim()[0]), max(plt.xlim()[1], plt.ylim()[1]))\nplt.plot(plt.xlim(), plt.ylim(), color='tab:red', alpha=0.5, linestyle='--');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 6))\nplt.subplot(1, 2, 1)\nplt.scatter(df_train.Weeks_target[valid_mask], pred_mean[valid_mask], alpha=0.5)\nplt.subplot(1, 2, 2)\nplt.scatter(df_train.Weeks_target[valid_mask], pred_sigma[valid_mask], alpha=0.5);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 8))\nfor i, pid in enumerate(df_train[valid_mask].Patient.unique()[:12]):\n    plt.subplot(3, 4, i + 1)\n    idx = (df_train.Patient == pid) & (df_train.n_obs_base == 0)\n    plt.fill_between(\n        df_train[idx].Weeks_target,\n        pred_mean[idx] - 1 * pred_sigma[idx],\n        pred_mean[idx] + 1 * pred_sigma[idx],\n        color='tab:blue', alpha=0.1\n    )\n    plt.fill_between(\n        df_train[idx].Weeks_target,\n        pred_mean[idx] - 2 * pred_sigma[idx],\n        pred_mean[idx] + 2 * pred_sigma[idx],\n        color='tab:blue', alpha=0.1\n    )\n    plt.plot(df_train[idx].Weeks_target, pred_mean[idx], marker='x', color='tab:blue')\n    plt.plot(df_train[idx].Weeks_target, df_train[idx].FVC_target, marker='o', color='tab:red')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i, pid in enumerate(df_train[valid_mask].Patient.unique()):\n    idx = (df_train.Patient == pid) & (df_train.n_obs_base == 0)\n    plt.plot(df_train[idx].n_weeks, pred_mean[idx] / df_train[idx].FVC_base, color='tab:blue', alpha=0.2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = df_test.copy()[['Patient', 'Weeks_target']]\nsubmission['Patient_Week'] = submission['Patient'] + \"_\" + submission['Weeks_target'].astype('str')\nsubmission['FVC'] = test_mean\nsubmission['Confidence'] = test_sigma\nsubmission = submission.sort_values(['Weeks_target', 'Patient'])[['Patient_Week', 'FVC', 'Confidence']]\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 8))\nfor i, pid in enumerate(df_test[valid_mask].Patient.unique()):\n    plt.subplot(2, 3, i + 1)\n    idx = (df_test.Patient == pid)\n    plt.fill_between(\n        df_test[idx].Weeks_target,\n        test_mean[idx] - 1 * test_sigma[idx],\n        test_mean[idx] + 1 * test_sigma[idx],\n        color='tab:blue', alpha=0.1\n    )\n    plt.fill_between(\n        df_test[idx].Weeks_target,\n        test_mean[idx] - 2 * test_sigma[idx],\n        test_mean[idx] + 2 * test_sigma[idx],\n        color='tab:blue', alpha=0.1\n    )\n    plt.plot(df_test[idx].Weeks_target, test_mean[idx], color='tab:blue')\n    plt.plot(df_test[idx].Weeks_base, df_test[idx].FVC_base, marker='o', color='tab:red')","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}