{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook takes a look at the post processing from [this](https://www.kaggle.com/cdeotte/triple-stratified-kfold-with-tfrecords) notebook and also explores the use of test time augmentation (TTA).\n\nWe will explore the ways that the predictions from each model are combined, the use of predictions from TTA and whether we extract all possible information from the models that we spend many hours training. \n\nThe main point that will be pursued is a method for encoding the uncertainty each model has in a given data point and using this information while ensembling.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Import everything\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport scipy.stats\nimport matplotlib.pyplot as plt\nfrom scipy.stats import spearmanr\nfrom matplotlib.colors import LogNorm\n\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.model_selection import StratifiedKFold\nimport tensorflow as tf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\n!cp ../input/rapids/rapids.0.14.0 /opt/conda/envs/rapids.tar.gz\n!cd /opt/conda/envs/ && tar -xzvf rapids.tar.gz > /dev/null\nsys.path = [\"/opt/conda/envs/rapids/lib/python3.7/site-packages\"] + sys.path\nsys.path = [\"/opt/conda/envs/rapids/lib/python3.7\"] + sys.path\nsys.path = [\"/opt/conda/envs/rapids/lib\"] + sys.path \n!cp /opt/conda/envs/rapids/lib/libxgboost.so /opt/conda/lib/\nimport cuml #, cupy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Setting color palette.\norange_black = [\n    '#fdc029', '#df861d', '#FF6347', '#aa3d01', '#a30e15', '#800000', '#171820'\n]\n\n# Setting plot styling.\nplt.style.use('ggplot')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I have trained a network using [this](https://www.kaggle.com/cdeotte/triple-stratified-kfold-with-tfrecords) notebook and saved all of the TTA predictions for both out of fold (OOF) points and the test set that each model makes.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# For this model I augmented 11 images at test time\nn_tta = 1500\nfile = 'tta-exploration-128x128-b0'\n\n# Load all of the TTA predictions for the test set.\nall_tta_test = pd.read_csv('../input/'+file+'/all_tta_test.csv')\nall_tta_test.columns.values[0] = 'image_name'\nsubmission = pd.read_csv('../input/'+file+'/submission.csv')\n\n# Load the oof df where the TTA predictions are contained in the numbered columns.\ntta_keys = [str(i) for i in range(n_tta)]\noof = pd.read_csv('../input/'+file+'/oof.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first thing we will look at is the distribution of the standard deviation across the mean TTA predictions. We take the mean of the TTA predictions for each example, as well as the standard deviation. \n\nExamples with a small standard deviation across TTA predictions are points that we can consider to sit in an area of the input distribution that each model is confident about. Those that have a large standard deviation are unstable to changes (augmentations) that we have decided do not change the class of that example, and therefore come from a part of the input distribution that is close to the decision boundary. This uncertainity is what we will try to leverage.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def IQR(data,ax=1):\n    return np.subtract(*np.percentile(data, [75, 25],axis=ax))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tta_preds = oof[tta_keys]\nmn_pred = np.mean(tta_preds,axis=1)\nstd_pred = np.std(tta_preds,axis=1)\n\n\nplt.figure(figsize=(15,5))\nplt.subplot(1,2,1)\nmxb = oof['target']==0\nh = plt.hist2d(mn_pred[mxb],std_pred[mxb],norm=LogNorm(),bins=100)\nplt.colorbar(h[3])\nplt.xlabel('Arithmetic mean')\nplt.ylabel('Standard deviation')\nplt.title('Benign')\nplt.subplot(1,2,2)\nmxm = oof['target']==1\nh = plt.hist2d(mn_pred[mxm],std_pred[mxm],norm=LogNorm(),bins=30)\nplt.colorbar(h[3])\nplt.xlabel('Arithmetic mean')\nplt.ylabel('Standard deviation')\nplt.title('Malignant')\nplt.show()\n\nplt.figure(figsize=(8,8))\nplt.scatter(mn_pred[mxb],std_pred[mxb],color='orange',s=10,label='Benign')\nplt.scatter(mn_pred[mxm],std_pred[mxm],color='blue',s=10,label='Malignant')\nplt.xlabel('Arithmetic mean')\nplt.ylabel('Standard deviation')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These plots show that there is a wide range in the uncertainties that each model has for different means, if the above points traced out a parabola with no width then this would be of no use because the relationship between mean and standard deviation would be deterministic - this would mean that all of the information in the standard deviation is captured by the mean, and vice versa. As this relationship is not deterministic, we might be able to gain some leverage from it, and we will explore this.\n\nI would like to note that I do not think this extra information will be useful on its own, but rather that it will be useful when ensembling multiple different models. The reason that I think this is the case can be explained through an example. The model shown above is relatively confident about the points that run along the bottom of the bell. Now imagine we have another model, where those same points run along the top of the bell, but the predictions for each point differ. When blending these models I would like to place more weight on the confident predictions from the former model than on the unconfident predictions of the latter (with some weighting also given to each models accuracy).\n\nThe incorporation of the standard deviation into an ensemble prediction can be done in many ways, and we will explore just one below.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Distribution and error of TTA predictions\nThe first thing we need to know is how the TTA predictions are distributed, how many samples we need to gather to get accurate statistics on that distribution, and what measures we should use to quantify the uncertainty in this distribution. \n\nNote that to be safe I have not assumed I have enough examples to do bootstrap sampling, even though I do for low TTA numbers.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"nplot = 5\n\nfig, axs = plt.subplots(nplot, nplot,figsize=(10,10))\nfig.suptitle(\"Distribution of TTA predictions per example\", fontsize=14)\nfor i,ax in enumerate(axs.flatten()):\n    ax.hist(tta_preds.iloc[i,:],bins=50,color=orange_black[3])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of these distributions are skewed, but I do not think that looking at quartiles will provide any more information because the data is always skewed towards either zero or one, and the amount that it is skewed depends on its distance from zero or one and so the information about how the predictions are skewed is already contained in the mean.\n\nThe use of quartiles would be interesting to check. For now I will keep things simple and look at the standard deviation.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Look at the range and standard deviation for the TTA predictions.\nranges = np.max(tta_preds,axis=1)-np.min(tta_preds,axis=1)\niqr = IQR(tta_preds)\n\nfig, axs = plt.subplots(1, 3,figsize=(20,5))\naxs[0].hist(ranges,bins=100)\naxs[0].set_title('Range')\naxs[1].hist(std_pred,bins=100)\naxs[1].set_title('Standard deviations')\naxs[2].hist(iqr,bins=100)\naxs[2].set_title('Interquartile range')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The tail of these distributions are quite long, suggesting that they should be more heavily sampled to get a more accurate prediction for these data points. The other data points do not need as many samples for an accurate mean to be calculated. Such a method would be useful for evaluation even if the standard deviation proves to be uninteresting.\n\nNext I want to look at how many TTA predictions I need to get an accurate estimate of the mean and standard deviation for each data point. To do so we look at the standard errors for both of these statistics.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"tta_arr = np.array(tta_preds)\nn_samp = tta_arr.shape[1]\n# mn_pred, std_pred\naverage_t=[];max_t=[];std_t=[]\nerror_std=[]\n# If you drop the relative part you can run this to n_samp, but that is less informative.\nfor TTA in range(2,100):\n    end = int(n_tta//TTA)\n    samples = tta_arr[:,:TTA*end].reshape(-1,TTA)\n    all_avs = np.mean(samples,axis=1).reshape(len(tta_arr),end)\n    if all_avs.shape[1]>10:\n        all_stds = np.std(samples,axis=1).reshape(len(tta_arr),end)\n        all_e_std = np.std(all_stds,axis=1)\n        # error_std+=[np.mean(all_e_std)]\n        # Make it the relative standard error of the standard deviations\n        error_std+=[np.mean(all_e_std/std_pred)]\n    if all_avs.shape[1] > 30:\n        # Calculate the standard errors directly.\n        std_errors = np.std(all_avs,axis=1)\n    else:\n        # Estimate the standard error per example and then take the average as the best estimate.\n        std_errors = np.mean(all_stds/np.sqrt(TTA),axis=1)\n    # Take the relative standard error of the mean\n    std_errors = std_errors/mn_pred\n    std_t+=[np.std(std_errors)]\n    average_t+=[np.mean(std_errors)]\n    max_t+=[np.max(std_errors)]\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(1, 4,figsize=(20,5))\naxs[0].plot(average_t)\naxs[0].set_title('Mean standard error')\naxs[0].set_xlabel('TTA')\naxs[1].plot(std_t)\naxs[1].set_title('Std of standard errors')\naxs[1].set_xlabel('TTA')\naxs[2].plot(max_t)\naxs[2].set_title('Max standard error')\naxs[2].set_xlabel('TTA')\n# axs[3].plot([int(n_tta//TTA) for TTA in range(2,len(max_t))])\n# axs[3].set_title('N samples')\n# axs[3].set_xlabel('TTA')\naxs[3].plot(error_std)\naxs[3].set_title('Mean standard error of std')\naxs[3].set_xlabel('TTA')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The estimate of the mean converges much faster than the estimate of the standard deviation. That leaves the question of whether the standard deviation can be successcully used for training when using models that do not have a large amount of TTA. An importance sampling algorithm would overcome this issue. Unfortunately I do not have this implemented at present.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Performance of different means\nBefore moving on to stacking we will look at the use of different means for extracting a prediction from the TTA predictions. The distributions for different TTA are heavily skewed and so the best summary statistic is not obvious.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"tta_preds = oof[tta_keys]\nmns = dict(\n    amean = np.mean(tta_preds,axis=1),\n    gmean = scipy.stats.gmean(tta_preds,axis=1),\n    median = np.median(tta_preds,axis=1)\n)\n\nfor mn in mns:\n    print('{} AUC = {}'.format(mn,roc_auc_score(oof['target'],mns[mn])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The median is the best summary here, which is expected as the distributions are skewed. But for small sample sizes (the test set only has 11 TTA predictions) the arithmetic mean is more stable, so we will use that.\n\nTo explore the methods of combining the predictions for submission from the different models the best we can do in this case is to look at the private leaderboard scores. As the competition is over we no longer have to worry about overfitting to this and it can be considered a held out test set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_sub(func,name='submission'):\n    submission = pd.DataFrame(dict(image_name=test_tta['0'], target=func(preds_all,axis=1)))\n    submission.to_csv(name+'.csv', index=False)\n    submission.head()\n\ntest_tta = pd.read_csv('../input/tta-exploration-128x128-b0/all_tta_test.csv')\nFOLDS=3; n_test_tta=11\npreds_all = [np.mean(test_tta[[str(i) for i in range(j*n_test_tta,(j+1)*n_test_tta)]],axis=1) for j in range(FOLDS)]\npreds_all = np.vstack(preds_all).T\n\nmake_sub(np.mean,name='amean_submission')\nmake_sub(scipy.stats.gmean,name='gmean_submission')\nmake_sub(np.median,name='median_submission')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Blend/Stack models\nThis is the key part of this notebook, can the incorporation of the uncertainty of base models increase the OOF CV score more than blending with linear weights?\n\nIn the following we will load the predictions from two models, find the optimal blend using weights and the means, and also train a neural network that takes as input their means and std.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load a new model\nn_tta = 11\n# n_tta = 1500\nFOLDS=3\nfiles = ['tta-exploration-128x128-b0','tta-exploration-128x128-b3']\n# files = ['tta-exploration-128x128-b0','tta-exploration-128x128-b3','meta-384-b5-coarse','meta-384-b3']\n\n# Load all of the TTA predictions for the test set.\nall_tta_list = [pd.read_csv('../input/'+file+'/all_tta_test.csv') for file in files]\nfor df in all_tta_list:\n    df.columns.values[0] = 'image_name' \nsub_list= [pd.read_csv('../input/'+file+'/submission.csv') for file in files]\n\n# Load the oof df where the TTA predictions are contained in the numbered columns.\ntta_keys = [str(i) for i in range(n_tta)]\noof_list = [pd.read_csv('../input/'+file+'/oof.csv') for file in files]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Add the std as features and then drop the tta keys\ndef set_right(df,tta_keys):\n    df['std']=np.std(df[tta_keys],axis=1)\n    return df.drop(tta_keys,axis=1)\n\noof1 = [set_right(df,tta_keys).set_index('image_name') for df in oof_list]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Assemble\ndf_oof = pd.concat(oof1,axis=1)\nlabels = df_oof['target'].iloc[:,0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualize the correlation between the two models.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_scatter(ax,x,y,mx1,mx2,title):\n    ax.scatter(x[mx1],y[mx1],color='orange',s=10,label='Benign')\n    ax.scatter(x[mx2],y[mx2],color='blue',s=10,label='Malignant')\n    ax.set_title(title)\n    ax.legend()\n    \n    \nfig, axs = plt.subplots(2, 2,figsize=(20,15))\nmxb = labels==0\nmxm = labels==1\nstd1 = df_oof['std'].iloc[:,0]; pred1 = df_oof['pred'].iloc[:,0]\nstd2 = df_oof['std'].iloc[:,1]; pred2 = df_oof['pred'].iloc[:,1]\nplot_scatter(axs[0,0],std1,std2,mxb,mxm,'Standard Deviations')\nplot_scatter(axs[0,1],pred1,pred2,mxb,mxm,'Predictions')\nplot_scatter(axs[1,0],pred1,std2,mxb,mxm,'Prediction 1 vs Standard Deviation 2')\nplot_scatter(axs[1,1],std1,pred2,mxb,mxm,'Standard Deviation 1 vs Prediction 2')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"corr_disp = df_oof[['pred','std']].copy()\ncorr_disp.columns = ['p1','p2','s1','s2']\n\ncorr=corr_disp.astype(float).corr()\ncorr.style.background_gradient(cmap='coolwarm').set_precision(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Look at the model spaces when including std as a feature with tsne\n\nmodel = cuml.TSNE()\nembed2D = model.fit_transform(df_oof[['pred','std']])\n\nfig, axs = plt.subplots(1,2,figsize=(20,8))\nmxb = labels==0\nmxm = labels==1\nplot_scatter(axs[0],embed2D[:,0],embed2D[:,1],mxb,mxm,'Predictions and Uncertainty')\n\n# Look at the model spaces when including std as a feature with tsne\n\nmodel = cuml.TSNE()\nembed2D = model.fit_transform(df_oof[['pred']])\nplot_scatter(axs[1],embed2D[:,0],embed2D[:,1],mxb,mxm,'Predictions Only')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The clustering in both examples looks approximately the same, so we cannot say much about which will be better without fitting a model and checking.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Linear Blend","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Find the best weights\n\nroc=[]\nWGTS=np.linspace(0,1,11)\nfor wgt in WGTS:\n    pred = wgt*df_oof['pred'].iloc[:,0] + (1-wgt)*df_oof['pred'].iloc[:,1]\n    roc+=[roc_auc_score(labels,pred)]\nbest_roc = max(roc)\nloc = np.where(roc==best_roc)[0][0]\nbest_weight = WGTS[loc]\nprint('The best roc = {:.6f} with weights : ({},{})'.format( best_roc, best_weight, (1-best_weight)) )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Stack Models ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from xgboost import XGBClassifier","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ids = df_oof.index\nfold_finder = df_oof['fold'].iloc[:,0]\n\ndef get_oof(DATA,LABELS,folds=10,PRINT=0,SEED=42):\n    skf = StratifiedKFold(n_splits=FOLDS,shuffle=True,random_state=SEED)\n    oof_pred = []; oof_tar = []; oof_val = []; oof_names = []; oof_folds = [] \n\n    for fold,(idxT,idxV) in enumerate(skf.split(DATA,LABELS)):\n    #  The below will use the same splits as the original data, but using larger folds reduces the dependence on SEED.\n    # for fold in range(FOLDS):\n    #     idxT = fold_finder!=fold; idxV = fold_finder==fold;\n\n        X = DATA[idxT]; y = LABELS[idxT]\n        X_val = DATA[idxV]; y_val = LABELS[idxV]\n\n        # Train the model\n        # Note that the optimal parameters for this model were found on the data that did not include the standard\n        # deviation\n        model = XGBClassifier(  objective='binary:logistic',\n                            seed=0,  \n                            nthread=-1, \n                            learning_rate=0.01, \n                            n_estimators=1200,\n                            scale_pos_weight=1,\n                            min_child_weight=1,\n                            max_depth=2,\n                            subsample=0.7,\n                            colsample_bytree=0.9,\n                            gamma=5,\n                            random_state=42,\n                         )\n        model.fit(X,y)\n\n    #     model.load_weights('fold-%i.h5'%fold)\n\n        # PREDICT OOF USING TTA\n        pred = model.predict_proba( X_val )[:,1]\n        oof_pred.append( pred )                \n\n        # GET OOF TARGETS AND NAMES\n        oof_tar.append( y_val )\n        oof_names.append( ids[idxV] )\n        oof_folds.append( np.ones_like(oof_tar[-1],dtype='int8')*fold )\n\n        # REPORT RESULTS\n        auc = roc_auc_score(oof_tar[-1],oof_pred[-1])\n        if PRINT:\n            print('#### FOLD %i OOF AUC = %.3f'%(fold+1,auc))\n\n\n    # COMPUTE OVERALL OOF AUC\n    oof = np.concatenate(oof_pred); true = np.concatenate(oof_tar);\n    names = np.concatenate(oof_names); folds = np.concatenate(oof_folds)\n    auc = roc_auc_score(true,oof)\n    print('Overall OOF AUC with TTA = %.5f'%auc)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# model.feature_importances_\nFOLDS = 20\nprint('Including confidence estimate:')\nget_oof(np.concatenate((df_oof['pred'],df_oof['std']),axis=1),labels,folds=FOLDS)\nprint('\\nOnly predictions:')\nget_oof(np.array(df_oof['pred']),labels,folds=FOLDS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we can see that including the standard deviation provides a slight improvement in the OOF AUC score when stacking, but it does not beat straightforward blending. The stacking and blending scores cannot be properly explored without further experiments. However, this idea can be explored more robustly in less expensive settings where statistics can be gathered easily and tests can be run where more models are included (the effect of including the standard deviation could grow or vanish as more models are added to the stacking). \n\nFor now, this is a moderate success that implies it is worth pursuing further.","execution_count":null}],"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}