{"cells":[{"metadata":{},"cell_type":"markdown","source":"First things first I want to give credit where it's due. Thank you Ulrich G. for providing the bulk of this notebook in https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter . The dataset is at https://www.kaggle.com/eladwar/conditionalrnn"},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\nsys.path.insert(0,\"../input/conditionalrnn\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport os\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold,GroupKFold\n\nimport tensorflow as tf\nfrom cond_rnn import ConditionalRNN\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\nfrom tensorflow.keras.optimizers import Adam,Nadam\nfrom tensorflow.keras import initializers","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\ndef 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(43)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NFOLD = 5\nkf = KFold(n_splits=NFOLD)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = 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])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(tr.shape, chunk.shape, sub.shape, data.shape)","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')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)","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":{},"cell_type":"markdown","source":"# Categorical Features for RNN:\n* The commented categorical feature generation works for better results in the competition (Notebook V4)\n* The uncommented categorical feature generation is used particularly for less repetition (in other words, I don't want to One-hot Encode)\n* You can make two models focusing on each set of generated features and ensemble \n* You can also make a model using all possible ways to generate the categorical features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create Alternate Categorical Features \nCOLS = ['Sex','SmokingStatus']\nFE = []\nfor col in COLS:\n    data[col] = pd.factorize(data[col])[0]\nFE.extend(COLS)\n#=================\nCOLS = ['Sex','SmokingStatus']\n# FE = []\nfor mod in data[col].unique():\n    FE.append(mod)\n    data[mod] = (data[col] == mod).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Scale features to quantiles \ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\ndata['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min() )\ndata['week'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\ndata['percent'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )\nFE += ['age','percent','week','BASE']","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":"print(tr.shape, chunk.shape, sub.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What do I do to get everything ready for a categorical RNN\n* Split categorical and measurement features\n* Measurement features are reshaped into Recurent Neural Network Standard Format: [batch, timesteps, features]\n* Category features are left in a normal shape"},{"metadata":{"trusted":true},"cell_type":"code","source":"tr[FE]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get Features\n# groups = tr['Patient']\ny = tr['FVC'].values\nz = tr[FE].values\nze = sub[FE].values\n# --------------------------------------------------------\n#Create oof and prediction arrays\npe = np.zeros((ze.shape[0], 3))\npred = np.zeros((z.shape[0], 3))\n# --------------------------------------------------------\n#Split the train set's categorical and measurement features\ntrain_categories = z[:,0:5]#normal shape\ntrain_measurements = z[:,5:].reshape(z.shape[0],1,-1)#[batch, timesteps, features]\n# --------------------------------------------------------\n#Split the test set's categorical and measurement features\ntest_categories = ze[:,0:5] #normal shape\ntest_measurements = ze[:,5:].reshape(ze.shape[0],1,-1)#[batch, timesteps, features]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Loss & Scoring Functions"},{"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    \n    #sigma_clip = sigma + C1\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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conditional RNN model"},{"metadata":{"trusted":true},"cell_type":"code","source":"class MySimpleModel(tf.keras.Model):\n    def __init__(self):\n        super(MySimpleModel, self).__init__()\n        self.cond = ConditionalRNN(100, cell='LSTM', dtype=tf.float32)\n        self.x = L.Dense(100, activation='elu',kernel_initializer='he_uniform', name=\"d2\")\n        self.p1 = L.Dense(3, activation=\"linear\", name=\"p1\")\n        self.p2 = L.Dense(3, activation=\"relu\",kernel_initializer='he_uniform', name=\"p2\")\n        self.preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                         name=\"preds\")\n\n    def call(self, inputs, **kwargs):\n        o = self.cond(inputs)\n        o = self.x(o)\n        linear = self.p1(o)\n        relu = self.p2(o)\n        o = self.preds([linear,relu])\n        return o","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"net = MySimpleModel()\nnet.call([train_measurements,train_categories])\nnet.compile(optimizer='adam', loss=mloss(0.8), metrics=[score])\n# optimizer = tf.keras.optimizers.Adam(learning_rate=0.001)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Just double checking model inpute\nprint(train_measurements.shape,train_categories.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting jit/xla for greater efficiency \ntf.keras.backend.clear_session()\ntf.config.optimizer.set_jit(True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ncnt = 0\nBatchSize = 200\nfor tr_idx, val_idx in kf.split(z):\n    cnt += 1\n    \n    print(f\"FOLD {cnt}\")\n#     with tf.device('/gpu:0'):\n    net.fit([train_measurements[tr_idx,:,:],train_categories[tr_idx,:]], pd.Series(y[tr_idx].astype(float).flatten()), batch_size=BatchSize, epochs=800, \n    validation_data=([train_measurements[val_idx,:,:],train_categories[val_idx,:]], pd.Series(y[val_idx].astype(float).flatten())), verbose=0) #\n    print(\"train\", net.evaluate([train_measurements[tr_idx,:,:],train_categories[tr_idx,:]], y[tr_idx], verbose=0, batch_size=BatchSize))\n    print(\"val\", net.evaluate([train_measurements[val_idx,:,:],train_categories[val_idx,:]], y[val_idx], verbose=0, batch_size=BatchSize))\n    print(\"predict val...\")\n    pred[val_idx] = net.predict([train_measurements[val_idx,:,:],train_categories[val_idx,:]], batch_size=BatchSize, verbose=0)\n    print(\"predict test...\")\n    pe += net.predict([test_measurements,test_categories], batch_size=BatchSize, verbose=0) / NFOLD\n# ==============","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:,1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This simply checks what quantile reveals the best results, if you want to do some extra data transformation toward the end to boost your score"},{"metadata":{"trusted":true},"cell_type":"code","source":"increment = 0.001\nerror = []\nfor i in np.arange(0,1,increment):\n    quant_5 = np.quantile(pred,i,axis=1)\n#     print(mean_absolute_error(y,quant_5))\n    error.append(mean_absolute_error(y,quant_5))\n    \nprint('Best Quantile:',np.arange(0,1,increment)[np.argmin(error)])\nprint('Best MAE of Optimized Quantile:',error[np.argmin(error)])\nprint('Baseline MAE:', mean_absolute_error(y,pred[:,1]))\n# Optimized vs Baseline Graphs\nplt.plot(quant_5,c='g')\nplt.title('Optimized vs Baseline')\nplt.plot(pred[:,1])\nplt.plot(y)\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import math\nimport scipy as sp\nfrom functools import partial\nscoring_df = pd.DataFrame(z)\nscoring_df['FVC_pred'] = pred[:,1]\nscoring_df['FVC'] = y\n# baseline score\nscoring_df['Confidence'] = 100\nscoring_df['sigma_clipped'] = scoring_df['Confidence'].apply(lambda x: max(x, 70))\nscoring_df['diff'] = abs(scoring_df['FVC'] - scoring_df['FVC_pred'])\nscoring_df['delta'] = scoring_df['diff'].apply(lambda x: min(x, 1000))\nscoring_df['score'] = -math.sqrt(2)*scoring_df['delta']/scoring_df['sigma_clipped'] - np.log(math.sqrt(2)*scoring_df['sigma_clipped'])\nscore = scoring_df['score'].mean()\nprint(score)\n\ndef loss_func(weight, row):\n    confidence = weight\n    sigma_clipped = max(confidence, 70)\n    diff = abs(row['FVC'] - row['FVC_pred'])\n    delta = min(diff, 1000)\n    score = -math.sqrt(2)*delta/sigma_clipped - np.log(math.sqrt(2)*sigma_clipped)\n    return -score\n\nresults = []\ntk0 = tqdm(scoring_df.iterrows(), total=len(scoring_df))\nfor _, row in tk0:\n    loss_partial = partial(loss_func, row=row)\n    weight = [100]\n    result = sp.optimize.minimize(loss_partial, weight, method='SLSQP')\n    x = result['x']\n    results.append(x[0])\n\n# optimized score\nscoring_df['Confidence'] = results\nscoring_df['sigma_clipped'] = scoring_df['Confidence'].apply(lambda x: max(x, 70))\nscoring_df['diff'] = abs(scoring_df['FVC'] - scoring_df['FVC_pred'])\nscoring_df['delta'] = scoring_df['diff'].apply(lambda x: min(x, 1000))\nscoring_df['score'] = -math.sqrt(2)*scoring_df['delta']/scoring_df['sigma_clipped'] - np.log(math.sqrt(2)*scoring_df['sigma_clipped'])\nscore = scoring_df['score'].mean()\nprint(score)","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()","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":"plt.hist(unc)\nplt.title(\"Difference between 20th and 80th Quantiles\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### PREDICTION"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC1'] = pe[:,1]\nsub['Confidence1'] = (pe[:, 2] - pe[:, 0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.loc[~subm.FVC1.isnull()].head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.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']","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":"subm[\"Confidence\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission.csv\", index=False)","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}