{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction"},{"metadata":{},"cell_type":"markdown","source":"\"Confidence\" is tuned to increase by weeks\n\nThis is one of the approaches that earned my team a silver medal."},{"metadata":{},"cell_type":"markdown","source":"# Acknowledgements\nThis notebook is forked from https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter by @ulrich07 \n\nif you are kind enough to upvote my notebook, please also upvote Ulrich's."},{"metadata":{},"cell_type":"markdown","source":"# Import libraries "},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport random\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, StratifiedKFold, GroupKFold\nfrom tqdm.notebook import tqdm\n\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\n\npd.set_option('display.max_columns', 60)\npd.set_option('display.max_rows', 100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    \nseed_everything(42)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing"},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\n\ntr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nchunk = pd.read_csv(f\"{ROOT}/test.csv\")\n\nprint(\"add infos\")\nsub = pd.read_csv(f\"{ROOT}/sample_submission.csv\")\nsub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient','Weeks','Confidence','Patient_Week']]\nsub = sub.merge(chunk.drop('Weeks', axis=1), on=\"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = tr.append([chunk, sub])\n\nprint(tr.shape, chunk.shape, sub.shape, data.shape)\nprint(tr.Patient.nunique(), chunk.Patient.nunique(), sub.Patient.nunique(), \n      data.Patient.nunique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')\n\nbase = (\n    data\n    .loc[data.Weeks == data.min_week][['Patient','FVC', 'Percent']]\n    .rename({'FVC': 'base_FVC', 'Percent':'base_Percent'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FE = list(data.Sex.unique()) + list(data.SmokingStatus.unique())\ndata = pd.concat([\n    data,\n    pd.get_dummies(data.Sex),\n    pd.get_dummies(data.SmokingStatus)\n], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def Normalization(df):\n    \n    def get_fillness(series):\n        return (series - series.min()) / (series.max() - series.min())\n\n    df['Age'] = get_fillness(df['Age'])\n    df['base_FVC'] = get_fillness(df['base_FVC'])\n    df['base_week'] = get_fillness(df['base_week'])\n    df['base_Percent'] = get_fillness(df['base_Percent'])\n    \n    return df\n\nFE += ['Age','base_FVC','base_week','base_Percent']\ndata = Normalization(data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FE","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = data.loc[data.WHERE=='train']\nchunk = data.loc[data.WHERE=='val']\nsub = data.loc[data.WHERE=='test']\ndel data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr.shape, chunk.shape, sub.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Quantile Regression"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Activation >> Mish\nfrom tensorflow.keras.layers import Activation\nfrom tensorflow.keras.utils import get_custom_objects\n\nclass Mish(Activation):\n    def __init__(self, activation, **kwargs):\n        super(Mish, self).__init__(activation, **kwargs)\n        self.__name__ = 'Mish'\n\ndef mish(inputs):\n    return inputs * tf.math.tanh(tf.math.softplus(inputs))\n\nget_custom_objects().update({'Mish': Mish(mish)})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n    return K.mean(metric)\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    qs = [0.2, 0.50, 0.8]\n    q = tf.constant(np.array([qs]), dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return K.mean(v)\n#=============================#\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss\n#=================\ndef make_model(nh):\n    z = L.Input((nh,), name=\"Patient\")\n    x = L.Dense(100, activation=\"Mish\", name=\"d1\")(z)\n    x = L.Dense(100, activation=\"Mish\", name=\"d2\")(x)\n    p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n    model = M.Model(z, preds, name=\"NN\")\n    model.compile(loss=mloss(0.8),\n                  optimizer=tf.keras.optimizers.Adam(lr=0.01, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_cv_score(y_true, y_pred):\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n    sigma_clip = np.maximum(sigma, 70)\n    delta = np.abs(y_true[:, 0] - fvc_pred)\n    delta = np.minimum(delta, 1000)\n    sq2 = np.sqrt(2.)\n    metric = (delta / sigma_clip)*sq2 + np.log(sigma_clip* sq2)\n    return -np.mean(metric)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ncnt = 0\nBATCH_SIZE =256\nEPOCHS = 1500\nNFOLD = 11\n\nkf = GroupKFold(n_splits=NFOLD)\n\ny = tr['FVC'].values.astype('float32')\nz = tr[FE].values\nze = sub[FE].values\nnh = z.shape[1]\npe = np.zeros((ze.shape[0], 3))\npred = np.zeros((z.shape[0], 3))\n\n\nfor tr_idx, val_idx in kf.split(z, y, tr['Patient']):\n    cnt += 1\n    print(f\"FOLD {cnt}\")\n    net = make_model(nh)\n    \n    es = tf.keras.callbacks.EarlyStopping(monitor='val_loss',patience=200, min_delta = 0.000001,\n                                          verbose=1, mode='min')\n    lr_sch = tf.keras.callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.4, patience=50,\n                                                  verbose=0, mode='min', min_delta=0.000001, cooldown=0, min_lr=0)\n    net.fit(z[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS,\n            callbacks = [es, lr_sch], \n            validation_data=(z[val_idx], y[val_idx]), verbose=0)\n    \n    print(\"train\", net.evaluate(z[tr_idx], y[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"val\", net.evaluate(z[val_idx], y[val_idx], verbose=0, batch_size=BATCH_SIZE))\n    \n    print(\"predict val...\")\n    pred[val_idx] = net.predict(z[val_idx], batch_size=BATCH_SIZE, verbose=0)\n    print(calc_cv_score(y[val_idx].reshape(-1, 1), pred[val_idx]))\n\n    print(\"predict test...\")\n    pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD\n\n    #==============\nprint(\"CV SCORE\", calc_cv_score(y.reshape(-1, 1), pred))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Confidence Tuning"},{"metadata":{"trusted":true},"cell_type":"code","source":"import optuna\nfrom functools import partial\n\ntr['FVC_pred'] = pred[:, 1]\ntr['Confidence_pred'] = pred[:, 2] - pred[:, 0]\n\ndf_last_3 = tr.groupby('Patient').tail(3).reset_index(drop=True)\nX = df_last_3[['Weeks', 'FVC', 'FVC_pred', 'Confidence_pred']].values\nC = 0\n\ndef calc_tunned_score(y_true, y_pred, Conf):\n    sigma = Conf\n    fvc_pred = y_pred\n    sigma_clip = np.maximum(sigma, 70)\n    delta = np.abs(y_true - fvc_pred)\n    delta = np.minimum(delta, 1000)\n    sq2 = np.sqrt(2.)\n    metric = (delta / sigma_clip)*sq2 + np.log(sigma_clip*sq2)\n    return -np.mean(metric)\n\ndef objective(trial, X, y):\n    # create hyperparameters\n    a = trial.suggest_uniform('a', 0, 15)\n    b = trial.suggest_uniform('b', -100, 100)\n    \n    y = a * X[:, 0] + b\n    New_Confidence = X[:, 3] + y\n\n    # calculate score\n    return calc_tunned_score(X[:, 1], X[:, 2], New_Confidence)\n\nn_trials = 500\nobj = partial(objective, X=X, y=C)\nstudy = optuna.create_study(direction=\"maximize\")\noptuna.logging.disable_default_handler()\nstudy.optimize(obj, n_trials=n_trials)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('last 3 score befor tuning', calc_tunned_score(X[:, 1], X[:, 2], X[:, 3]))\nprint('last 3 score after tuning', study.best_value)\nparam = {k:v for k,v in study.best_params.items()}\nprint('param', param)\nprint('Training data score', calc_tunned_score(tr['FVC'].values, tr['FVC_pred'].values, tr['Confidence_pred'].values+param['a']*tr['Weeks'].values+param['b']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nunc = unc + (param['a'] * tr['Weeks'] + param['b'])\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(unc.min(), unc.mean(), unc.max(), (unc>=0).mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()\n\nplt.hist(unc, bins=30)\nplt.title(\"uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"# PREDICTION\nsub['FVC1'] = pe[:, 1]\nsub['Confidence1'] = pe[:, 2] - pe[:, 0]\nsubm = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1', 'Weeks']].copy()\n\nsubm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']\nif sigma_mean<70:\n    subm['Confidence'] = sigma_opt\nelse:\n    subm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1'] + param['a'] * subm.loc[~subm.FVC1.isnull(),'Weeks'] + param['b']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"otest = pd.read_csv(f\"{ROOT}/test.csv\")\n\nfor i in range(len(otest)):\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'FVC'] = otest.FVC[i]\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'Confidence'] = 0.1\n\nsubm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/carlossouza/bayesian-experiments\n\nsubm['Patient'] = subm['Patient_Week'].apply(lambda x:x.split('_')[0])\nsubm['Weeks'] = subm['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\n\ndef chart(df, patient_id, ax):\n    \n    plot_data = df[df['Patient'] == patient_id]\n    x = plot_data['Weeks']\n    FVC_low = plot_data['FVC'] - plot_data['Confidence']\n    FVC_high = plot_data['FVC'] + plot_data['Confidence']\n    \n    plot_data_tr = tr[tr['Patient'] == patient_id]\n    ax.plot(plot_data_tr['Weeks'], plot_data_tr['FVC'], 'o')\n    ax.plot(x, plot_data['FVC'])\n    ax.fill_between(x.values, FVC_low.values, FVC_high.values,\n                        alpha=0.5, color='#ffcd3c')\n    ax.set_title(patient_id)\n    ax.set_ylabel('FVC')\n    ax.set_ylim(min(FVC_low)-100, max(FVC_high)+100)\n\nf, axes = plt.subplots(2, 3, figsize=(15, 10))\nchart(subm, 'ID00419637202311204720264', axes[0, 0])\nchart(subm, 'ID00421637202311550012437', axes[0, 1])\nchart(subm, 'ID00422637202311677017371', axes[0, 2])\nchart(subm, 'ID00423637202312137826377', axes[1, 0])\nchart(subm, 'ID00426637202313170790466', axes[1, 1])","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}