{"cells":[{"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 matplotlib.pyplot as plt\nimport tensorflow as tf\nimport glob\n\nfrom tqdm import tqdm\n\nfrom sklearn.model_selection import (\n    cross_val_score, GroupKFold, GridSearchCV, \n    cross_val_predict, RandomizedSearchCV)\n\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import make_scorer\nfrom sklearn.feature_selection import SelectKBest, f_regression\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom sklearn.base import RegressorMixin, BaseEstimator\nfrom sklearn.preprocessing import PolynomialFeatures, OneHotEncoder, OrdinalEncoder\n\nplt.style.use(\"dark_background\")\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# kaggle\nmain_dir = \"../input/osic-pulmonary-fibrosis-progression\"\n\n!ls {main_dir}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files = tf.io.gfile.glob(main_dir+\"/train/*/*\")\ntest_files = tf.io.gfile.glob(main_dir+\"/test/*/*\")\nsample_sub = pd.read_csv(main_dir + \"/sample_submission.csv\")\ntrain = pd.read_csv(main_dir + \"/train.csv\")\ntest = pd.read_csv(main_dir + \"/test.csv\")\n\nprint (\"Number of train patients: {}\\nNumber of test patients: {:4}\"\n       .format(train.Patient.nunique(), test.Patient.nunique()))\n\nprint (\"\\nTotal number of Train patient records: {}\\nTotal number of Test patient records: {:6}\"\n       .format(len(train_files), len(test_files)))\n\ntrain.shape, test.shape, sample_sub.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_log_likelihood(y_true, y_pred, sigma=70):\n    # values smaller than 70 are clipped\n    sigma_clipped = tf.maximum(sigma, 70)\n\n    # errors greater than 1000 are clipped\n    delta_clipped = tf.minimum(tf.abs(y_true - y_pred), 1000)\n    \n    # type cast them suitably\n    delta_clipped = tf.cast(delta_clipped, dtype=tf.float32)\n    sigma_clipped = tf.cast(sigma_clipped, dtype=tf.float32)\n    \n    # score function\n    score = - tf.sqrt(2.0) * delta_clipped / sigma_clipped - tf.math.log(tf.sqrt(2.0) * sigma_clipped)\n    \n    return tf.reduce_mean(score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create a simple scorer function that can alter the value of sigma\ndef l1(s):\n    def scorer_func(x, y, sigma=s):\n        return laplace_log_likelihood(x, y, sigma=s).numpy()\n    \n    return make_scorer(scorer_func, greater_is_better=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def base_shift(data, q=50):\n    x = data.copy()\n\n    # Create base_Week, Base_FVC and Base_Percent for train\n    temp = (x.groupby(\"Patient\")\n            .apply(lambda x: x.loc[int(\n                np.percentile(x['Weeks'].index, q=q)\n            ), [\"Weeks\", \"FVC\", \"Percent\"]]))\n\n    temp.rename(\n        {\"Weeks\": \"Base_Week\", \n         \"FVC\": \"Base_FVC\", \n         \"Percent\": \"Base_Percent\"}, \n        axis=1, inplace=True)\n\n    # merge it with train data\n    x = x.merge(temp, on='Patient')\n\n    # create week offsets\n    x['Week_Offset'] = x['Weeks'] - x['Base_Week']\n    \n    return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def multi_baseweek_frame(data, display=True):\n    '''Function to return multiple base week frames -> instead of creating one base week,\n    creates several to help our model learn better and predict past and future data better.'''       \n    \n    op = data.merge(\n        data[['Patient', 'Weeks', 'FVC', 'Percent']].rename(\n            {\"Weeks\": \"Base_Week\", \n             \"FVC\": \"Base_FVC\", \n             \"Percent\": \"Base_Percent\"}, axis=1), \n        on='Patient')\n\n    # create week offsets\n    op['Week_Offset'] = op['Weeks'] - op['Base_Week']\n\n    # only take those rows with offset other than 0\n    op = op[op['Week_Offset'] != 0]\n\n    if display:\n        # number of training samples\n        print (\"Number of Samples:{:5} -> {:5}\\nNumber of Columns:{:5} -> {:5}\".format(\n            data.shape[0], op.shape[0], data.shape[1], op.shape[1]))\n    \n    return op.sort_values(by=['Patient', 'Base_Week']).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model_data(data, cat_cols, num_cols, to_drop, cat_method='1h', transform_stats=None,\n                   age_bins=None, train=True, display_stats=True, math=None, factor=False):\n    \n    '''Our pipeline for this notebook. This portion is complex, could be written much more efficiently \n    using simple sklearn tools. I have this bad habit of reinventing the wheel from scratch.'''\n    \n    X = data.copy().reset_index(drop=True)\n    \n    ##########################################################################\n    \n    if age_bins:    \n        X['binned_age'] = pd.cut(X['Age'], bins=range(0, 101, 100//(age_bins-1))).cat.codes / age_bins\n        to_drop = to_drop + ['Age']   \n        \n    if math:\n        prod = np.ones(X.shape[0])\n        if 'sin' in math:\n            X['Sin_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.sin)\n            prod = prod * X['Sin_week']\n        if 'cos' in math:\n            X['Cos_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.cos)\n            prod = prod * X['Cos_week']\n        if 'tan' in math:\n            X['Tan_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.tan)\n            prod = prod * X['Cos_week']\n            \n        if len(math) > 1:\n            X['Math_Prod'] = prod\n            \n    if factor:\n        X['factor'] = X['Base_FVC'] / X['Base_Percent']\n        \n    ##########################################################################\n    \n    if cat_cols != []:\n    \n        if cat_method == 'ord': # ordinal encoding for tree based models\n            global ordenc\n            if train:\n                from sklearn.preprocessing import OrdinalEncoder\n                ordenc = OrdinalEncoder()\n                X = X.merge(\n                    pd.DataFrame(\n                        ordenc.fit_transform(X[cat_cols]).astype(int),\n                        columns=map(lambda x: x+\"_ord\", cat_cols)),\n                    left_index=True, right_index=True)\n\n            else:\n                X = X.merge(\n                    pd.DataFrame(\n                        ordenc.transform(X[cat_cols]).astype(int),\n                        columns=map(lambda x: x+\"_ord\", cat_cols)),\n                    left_index=True, right_index=True)           \n\n        elif cat_method == '1h': # one hot encoding\n            global onehenc\n            if train:\n                onehenc = OneHotEncoder()\n                X = X.merge(\n                    pd.DataFrame(\n                        onehenc.fit_transform(X[cat_cols]).todense(),\n                        columns=[*np.concatenate(onehenc.categories_)]),\n                    left_index=True, right_index=True)\n\n            else:\n                X = X.merge(\n                    pd.DataFrame(\n                        onehenc.transform(X[cat_cols]).todense(),\n                        columns=[*np.concatenate(onehenc.categories_)]),\n                    left_index=True, right_index=True)\n\n        elif cat_method == 'poly': # polynomial feature encoding\n            global cat_comb\n            if train:\n                cat_comb = np.array(\n                    np.meshgrid(*[X[cat].unique() for cat in cat_cols])\n                ).T.reshape(-1, len(cat_cols))\n\n            for combination in cat_comb:\n                name = \"_\".join(map(str, combination))\n                X[name] = 1\n                for i in range(len(cat_cols)):\n                    X[name] = X[name] & (X[cat_cols[i]] == combination[i]).astype(int)\n\n    # drop the columns after they have been encoded\n    to_drop = to_drop + cat_cols\n                \n    ##########################################################################\n    \n    global stats\n    if train:\n        if transform_stats is None:\n            # saving stats for the future\n            stats = X.describe().T\n        else:\n            stats = transform_stats\n\n    # lets scale the numeric columns (We scale it with max possible values)\n    for col in num_cols:\n        \n        if (not train) and (col not in X.columns):\n            continue\n            \n        X[col] = (X[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n        \n    ##########################################################################\n\n    global x_cols\n    if train:\n        \n        Y = X['FVC'].dropna()\n        \n        if display_stats:\n        \n            # print out how well our features would do\n            print (X.corr()['FVC'].abs().sort_values(ascending=False)[1:])\n        \n        X = X.drop(to_drop, axis=1)\n        x_cols = X.columns\n        \n        return X, Y\n    \n    else:\n        \n        X = X.drop(to_drop, axis=1, errors='ignore')\n        X = X[x_cols]\n        \n        return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def augment_train_cosine(data, n_similar=3, threshold=0.25, display_sample=True):\n    \n    '''\n    - `n_similar` is number of patients we cluster at a time more the cluster, more eratic it gets.\n    - `threshold` is used for as a measure to counter the influence of \"outliers\"\n    '''\n    \n    from sklearn.metrics.pairwise import cosine_similarity\n\n    temp = base_shift(data.copy(), q=0)\n\n    # feature engineering -> simiarity of percentage changes across patients\n    # we use slopes, sex, SmokingStatus and the way the Percent features vary \n    # from Base_Percent as preditive features for clustering patients together.\n    \n    temp['present_minus_past'] = temp.groupby(\"Patient\")['FVC'].transform('diff').fillna(0)\n    temp['Week_diff'] = temp.groupby(\"Patient\")['Weeks'].transform('diff').fillna(0)\n    temp['pms'] = (temp['present_minus_past'] / temp['Week_diff']).replace([np.inf, -np.inf]).fillna(0)\n    temp['pms_min']  = temp.groupby(\"Patient\")['pms'].transform('min')\n    temp['pms_25']   = temp.groupby(\"Patient\")['pms'].transform(lambda x: np.percentile(x, q=25))\n    temp['pms_mean'] = temp.groupby(\"Patient\")['pms'].transform('mean')\n    temp['pms_75']   = temp.groupby(\"Patient\")['pms'].transform(lambda x: np.percentile(x, q=75))\n    temp['pms_max']  = temp.groupby(\"Patient\")['pms'].transform('max')\n    temp['pms_sum']  = temp.groupby(\"Patient\")['pms'].transform('sum')\n    \n    temp['pmb_avg'] = (temp['Percent'] - temp['Base_Percent']).groupby(temp.Patient).transform(\"mean\")\n    temp['p_std']    = temp.groupby(\"Patient\")['Percent'].transform('std')\n\n    temp = temp.merge(\n        pd.concat([\n            temp.groupby(\"Patient\").apply(\n                lambda x: (x['Percent'].values[-1] - x['Percent'].values[0]) / \n                (x['Weeks'].values[-1] - x['Weeks'].values[0])).rename(\"Slope\"),\n            \n            temp.groupby(\"Patient\").apply(\n                lambda x: x['Week_Offset'].iloc[np.argmax(x['pms'])]).rename(\"pmsw_max\"),    \n            temp.groupby(\"Patient\").apply(\n                lambda x: x['Week_Offset'].iloc[np.argmin(x['pms'])]).rename(\"pmsw_min\")\n            \n        ], axis=1, ignore_index=False), on='Patient')\n\n    # we take just the head row for each patient for comparison\n    temp = temp.groupby('Patient').head(1).reset_index(drop=True)\n\n    # another measure of similartiy\n    temp['factor'] = temp['Base_FVC'] / temp['Base_Percent']\n\n    # features to use for similarity clustering\n    train_cols = [\n        # compulsary features\n        'Patient', #'Sex', 'SmokingStatus'\n        \n        # optionally added features for clustering\n        # 'pms_min', 'pms_max', 'p_std',\n        'pms_25', 'pms_mean', 'pms_75', \n        'pms_sum', 'pmsw_min', 'pmsw_max', \n        'pmb_avg', 'Slope', 'factor'\n    ]\n    \n    cat_cols = np.intersect1d(['Sex', 'SmokingStatus'], train_cols)\n\n    temp = temp[train_cols]\n    temp = pd.get_dummies(temp, columns=cat_cols, drop_first=True, prefix='', prefix_sep='')\n\n    # including that patient, find n_similar more patients\n    n_similar += 1\n    groups = (pd.DataFrame(\n        np.argsort(\n            # cosine similarity to get their similarity scores\n            cosine_similarity(temp.drop(\"Patient\", 1), temp.drop('Patient', 1)))\n        [:, -1:-n_similar-1:-1]))\n\n    # cosine similarity is symmetric, so we remove the redundant ones\n    groups = groups[~pd.DataFrame(np.sort(groups.values, axis=1)).duplicated()]\n    \n    # convert the indices to patient ids\n    groups = groups.applymap(lambda x: temp.Patient.to_dict()[x]).apply(list, axis=1).to_dict()\n\n    # the bottle neck of this function\n    aug_data = []\n    for group in tqdm(groups.values(), disable=not display_sample):\n\n        temp = base_shift(train[train.Patient.isin(group)], q=0)\n\n        temp['base_per_diff_from_mean'] = temp['Base_Percent'] - temp['Base_Percent'].unique().mean()\n        temp['Percent_shifted'] = temp['Percent'] - temp['base_per_diff_from_mean']\n\n        temp['base_week_diff_from_mean'] = temp['Base_Week'] - temp['Base_Week'].unique().mean()\n        temp['Week_shifted'] = temp['Weeks'] - temp['base_week_diff_from_mean']\n\n        temp = pd.merge_ordered(\n            temp.drop(['Age', 'Sex', 'SmokingStatus', 'Base_Week', 'Week_Offset'], axis=1),\n\n            # we obtain the mean only for those with samples greater than threshold %\n            # these mean values are fit as such as the expected Percent\n            (temp.groupby(\"Week_shifted\")['Percent_shifted']\n             .agg(['mean', 'count']).query(f'count > {n_similar * threshold}')\n             .drop(\"count\", 1)),\n\n            on='Week_shifted', left_by='Patient'\n        )\n\n        aug_data.append(temp)\n\n    temp = pd.concat(aug_data).reset_index(drop=True)\n\n    # recuring features can simply be padded\n    temp[['base_per_diff_from_mean', 'base_week_diff_from_mean', 'Base_Percent', 'Base_FVC']] = (\n        temp.groupby(\"Patient\")[['base_per_diff_from_mean', 'base_week_diff_from_mean', \n                                 'Base_Percent', 'Base_FVC']].fillna(method='ffill'))\n\n    # get back the weeks from shifted weeks\n    temp['Week_aug'] = temp['Week_shifted'] + temp['base_week_diff_from_mean']\n\n    # For those percent values already present copy them, augment the rest\n    temp['Percent_aug'] = np.where(\n        temp['Percent'].isna(), \n        temp['mean'] + temp['base_per_diff_from_mean'],\n        temp['Percent'])\n\n    # drop/clean those which have neither \n    temp = temp[~temp['Percent_aug'].isna()]\n\n    # fill the FVC values using percent -> FVC corelation\n    test_ids = temp['FVC'].isna()\n    temp.loc[test_ids, 'FVC'] = LinearRegression().fit(\n        temp.loc[~test_ids, ['Base_Percent', 'Percent_aug', 'Base_FVC']], \n        temp.loc[~test_ids, 'FVC']).predict(\n        temp.loc[test_ids, ['Base_Percent', 'Percent_aug', 'Base_FVC']])\n\n    temp = temp.groupby([\"Patient\", 'Week_aug']).mean().reset_index()\n\n    # retain only those columns that we need & rename them to match train\n    temp = temp[['Patient', 'Week_aug', 'Percent_aug', 'FVC']]\n    temp = temp.rename({'Week_aug': 'Weeks', 'Percent_aug': 'Percent'}, axis=1)\n    \n    # some weeks may end up slightly shifted to left or right\n    temp['Weeks'] = temp['Weeks'].astype(int)\n    \n    # add these columns to be able to fit inside multibaseweek pipeline\n    temp = temp.merge(\n        (data[['Patient', 'Sex', 'Age', 'SmokingStatus']]\n         .groupby('Patient').head(1).reset_index(drop=True)), \n        \n        on='Patient')\n    \n    if display_sample:\n        \n        print (\"Data augmented by factor: {:.2f}x\".format(1 + (\n            temp.shape[0] - data.shape[0]) / data.shape[0]))\n\n        f, ax = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))\n        for i, pat in enumerate(temp.Patient.unique()[:4]):\n\n            ax[i][0].plot(*list(zip(*data[data.Patient == pat][['FVC', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n\n            ax[i][1].plot(*list(zip(*data[data.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n\n            ax[i][0].scatter(*list(zip(*temp[temp.Patient == pat][['FVC', 'Weeks']].values))\n                          [::-1], c='r')\n\n            ax[i][1].scatter(*list(zip(*temp[temp.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='r')\n            \n\n            ax[i][0].set(xlabel='Weeks', ylabel='FVC')\n            ax[i][1].set(xlabel='Weeks', ylabel='Percent')\n            f.suptitle(\"FVC & Percent Augmentation\", y=.9)\n    \n    return temp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def augment_train_naive(\n    data, steps=5, method='index', noise=.25, val_split=0.25, \n    end_pts=[None, None], display_sample=True):\n    \n    '''\n    end_pts -> start and end of augmentation, if None defaults to min/max for that patient\n    '''\n    \n    temp = data[['Patient','Weeks', 'FVC', 'Percent']].merge(\n        \n        ((data.groupby(\"Patient\")['Weeks']\n         .apply(lambda x: pd.Series(\n             np.union1d(np.arange(\n                 end_pts[0] if end_pts[0] else x.min(), \n                 end_pts[1] if end_pts[1] else x.max(), \n                 step=steps), x))\n               ).reset_index(level=0))),\n\n        on=['Patient', 'Weeks'], how='right')\n\n    temp.loc[:, ['FVC', 'Percent']] = (\n        temp.groupby(\"Patient\")[['FVC', 'Percent']]\n        .apply(lambda x: (\n            # interpolate \n            x.interpolate(method=method, limit_direction='both') + \n            \n            # noise factor: Gaussian noice + standard deviation of resp features\n            # (we assume std of percent is scaled up version of std of FVC)\n            (x.std().values * np.random.uniform(-noise, noise, [len(x), 1])))))\n\n    temp = temp.merge(\n        data.groupby(\"Patient\")[['Patient', 'Age', 'Sex', 'SmokingStatus']].head(1),\n        on='Patient')\n    \n    if display_sample:\n        f, ax = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))\n\n        for i, pat in enumerate(temp.Patient.unique()[:4]):\n\n            ax[i][0].plot(*list(zip(*data[data.Patient == pat][['FVC', 'Weeks']].values))\n                [::-1], c='g', alpha=0.7)\n            \n            ax[i][1].plot(*list(zip(*data[data.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n            \n            ax[i][0].scatter(*list(zip(*temp[temp.Patient == pat][['FVC', 'Weeks']].values))[::-1], c='r')\n            \n            ax[i][1].scatter(*list(zip(*temp[temp.Patient == pat][['Percent', 'Weeks']].values))[::-1], c='r')\n            \n            ax[i][0].set(xlabel='Weeks', ylabel='FVC')\n            ax[i][1].set(xlabel='Weeks', ylabel='Percent')\n            f.suptitle(\"FVC & Percent Augmentation\", y=.9)\n            \n        print (\"Data augmented by factor: {:.2f}x\".format(1 + (\n            temp.shape[0] - data.shape[0]) / data.shape[0]))\n    \n    return temp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = sample_sub.Patient_Week.str.extract(\"(ID\\w+)_(\\-?\\d+)\").rename({0: \"Patient\", 1: \"Weeks\"}, axis=1)\nsub['Weeks'] = sub['Weeks'].astype(int)\nsub = pd.merge(sub, test[['Patient', 'Sex', 'SmokingStatus']], on='Patient')\nsub[\"Patient_Week\"] = sub.Patient + \"_\" + sub.Weeks.astype(str)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class GBR(RegressorMixin, BaseEstimator):\n    def __init__(self, alpha=.75, **params):\n        self.alpha = alpha\n        self.umodel = self._create_model(loss='quantile', q=self.alpha, **params)\n        self.mmodel = self._create_model(loss='lad', **params)\n        self.lmodel = self._create_model(loss='quantile', q=1-self.alpha, **params)\n              \n    def _create_model(self, loss, q=.75, **params):\n        model = GradientBoostingRegressor(\n            init=LinearRegression(),\n            criterion='friedman_mse',\n            n_estimators=50, max_depth=2, \n            loss=loss, alpha=q, **params)\n        \n        return model\n        \n    def fit(self, x, y):\n        \n        self.umodel.fit(x, y)\n        self.mmodel.fit(x, y)\n        self.lmodel.fit(x, y)\n        return self\n    \n    def predict(self, X):\n        \n        return self.mmodel.predict(X)\n    \n    def predict_forecast(self, X, return_bounds=False):\n        \n        preds = self.mmodel.predict(X)\n        upper = self.umodel.predict(X)\n        lower = self.lmodel.predict(X)\n        \n        if return_bounds:\n            return preds, upper, lower\n        else:\n            return preds, (upper - lower)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# saving stats for feeding to pipe\nstats = multi_baseweek_frame(pd.concat([train, test]), display=False).describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"op = multi_baseweek_frame(\n    augment_train_cosine(train, display_sample=False, n_similar=3, threshold=.25)\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"alpha = .85","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# preprocessing pipe essentials\ncat_cols = ['Sex', 'SmokingStatus']\nto_drop = ['FVC', 'Percent', 'Weeks', 'factor', 'Base_Percent']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\ncat_method = 'ord'\nmath = []\nage_bins = 5\n\n# cross validation \nfolds = 7\ntotal_patients = train.Patient.unique()\nnp.random.shuffle(total_patients)\nval_len = len(total_patients) // folds\n\n# data frame to hold the predictions on train & test data\ntemp = pd.DataFrame()\npreds = pd.DataFrame()\n\n# creating data suitable for model fitting and predictions\nX, Y = get_model_data(\n    op, factor=True, num_cols=num_cols, cat_cols=cat_cols, \n    age_bins=age_bins, to_drop=to_drop, display_stats=False, \n    cat_method=cat_method, transform_stats=stats, math=math)\n\nX_VAL = base_shift(train, q=0)\nY_VAL = X_VAL['FVC'].dropna()\n\n# percent value as the base percent\nX_VAL['Percent'] = X_VAL['Base_Percent']\n\nX_VAL = get_model_data(\n    X_VAL, num_cols=num_cols, cat_cols=cat_cols, to_drop=to_drop, factor=True,\n    cat_method=cat_method, train=False, age_bins=age_bins, math=math)\n\n# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n\n# percent value as the base percent\nx_test['Percent'] = x_test['Base_Percent']\n\nx_test = get_model_data(\n    x_test, cat_cols=cat_cols, num_cols=num_cols, to_drop=to_drop, math=math,\n    factor=True, train=False, cat_method=cat_method, \n    age_bins=age_bins).drop(\"Patient\", 1)\n\nfor i in range(folds):\n   \n    val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n    train_patients = np.setdiff1d(total_patients, val_patients)\n    \n    assert len(np.intersect1d(val_patients, train_patients)) == 0\n    \n    x, y, = (X[X.Patient.isin(train_patients)].drop(\"Patient\", 1), Y[X.Patient.isin(train_patients)])\n    x_val, y_val = (X_VAL[X_VAL.Patient.isin(val_patients)].drop(\"Patient\", 1), \n                    Y_VAL[X_VAL.Patient.isin(val_patients)])\n        \n    # creating base model, no parameter tweakingparam_distributions\n    model = GBR(alpha=alpha)\n\n    model.fit(x, y)\n    y_middle, y_upper, y_lower = model.predict_forecast(x_val, return_bounds=True)\n    y_middle_pred, y_test_conf = model.predict_forecast(x_test)\n    \n    print (\"For Fold #{} Val Score: {:.2f} @ 70 Confidence | {:.2f} @ Pred Confidence\".format(\n        i+1, - laplace_log_likelihood(y_middle, y_val), \n        - laplace_log_likelihood(y_middle, y_val, y_upper - y_lower)))\n    \n    temp = temp.append(pd.DataFrame(\n        data=np.stack([y_upper, y_lower, y_middle, y_val], axis=1),\n        columns=['upper', 'lower', 'pred', 'actual']\n    ))\n    \n    preds = preds.append(pd.DataFrame(\n        data=np.stack([y_middle_pred, y_test_conf], axis=1) / folds,\n        columns=['pred', 'Confidence']\n    ))\n    \npreds = preds.groupby(preds.index).sum()\ntemp['Confidence'] = temp['upper'] - temp['lower']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pat_scores = (\n    temp.reset_index(drop=True)\n    .groupby(train.Patient)\n     .apply(lambda x: -laplace_log_likelihood(x['actual'], x['pred'], x['Confidence']).numpy())\n).rename('scores').reset_index().sort_values(\"scores\", ascending=False)\n\npat_scores = pat_scores.head(8)\nprint (\"Worst Patient-Mean-Score: {:.2f}\".format(pat_scores.scores.mean()))\npat_scores = pat_scores.Patient.values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (\"\\n|================== Summary ==================|\\n\\\nScore on Total Dataset: {:.3f} @   70 Confidence\\n\\\nScore on Total Dataset: {:.3f} @  225 Confidence\\n\\\nScore on Total Dataset: {:.3f} @ Pred Confidence\".format(\n    -laplace_log_likelihood(temp['actual'], temp['pred'], 70),\n    -laplace_log_likelihood(temp['actual'], temp['pred'], 225),\n    -laplace_log_likelihood(temp['actual'], temp['pred'], temp['Confidence'])\n))\n\nf, ax = plt.subplots(figsize=(40, 40), nrows=4, ncols=2)\nax = ax.ravel()\nfor i, pat in enumerate(pat_scores):\n    (temp.reset_index(drop=True).loc[train.Patient == pat]\n     .drop([\"Confidence\"], 1).plot(ax=ax[i], legend=False))\nf.suptitle(\"Model's Worst Predictions\", size=30)\nf.tight_layout(rect=[0, 0.03, 1, 0.95]);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC'] = preds['pred']\nsub['Confidence'] = preds['Confidence']\n\n# final touches before submission\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"quant_submission.csv\", index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_cols = ['Sex', 'SmokingStatus']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\nto_drop = [\"FVC\", 'Percent', 'Weeks', 'Base_Week', 'Age', 'Base_Percent', 'factor']\nmath = []\n\nmulti_method = {}\nmulti_data = {}\nfolds = 7\n\nmethods = ['cubic', 'quadratic', 'cubicspline', \n           'pchip', 'akima', 'nearest', 'zero', \n           'slinear', 'linear', 'SIMILARITY_AUG']\n\n# train patients for CV\ntotal_patients = train.Patient.unique()\nval_len = len(total_patients) // folds\n\nfor method in methods:\n    \n    # create the agumented dataset\n    if method == 'SIMILARITY_AUG':\n        temp = augment_train_cosine(train, n_similar=3, threshold=0.25, display_sample=False)\n    else:\n        temp = augment_train_naive(data=train, method=method, steps=15, noise=0., display_sample=False)\n        \n    # create multi base week data\n    temp = multi_baseweek_frame(temp, display=False)\n    # save the augmented dataframe\n    multi_data[method] = temp\n    \n    # creating the data (processed for fitting)\n    X, Y = get_model_data(temp, num_cols=num_cols, cat_cols=cat_cols, to_drop=to_drop, \n        display_stats=False, cat_method='1h', math=math, factor=True)\n        \n    X_VAL = base_shift(train, q=0)\n    Y_VAL = X_VAL['FVC'].dropna()\n    X_VAL['Percent'] = X_VAL['Base_Percent']\n    X_VAL = get_model_data(\n        X_VAL, num_cols=num_cols, cat_cols=cat_cols, factor=True,\n        to_drop=to_drop, train=False, cat_method='1h', math=math)\n    \n    # shuffle the patients for cv\n    np.random.shuffle(total_patients)\n    scores = {}\n    \n    for i in range(folds):\n   \n        val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n        train_patients = np.setdiff1d(total_patients, val_patients)\n\n        x, y = X[X.Patient.isin(train_patients)], Y[X.Patient.isin(train_patients)]\n        x_val, y_val = X_VAL[X_VAL.Patient.isin(val_patients)], Y_VAL[X_VAL.Patient.isin(val_patients)]\n        \n        assert len(np.intersect1d(x_val.Patient.unique(), x.Patient.unique())) == 0\n        \n       # how does it perform on train\n        scores['Train'] = scores.get('Train', []) + [cross_val_score(\n            GBR(alpha=alpha), x.drop(\"Patient\", 1), y, \n            scoring=l1(70), cv=GroupKFold(5), groups=x.Patient).mean()]\n\n        # fit to measure model's performance\n        lr = GBR(alpha=alpha).fit(x.drop(\"Patient\", 1), y) \n        temp = lr.predict_forecast(x_val.drop(\"Patient\", 1))\n        temp = pd.DataFrame(np.stack(temp, 1), columns=['pred', 'conf'])\n        temp['actual'] = y_val.reset_index(drop=True)\n        \n        # performance on validation data\n        scores['Val'] = scores.get('Val', []) + [-laplace_log_likelihood(\n            temp['actual'], temp['pred'], 70\n        ).numpy()]\n        \n        # performance with confidence        \n        scores['ValC'] = scores.get('ValC', []) + [-laplace_log_likelihood(\n            temp['actual'], temp['pred'], temp['conf']\n        ).numpy()]\n        \n        # Worst Performing mean scores\n        scores['ValW'] = scores.get('ValW', []) + [temp.apply(lambda x: -laplace_log_likelihood(\n            x['actual'], x['pred'], x['conf']).numpy(), 1).nlargest(25).mean()]\n        \n    print (\"Method: {}\\nTrain score: {:5.2f} @ {:.2f} Variance \\\n    \\nVal Score: {:7.2f} @ {:.2f} Variance\\\n    \\nValC Score: {:6.2f} @ {:.2f} Variance\\\n    \\nWorst Score: {:5.2f} @ {:.2f} Variance\\n{}\\n\".format(\n        method.upper(), np.mean(scores['Train']), np.std(scores['Train']), \n        np.mean(scores['Val']), np.std(scores['Val']), np.mean(scores['ValC']), \n        np.std(scores['ValC']), np.mean(scores['ValW']), np.std(scores['ValW']), \"=\" * 35\n    ))\n    \n    multi_method[method] = scores","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"val_scores = [np.mean(multi_method[i]['ValW']) for i in methods]\ncut_off = np.percentile(val_scores, 50)\n        \nprint(\"At {:.1f} cutoff, the Worst Score would be {:.3f} @ Pred Conf\".format(\n    cut_off, np.mean(list(filter(lambda x: x < cut_off, val_scores)))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n\n# percent for test\nx_test['Percent'] = x_test['Base_Percent']\n\n# predictions dataframe\npreds = {i: None for i in methods}\n\n# train each model on all the saved augmented frames\nfor method in methods:\n    \n    if np.mean(multi_method[method]['ValW']) >= cut_off:\n        preds.pop(method)\n        continue\n    \n    x, y = get_model_data(\n        multi_data[method], num_cols=num_cols, cat_cols=cat_cols, \n        to_drop=to_drop, train=True, cat_method='1h', \n        display_stats=False, factor=True,\n    )\n\n    lr = GBR(alpha=alpha).fit(x.drop(\"Patient\", 1), y) \n\n    preds[method] = lr.predict_forecast(get_model_data(\n        x_test, cat_cols=cat_cols, num_cols=num_cols, to_drop=to_drop, train=False, factor=True,\n    ).drop(\"Patient\", 1), return_bounds=True)\n\npreds.keys()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = {}\nfor i in multi_method:\n    if i in preds.keys():\n        temp['Val'] = temp.get('Val', []) + [np.mean(multi_method[i]['Val'])]\n        temp['ValC'] = temp.get('ValC', []) + [np.mean(multi_method[i]['ValC'])]\n        \nprint (\"At same cutoff:\\n\\nBest val score: {:6.3f} @   70 Conf\\n\\\nBest Val score: {:6.3f} @ Pred Conf\".format(np.mean(temp['Val']), np.mean(temp['ValC'])))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(preds)\n\nmean = np.mean(np.stack(temp.iloc[0].values), 0)\nupper = np.mean(np.stack(temp.iloc[1].values), axis=0)\nlower = np.mean(np.stack(temp.iloc[2].values), axis=0)\n\nsub['FVC'] = mean\nsub['Confidence'] = upper - lower\n\n# save to csv file\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"submission.csv\", index=False)\n\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# max, min\nsub.Confidence.hist()","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}