{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport lightgbm as lgbm\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nfrom sklearn.metrics import roc_curve, auc, classification_report, confusion_matrix\n\ntrain = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/train.csv')\ntest = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Adversarial Validation\n\nIt has come to some Kagglers attention that the CV scores do not match the LB scores, this is the main discussion thread (https://www.kaggle.com/c/predict-volcanic-eruptions-ingv-oe/discussion/192766). There are couple reasons this might be: many models are overfit, the train sets and test sets are not a good match, there are some outlier examples in train or test that might skew our models, there are some specific features generated that do match between train and test. Here we just try to determine/confirm a couple of these hypothesis using Adversarial Validation. We are simply trying to determine if train and test are good matches, to do so we are going to \"hide\" our test segments inside the train segments and have a model try to \"find them\". So, to clarify, if the two datasets are statistically similar it should be difficult for a model to find the test segments within train, if they are different it should be much easier.\n\n#### Methodology\nMy previous notebook has been a popular choice for initial starter, some people spinning off additional notebooks and adding features and oversampling to improve results, but for this excercise I will just take the features generated in that notebook (https://www.kaggle.com/ajcostarino/ingv-volcanic-eruption-prediction-lgbm-baseline). I will then discard the `time_to_eruption` column, as that is not important for this task. I set a new target column, called `target`, to `0` for the train set and `1` for the test set. We then build a model to try and see if the the sets are distinguishable.\n\n#### Calibration\nHere I decided to use `binary_logloss` for my evaluation metric because I would prefer if the model was well calibrated. What does that mean? It means I want to make sure that the output probabilities match real-life, that is, if the model says a particular segment has `0.57` chance of being a test segment I would like that to be true 57% of the time. Logloss penalizes on confidence, so we can look at our segments and see how confident we are that they are test or train. The test and train sets match if the model spits out roughly `0.5` probability for every segment, if that happens then we are validated that the test and train sets are a good match. The classes are balanced, there are an equal amount of each. For our `binary_logloss` metric a model that cannot determine test or train, i.e a model that is random, will return a `binary_logloss` of about `0.7` (`0.69` to be more exact). So if our model is able to get a `binary_logloss` score lower than that we know that the model is not random and that there are differences between the test and train set that might be decisive and might lead to our CV LB differences.\n\nUpdated by making dataset public."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = train.drop(['time_to_eruption'], axis = 1)\n\n# Set target of train to 0 and test to 1\n# We are trying to find the test segments if we hide them in train\ntrain['target'] = 0\ntest['target'] = 1\n\nall_segments = pd.concat([train, test])\nall_segments = all_segments.set_index(['segment_id'])\nall_segments.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_roc_auc(cv_scores, cutoff):\n    cv_scores = cv_scores.iloc[:cutoff]\n    fpr0, tpr0, thresholds1 = roc_curve(cv_scores['fold_0_val'], cv_scores['fold_0_predict_proba'])\n    fpr1, tpr1, thresholds2 = roc_curve(cv_scores['fold_1_val'], cv_scores['fold_1_predict_proba'])\n    fpr2, tpr2, thresholds3 = roc_curve(cv_scores['fold_2_val'], cv_scores['fold_2_predict_proba'])\n    fpr3, tpr3, thresholds4 = roc_curve(cv_scores['fold_3_val'], cv_scores['fold_3_predict_proba'])\n\n    roc_auc0 = auc(fpr0, tpr0)\n    roc_auc1 = auc(fpr1, tpr1)\n    roc_auc2 = auc(fpr2, tpr2)\n    roc_auc3 = auc(fpr3, tpr3)\n\n    plt.figure()\n    lw = 2\n\n    fig, ax = plt.subplots(figsize = (11, 10))\n    plt.plot(fpr0, tpr0, color='red',\n             lw=lw, label='Fold 0 ROC curve (area = %0.2f)' % roc_auc0)\n    plt.plot(fpr1, tpr1, color='green',\n             lw=lw, label='Fold 1 ROC curve (area = %0.2f)' % roc_auc1)\n    plt.plot(fpr2, tpr2, color='orange',\n             lw=lw, label='Fold 2 ROC curve (area = %0.2f)' % roc_auc2)\n\n    plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\n    plt.xlim([-0.01, 1.0])\n    plt.ylim([0.0, 1.05])\n    plt.xlabel('False Positive Rate')\n    plt.ylabel('True Positive Rate')\n    plt.title('Receiver operating characteristic example')\n    plt.legend(loc=\"lower right\")\n    plt.show()\n\n    \ndef feature_importances(model, index):\n    feature_importances = pd.DataFrame(model.feature_importances_,\n                                       index = index,\n                                       columns=['importance']).sort_values('importance',\n                                       ascending=False)\n\n    feature_importances = feature_importances.reset_index()\n    feature_importances.columns = ['feature', 'importance']\n\n    fig, ax = plt.subplots(figsize = (18, 24))\n    sns.set()\n    plt.subplot(1, 1, 1);\n    sns.barplot(x=\"importance\", y=\"feature\", orient='h', data=feature_importances.head(50));\n    plt.title('Feature Importance');\n    \n\ndef plot_lift_chart(df,\n                    actual_col,\n                    predicted_col,\n                    probability_col,\n                    ntiles):\n\n    df = df.sort_values(by=probability_col, ascending=False)\n\n    rows = []\n    for group in np.array_split(df, ntiles):\n        score = group[actual_col].sum()\n        rows.append({'cases': len(group), 'correct': score})\n\n    lift = pd.DataFrame(rows)\n\n    #Cumulative Gains Calculation\n    lift['cumcorrect'] = lift['correct'].cumsum()\n    lift['ntile'] =  (float(len(df)) / float(ntiles)) / float(len(df)) * 100\n    lift['cumntile'] = lift['ntile'].cumsum()\n    lift['avgcorrect'] = df[actual_col].sum() / ntiles\n    lift['cumavgcorrect'] = lift['avgcorrect'].cumsum()\n\n    #Lift Chart\n    lift['normalisedpercentavg'] = 1\n    lift['normalisedpercentmodel'] = lift['cumcorrect'] / lift['cumavgcorrect']\n\n    lift = lift.set_index('cumntile')\n\n    #fig, ax = plt.subplots()\n    \n    return lift\n\n#     ax.plot(lift['normalisedpercentavg'], 'r-', label='Normalised \\'response rate\\' with no model')\n#     ax.plot(lift['normalisedpercentmodel'], 'g-', label='Normalised \\'response rate\\' with using model')\n#     ax.set_xlim((-0.5,100))\n#     ax.legend()\n#     return ax","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_lightgbm(X, y, n_fold):\n    folds = KFold(n_splits=n_fold,\n                  shuffle=True,\n                  random_state=42)\n    \n    cv_scores = pd.DataFrame()\n    \n    params = {\n        'num_leaves': 85,\n        'min_data_in_leaf': 10, \n        'objective':'regression',\n        'max_depth': -1,\n        'learning_rate': 0.001,\n        'max_bins': 2048,\n        'boosting': \"gbdt\",\n        'feature_fraction': 0.91,\n        'bagging_freq': 1,\n        'bagging_fraction': 0.91,\n        'bagging_seed': 42,\n        'metric': 'binary_logloss',\n        'lambda_l1': 0.1,\n        'verbosity': -1,\n        'nthread': -1,\n        'random_state': 42\n    }\n\n    fold = 0\n    for fold_, (trn, val) in enumerate(folds.split(X, y)):\n        \n        strLog = \"fold {}\".format(fold_)\n        print(strLog)\n        \n        X_trn, X_val = X.iloc[trn], X.iloc[val]\n        y_trn, y_val = y.iloc[trn], y.iloc[val]\n        \n        lgbmc = lgbm.LGBMClassifier(**params, n_estimators = 2500, n_jobs = -1)\n        \n        eval_set = [\n            (X_trn, y_trn),\n            (X_val, y_val)\n        ]\n        \n        lgbmc.fit(X_trn, y_trn, eval_set=eval_set, verbose=200, early_stopping_rounds=400)\n        \n        \n        y_predict       = lgbmc.predict(X_val)\n        y_predict_proba = lgbmc.predict_proba(X_val)[:,1]\n\n        segment_col       = 'fold_{}_segment_id'.format(fold)\n        predict_col       = 'fold_{}_predict'.format(fold)\n        predict_proba_col = 'fold_{}_predict_proba'.format(fold)\n        val_col           = 'fold_{}_val'.format(fold)\n    \n        cv_scores[segment_col]       = pd.Series(list(X_val.index))\n        cv_scores[predict_col]       = pd.Series(y_predict)\n        cv_scores[predict_proba_col] = pd.Series(y_predict_proba)\n        cv_scores[val_col]           = pd.Series(list(y_val))\n    \n        fold += 1\n        \n    return cv_scores, lgbmc, X_trn.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model\nSo we can see here clearly that our model can determine some differences between test and train. However, `0.57` is not a great `binary_logloss` score so many of the train and test samples must be similar. There might be a subsegment of train and test that are recognizable. We will try and determine which ones those are."},{"metadata":{"trusted":true},"cell_type":"code","source":"cv_scores, lgbmc, columns = build_lightgbm(all_segments.drop(['target'], axis = 1), all_segments['target'], 5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## ROC/AUC and Lift Curves\nThe ROC/AUC curve allows us to visualize what random looks like as well as what the output probabilities look like. We can see an extremely steep initial curve meaning that there must be a subset of the test set that are clearly distibguishable from the train set. After that the curve looks relatively linear meaning that the rest of the test set is indisinguishable from train. In the lift charts we see similar behavior. The highest confident predictions show that the model is able to disinguish a subset of the test observations, the must be much different from train. However it is good that after that small subset our train and test are more distinguishable."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_roc_auc(cv_scores, 1789)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(20, 10))\n\ndef plot_lift(ax, lift):\n    ax.plot(lift['normalisedpercentavg'], 'r-', label='Normalised \\'response rate\\' with no model')\n    ax.plot(lift['normalisedpercentmodel'], 'g-', label='Normalised \\'response rate\\' with using model')\n    ax.set_xlim((-0.5,100))\n    ax.legend()\n\nplot_lift(axes[0, 0], plot_lift_chart(cv_scores,'fold_0_val','fold_0_predict','fold_0_predict_proba',100))\nplot_lift(axes[0, 1], plot_lift_chart(cv_scores,'fold_1_val','fold_1_predict','fold_1_predict_proba',100))\nplot_lift(axes[1, 0], plot_lift_chart(cv_scores,'fold_2_val','fold_2_predict','fold_2_predict_proba',100))\nplot_lift(axes[1, 1], plot_lift_chart(cv_scores,'fold_3_val','fold_3_predict','fold_3_predict_proba',100))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Important Adversarial Features\n\nWhich features are the most important in descriminating between the datasets. These features are potentially bad features as they might lead to some inferences in the train set that might not be true in the test set. We can also plot the distributions of the top features to see what the differences are between the test and train set."},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importances(lgbmc, columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_distribution_curves(ax, all_segments, feature, xlim):\n    sns.distplot(all_segments[all_segments['target'] == 0][feature], 30, ax = ax)\n    sns.distplot(all_segments[all_segments['target'] == 1][feature], 30, ax = ax)\n    ax.set_title(f'{feature} Distributions Train vs. Test')\n    if len(xlim) != 0:\n        ax.set_xlim(xlim)\n    \nfig, axes = plt.subplots(ncols=3, nrows=2, figsize=(20, 12))\nplot_distribution_curves(axes[0, 0], all_segments,'sensor_6_fft_real_median', (-2000, 2000))\nplot_distribution_curves(axes[0, 1], all_segments,'sensor_10_sum', (-5, 5))\nplot_distribution_curves(axes[0, 2], all_segments,'sensor_3_fft_imag_skew', ())\nplot_distribution_curves(axes[1, 0], all_segments,'sensor_2_sum', (-1, 1))\nplot_distribution_curves(axes[1, 1], all_segments,'sensor_1_fft_imag_median', (-1, 1))\nplot_distribution_curves(axes[1, 2], all_segments,'sensor_1_fft_imag_skew', (-1, 1))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Important Adversarial Sensors\nIs there a particular sensor that is responsible for most of the descrimination between train and test. We also want to see the sensors which create the least important features because those sensors might have less differences between train and test. Here it looks like there is pretty even between sensors with **sensor 10** being the only outstanding problematic sensor "},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importances_by_sensor = pd.DataFrame()\nfeature_importances_by_sensor['features'] = columns\nfeature_importances_by_sensor['feature_importances'] = lgbmc.feature_importances_\nfeature_importances_by_sensor['sensor'] = feature_importances_by_sensor['features'].apply(lambda r: r.split('_')[1])\nfeature_importances_by_sensor = feature_importances_by_sensor[feature_importances_by_sensor['sensor'] != 'A0']\nfeature_importances_by_sensor['sensor'] = feature_importances_by_sensor['sensor'].astype(int)\nfeature_importances_by_sensor = feature_importances_by_sensor.groupby(['sensor'], as_index=False).agg({'feature_importances' : 'sum'})\nfig, ax = plt.subplots(figsize = (20, 10))\nsns.set()\nplt.subplot(1, 1, 1);\nsns.barplot(x=\"sensor\", y=\"feature_importances\", orient='v', data=feature_importances_by_sensor);\nplt.title('Feature Importances By Sensor'); ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Important Adversarial Feature Functions\nAre there aggregation functions that were used to create the features. It appears that sum of the fast-fourier transform features are somewhat problematic in this particular instance."},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importances_by_function = pd.DataFrame()\nfeature_importances_by_function['features'] = columns\nfeature_importances_by_function['feature_importances'] = lgbmc.feature_importances_\nfeature_importances_by_function['function'] = feature_importances_by_function['features'].apply(lambda r: '_'.join(r.split('_')[2:]))\nfeature_importances_by_function = feature_importances_by_function.groupby(['function'], as_index=False).agg({'feature_importances' : 'sum'})\nfeature_importances_by_function = feature_importances_by_function.sort_values(['feature_importances'], ascending=False).head(30)\nfig, ax = plt.subplots(figsize = (20, 10))\nsns.set()\nplt.subplot(1, 1, 1);\nsns.barplot(x=\"function\", y=\"feature_importances\", orient='v', data=feature_importances_by_function);\nplt.xticks(rotation=45, ha='right')\nplt.title('Feature Importances By Function'); ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## High confidence test predictions\n\nHere we look at the distribution of probabilities output from the model. We can see that there is a gaussian distribution around `0.4` and a spike at around `0.95`. Since our model is calibrated from using the `binary_logloss` metric we know that if we look at the subset of 95% confident predictions roughly 95/100 of those segments will be true test segments. We should try and determine what make those segments similar, why are they different from the others. **Are these the segments that are the key to the competition?** Which segments are these? What can we do about them? **And what can we do with the train segments to account for these?**"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(20, 12))\nall_probabilities = np.append(\n    cv_scores['fold_0_predict_proba'].to_numpy(),\n    [cv_scores['fold_1_predict_proba'].to_numpy(),\n    cv_scores['fold_2_predict_proba'].to_numpy(),\n    cv_scores['fold_3_predict_proba'].to_numpy(),\n    cv_scores['fold_4_predict_proba'].to_numpy()]\n)\n\nsns.distplot(all_probabilities, 100)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def append_cv_scores(cv_scores, postfix):\n    return np.append(\n    cv_scores[f'fold_0_{postfix}'].to_numpy(),\n    [cv_scores[f'fold_1_{postfix}'].to_numpy(),\n    cv_scores[f'fold_2_{postfix}'].to_numpy(),\n    cv_scores[f'fold_3_{postfix}'].to_numpy(),\n    cv_scores[f'fold_4_{postfix}'].to_numpy()]\n)\n\nall_segment_probabilities = pd.DataFrame()\nall_segment_probabilities['segment_id'] = append_cv_scores(cv_scores, 'segment_id')\nall_segment_probabilities['predict_proba'] = append_cv_scores(cv_scores, 'predict_proba')\nall_segment_probabilities['val'] = append_cv_scores(cv_scores, 'val')\nall_segment_probabilities = all_segment_probabilities.sort_values(['predict_proba'], ascending=False)\nall_segment_probabilities.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## High Confidence test predictions similarities\nSo there are no train segments that fall into the high confidence bucket. It looks like the defining characteristic is that `sensor_10` is missing in those training sets, and `sensor_5`, `sensor_9`, and `sensor_2` is often missing in these training sets as well. Not the most interesting insight, but still a useful observation none the less.\n\n- `sensor_10` missing `100%` of the time\n- `sensor_9` missing `47.7%` of the time\n- `sensor_5` missing `63.3%` of the time\n- `sensor_2` missing `45.3%` of the time"},{"metadata":{"trusted":true},"cell_type":"code","source":"all_segment_probabilities_gt_80 = all_segment_probabilities[all_segment_probabilities['predict_proba'] > .80]\nall_segment_probabilities_gt_80.groupby(['val'], as_index=False).agg({'segment_id' : 'count'})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_sensors_gt_80 = test[test.segment_id.isin(all_segment_probabilities_gt_80.segment_id)].copy()\nfor i in range(1, 11):\n    missing_sensors_gt_80[f'sensor_{i}_missing'] = np.where(missing_sensors_gt_80[f'sensor_{i}_sum'] == 0, 1, 0)\n\nmissing_sensors_gt_80.agg({\n    'sensor_1_missing'  : 'sum',\n    'sensor_2_missing'  : 'sum',\n    'sensor_3_missing'  : 'sum',\n    'sensor_4_missing'  : 'sum',\n    'sensor_5_missing'  : 'sum',\n    'sensor_6_missing'  : 'sum',\n    'sensor_7_missing'  : 'sum',\n    'sensor_8_missing'  : 'sum',\n    'sensor_9_missing'  : 'sum',\n    'sensor_10_missing' : 'sum',\n    'segment_id' : 'count'})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\n\nThere are a set of segments in the test segment that are clearly identifiable. They don't have corresponding train segments that are similar, they might be driving the CV/LB differences, we don't have good proxies in train, we might have to create them. I have included a list of the segment_ids in the file `test_segments_adversarial_validation_prob_gt_80.csv` in the output of this notebook, so that others might explore."},{"metadata":{"trusted":true},"cell_type":"code","source":"all_segment_probabilities_gt_80.to_csv('./test_segments_adversarial_validation_prob_gt_80.csv', header=True, index=False)\nall_segment_probabilities.to_csv('./all_segments_adversarial_validation_prob.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Check the distributions of our previous predictions\nSo let's check what the distribution of predictions look like for the segments that look like train, and ones that don't look like train, and overall. We see that the unidentifiable is quite concentrated at around `27540221.376374934` (the median), with a mean of `26917173.260142826`. These are quite close to the predicted test `time_to_eruption`mean `24073643.301489953` and the train `time_to_eruption` mean `22848906.832769126`."},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/submission_recent.csv')\ntrain = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/train.csv')\npredictions_gt_80 = predictions[predictions.segment_id.isin(all_segment_probabilities_gt_80.segment_id)]\npredictions_lt_80 = predictions[~predictions.segment_id.isin(all_segment_probabilities_gt_80.segment_id)]\n\npredictions_gt_80_sensor10_null = predictions[predictions.segment_id.isin(missing_sensors_gt_80[missing_sensors_gt_80[f'sensor_10_sum'] == 0].segment_id)]\npredictions_gt_80_sensor9_null = predictions[predictions.segment_id.isin(missing_sensors_gt_80[missing_sensors_gt_80[f'sensor_9_sum'] == 0].segment_id)]\npredictions_gt_80_sensor5_null = predictions[predictions.segment_id.isin(missing_sensors_gt_80[missing_sensors_gt_80[f'sensor_5_sum'] == 0].segment_id)]\npredictions_gt_80_sensor2_null = predictions[predictions.segment_id.isin(missing_sensors_gt_80[missing_sensors_gt_80[f'sensor_2_sum'] == 0].segment_id)]\n\ntrain_lt_80_sensor10_null = train[train[f'sensor_10_sum'] == 0]\ntrain_lt_80_sensor9_null = train[train[f'sensor_9_sum'] == 0]\ntrain_lt_80_sensor5_null = train[train[f'sensor_5_sum'] == 0]\ntrain_lt_80_sensor2_null = train[train[f'sensor_2_sum'] == 0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\n\nsns.distplot(predictions_gt_80['time_to_eruption'], 50)\nsns.distplot(predictions['time_to_eruption'], 50)\nsns.distplot(predictions_lt_80['time_to_eruption'], 50)\nsns.distplot(train['time_to_eruption'], 50)\nplt.title('Distribution of predictions for Test segments - Blue: Test identifiable, Orange: Test all, Green: Test unidentifiable, Pink: Train time to eruption')\nplt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# A Guess...\nWhat is happening here... Well I'm guessing it's one of a couple things. Some hyptheses\n\n1. Sensor 10, 9, and 5 are sensitive. They fail intermittently and the researchers are going out to fix them or put them online again\n2. Sensors 10, 9, 5 failed at different points in time and were fixed once or twice.\n3. Sensors 10, 5 failed once for a while, sensors 9 and 2 failed intermittently once or twice\n\nWe can test these hypotheses by looking at the distributions of the predicted `time_to_eruption` of the test versus the predicted `time_to_eruption` of the train. In the below graph, Blue is the unidentifiable test, and orange is train. Looking at those side by side it seems like our theory number 1 and 3 are correct. Sensors 9 and 10 seem to fail intermittently while sensors 5 and 2 fail once for a while.\n\nUsing the below method let's assume our predictions are directionally correct. The distributions do not match however, let's see if we can build a regression model to remap these predictions so that the failure timelines match between test and train."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(20, 12))\n\nsns.distplot(predictions_gt_80_sensor10_null['time_to_eruption'], 50, ax = axes[0,0])\nsns.distplot(train_lt_80_sensor10_null['time_to_eruption'], 50, ax = axes[0,0])\naxes[0,0].set_title('Identifiable Test vs. Train time_to_eruption sensor 10 null')\nsns.distplot(predictions_gt_80_sensor9_null['time_to_eruption'], 50, ax = axes[0,1])\nsns.distplot(train_lt_80_sensor9_null['time_to_eruption'], 50, ax = axes[0,1])\naxes[0,1].set_title('Identifiable Test vs. Train time_to_eruption sensor 9 null')\nsns.distplot(predictions_gt_80_sensor5_null['time_to_eruption'], 50, ax = axes[1,0])\nsns.distplot(train_lt_80_sensor5_null['time_to_eruption'], 50, ax = axes[1,0])\naxes[1,0].set_title('Identifiable Test vs. Train time_to_eruption sensor 5 null')\nsns.distplot(predictions_gt_80_sensor2_null['time_to_eruption'], 50, ax = axes[1,1])\nsns.distplot(train_lt_80_sensor2_null['time_to_eruption'], 50, ax = axes[1,1])\naxes[1,1].set_title('Identifiable Test vs. Train time_to_eruption sensor 2 null')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Simple regression to remap distributions\n\nLet's just reset here. I want to use the `time_to_eruption`s from my best submission `6e6` to do this. We are going to use the distributions from the train datasets to remap the distributions for the identifiable train `time_to_eruption` predictions. We are going to go backwards by the number observations we have in train: `sensor_2`, `sensor_5`, `sensor_9`, `sensor_10`.\n\nWe can use a simple neural network to build these non-linear mappings, we then substitute the new predictions back into the predictions dataset. We can see finally that our distributions of `time_to_eruption` better matches the unidentifiable test dataset and the train dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get the best scoring predictions\npredictions = pd.read_csv('../input/ignv-submission-6e6/submission_best.csv')\n\n# Get train and test again\ntrain = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/train.csv')\ntest = pd.read_csv('../input/ingv-volcanic-eruption-prediction-lgbm-baseline/test.csv')\n\n# Split out identifiable and unidentifiable test segments and their predictions\npredictions_gt_80 = predictions[predictions.segment_id.isin(all_segment_probabilities_gt_80.segment_id)]\npredictions_lt_80 = predictions[~predictions.segment_id.isin(all_segment_probabilities_gt_80.segment_id)]\n\n# Build function to return to us the identifiable test segments by null sensors, and the train segments by null sensors\ndef segments_by_null_sensors(predictions_gt_80, test, train, sensor):\n    predictions_null_sensors = predictions_gt_80[predictions_gt_80.segment_id.isin(test[test[f'{sensor}_sum'] == 0].segment_id)]\n    train_predictions_null_sensors = train[train[f'{sensor}_sum'] == 0]\n    train_predictions_null_sensors = train_predictions_null_sensors[['segment_id', 'time_to_eruption']].copy()\n    return predictions_null_sensors, train_predictions_null_sensors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.neighbors import KernelDensity\n\ndef get_regression_dataset(predictions_null_sensors, train_predictions_null_sensors):\n\n    # Build Kernel density estimators\n    kde_predictions = KernelDensity(kernel='gaussian', bandwidth=0.75) \\\n                                   .fit(predictions_null_sensors['time_to_eruption'].to_numpy()[:, np.newaxis])\n    kde_train = KernelDensity(kernel='gaussian', bandwidth=0.75) \\\n                             .fit(train_predictions_null_sensors['time_to_eruption'].to_numpy()[:, np.newaxis])\n    \n    # Build KDE samples\n    kde_predictions_sample = kde_predictions.sample(100000, 42).flatten()\n    kde_train_sample = kde_train.sample(100000, 42).flatten()\n    \n    quantiles = []\n    kde_predictions_sample_quantile = []\n    kde_train_sample_quantile = []\n    for quantile in np.linspace(0, 1, 1000):\n        quantiles.append(quantile)\n        kde_predictions_sample_quantile.append(np.quantile(kde_predictions_sample, quantile))\n        kde_train_sample_quantile.append(np.quantile(kde_train_sample, quantile))\n    \n    pdf_mapping_df = pd.DataFrame()\n    pdf_mapping_df['quantiles'] = quantiles\n    pdf_mapping_df['predictions_tte'] = kde_predictions_sample_quantile\n    pdf_mapping_df['train_tte'] = kde_train_sample_quantile\n    \n    return pdf_mapping_df\n\npredictions_null_sensor_2, train_predictions_null_sensor_2 = segments_by_null_sensors(predictions_gt_80, test, train, 'sensor_2')\npdf_mapping_df_sensor_2 = get_regression_dataset(predictions_null_sensor_2, train_predictions_null_sensor_2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.scatterplot(x= pdf_mapping_df_sensor_2['predictions_tte'], y = pdf_mapping_df_sensor_2['train_tte'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom collections import OrderedDict\n\ndef training_loop(n_epochs, optimizer, model, loss_fn, t_u_train, t_u_val,\n                  t_c_train, t_c_val):\n    for epoch in range(1, n_epochs + 1):\n        t_p_train = model(t_u_train)\n        loss_train = loss_fn(t_p_train, t_c_train)\n        t_p_val = model(t_u_val)\n        loss_val = loss_fn(t_p_val, t_c_val)\n        optimizer.zero_grad()\n        loss_train.backward()\n        optimizer.step()\n\n        if epoch == 1 or epoch % 10000 == 0:\n            print(f\"Epoch {epoch}, Training loss {loss_train.item():.4f},\"f\" Validation loss {loss_val.item():.4f}\")\n            \n\ndef pdf_mapping_model(pdf_mapping_df, predictions_null_sensor, sensor):\n    pdf_mapping_df_pred_tte_max_scaled = pdf_mapping_df['predictions_tte'].copy().to_numpy() / pdf_mapping_df['predictions_tte'].copy().max()\n    pdf_mapping_df_train_tte_max_scaled = pdf_mapping_df['train_tte'].copy().to_numpy() / pdf_mapping_df['predictions_tte'].copy().max()\n\n    print (f'Building model for {sensor}')\n    # Neural Network to remap predictions based on PDF\n    seq_model = nn.Sequential(\n    OrderedDict([('hidden_linear1', nn.Linear(1, 8)),\n                 ('hidden_activation1', nn.Tanh()),\n                 ('hidden_linear2', nn.Linear(8, 8)),\n                 ('hidden_activation2', nn.Tanh()),\n                 ('output_linear', nn.Linear(8, 1))])\n    )\n    optimizer = optim.SGD(seq_model.parameters(), lr=1e-2)\n    \n    training_loop(\n        n_epochs = 20000,\n        optimizer = optimizer,\n        model = seq_model.double(),\n        loss_fn = nn.MSELoss(),\n        t_u_train = torch.tensor(pdf_mapping_df_pred_tte_max_scaled, dtype=torch.double).unsqueeze(1),\n        t_u_val = torch.tensor(pdf_mapping_df_pred_tte_max_scaled, dtype=torch.double).unsqueeze(1),\n        t_c_train = torch.tensor(pdf_mapping_df_train_tte_max_scaled, dtype=torch.double).unsqueeze(1),\n        t_c_val = torch.tensor(pdf_mapping_df_train_tte_max_scaled, dtype=torch.double).unsqueeze(1)\n    )\n    \n    # Results\n    predictions_null_sensor_max_scaled = predictions_null_sensor['time_to_eruption'].copy().to_numpy() / pdf_mapping_df['predictions_tte'].copy().max()\n    \n    predictions_null_sensor['time_to_eruption_remapped'] = seq_model(\n                    torch.tensor(predictions_null_sensor_max_scaled, dtype=torch.double).unsqueeze(1)\n                    ).detach().numpy() * pdf_mapping_df['predictions_tte'].copy().max()\n    \n    return predictions_null_sensor\n\npredictions_null_sensor_2 = pdf_mapping_model(pdf_mapping_df_sensor_2, predictions_null_sensor_2, 'sensor_2')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(predictions_null_sensor_2['time_to_eruption_remapped'])\nsns.distplot(predictions_null_sensor_2['time_to_eruption'])\nsns.distplot(train_predictions_null_sensor_2['time_to_eruption'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions_null_sensor_5, train_predictions_null_sensor_5 = segments_by_null_sensors(predictions_gt_80, test, train, 'sensor_5')\npdf_mapping_df_sensor_5 = get_regression_dataset(predictions_null_sensor_5, train_predictions_null_sensor_5)\npredictions_null_sensor_5 = pdf_mapping_model(pdf_mapping_df_sensor_5, predictions_null_sensor_5, 'sensor_5')\n\npredictions_null_sensor_9, train_predictions_null_sensor_9 = segments_by_null_sensors(predictions_gt_80, test, train, 'sensor_9')\npdf_mapping_df_sensor_9 = get_regression_dataset(predictions_null_sensor_9, train_predictions_null_sensor_9)\npredictions_null_sensor_9 = pdf_mapping_model(pdf_mapping_df_sensor_9, predictions_null_sensor_9, 'sensor_9')\n\npredictions_null_sensor_10, train_predictions_null_sensor_10 = segments_by_null_sensors(predictions_gt_80, test, train, 'sensor_10')\npdf_mapping_df_sensor_10 = get_regression_dataset(predictions_null_sensor_10, train_predictions_null_sensor_10.reset_index(drop=True).iloc[1:])\npredictions_null_sensor_10 = pdf_mapping_model(pdf_mapping_df_sensor_10, predictions_null_sensor_10, 'sensor_10')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(predictions_null_sensor_10['time_to_eruption_remapped'])\nsns.distplot(predictions_null_sensor_10['time_to_eruption'])\nsns.distplot(train_predictions_null_sensor_10.reset_index(drop=True).iloc[1:]['time_to_eruption'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions_null_sensor_2 = predictions_null_sensor_2.rename(columns={'time_to_eruption_remapped' : 'time_to_eruption_remapped_sensor_2'})\npredictions_null_sensor_5 = predictions_null_sensor_5.rename(columns={'time_to_eruption_remapped' : 'time_to_eruption_remapped_sensor_5'})\npredictions_null_sensor_9 = predictions_null_sensor_9.rename(columns={'time_to_eruption_remapped' : 'time_to_eruption_remapped_sensor_9'})\npredictions_null_sensor_10 = predictions_null_sensor_10.rename(columns={'time_to_eruption_remapped' : 'time_to_eruption_remapped_sensor_10'})\n\npredictions = pd.merge(predictions, predictions_null_sensor_2[['segment_id', 'time_to_eruption_remapped_sensor_2']], \n                       on=['segment_id'], how='left')\npredictions = pd.merge(predictions, predictions_null_sensor_5[['segment_id', 'time_to_eruption_remapped_sensor_5']], \n                       on=['segment_id'], how='left')\npredictions = pd.merge(predictions, predictions_null_sensor_9[['segment_id', 'time_to_eruption_remapped_sensor_9']], \n                       on=['segment_id'], how='left')\npredictions = pd.merge(predictions, predictions_null_sensor_10[['segment_id', 'time_to_eruption_remapped_sensor_10']], \n                       on=['segment_id'], how='left')\npredictions.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# predictions['substitute_time_to_eruption'] = np.where(predictions['time_to_eruption_remapped_sensor_10'].notnull(),\n#                                                       predictions['time_to_eruption_remapped_sensor_10'],\n#                                                       predictions['time_to_eruption'])\n\npredictions['new_time_to_eruption'] = np.where(predictions['time_to_eruption_remapped_sensor_10'].isnull(), predictions['time_to_eruption'], \n                                               predictions['time_to_eruption_remapped_sensor_10'])\n\npredictions.to_csv('./predictions_all_mappings.csv', header=True, index=False)\npredictions[['segment_id', 'new_time_to_eruption']].rename(columns={'new_time_to_eruption' : 'time_to_eruption'}) \\\n            .to_csv('./new_predictions.csv', header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\n\nsns.distplot(predictions['new_time_to_eruption'], 50)\nsns.distplot(predictions['time_to_eruption'], 50)\nsns.distplot(predictions_lt_80['time_to_eruption'], 50)\nsns.distplot(train['time_to_eruption'], 50)\nplt.title('Distribution of predictions for Test segments - Blue: Test identifiable, Orange: Test all, Green: Test unidentifiable, Pink: Train time to eruption')\nplt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Final Results\nNote that our preidction distribution above is much flatter and closer to the `time_to_eruption` distribution for the unidentifiable train. Using this distribution remapping for the indentifiable train parameters I reduce the score from `6e6` -> `5.8e6`. Again, please note I am only using the features from my original notebook."},{"metadata":{},"cell_type":"markdown","source":"# Euclidean Distances between Train Segments\n\nIt is clear that there are overlapping segments that are not similar. Let's assume that there are then multiple eruptions, if this is true, we should be able to find segments that came from the same eruption and segments that do not come from the same eruption. We can try and find these segments by first finding segments that overlap, these cannot be from the same eruption. Second we calculate the 10 nearest (in time to eruption) of the non-overlapping segments. From there we calculate the euclidean distance between those segments and finally we construct a wieghted graph of the train segments, where the weight is the euclidean distance between two segments.\n\nThe graph is saved in a pickle file in the output of this notebook."},{"metadata":{"trusted":true},"cell_type":"code","source":"train['start_time'] = train['time_to_eruption'] + 60000","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlapping_segments = {}\nten_closest_non_overlapping_segments = {}\n\ntrain = train.sort_values(['time_to_eruption']).reset_index(drop=True)\n\ndef go_backwards(\n    i, train, \n    segments_overlapping_segments):\n    \n    j = i-1\n    while(j >= 0):\n        if (abs(train.iloc[j]['time_to_eruption'] - train.iloc[i]['time_to_eruption']) <= 60001):\n            segments_overlapping_segments.append(train.iloc[j]['segment_id'])\n            j = j-1\n        else:\n            break\n    return segments_overlapping_segments\n\n\ndef go_forwards(\n    i, train, \n    segments_overlapping_segments):\n    \n    j = i+1\n    while(j < train.shape[0]):\n        if (abs(train.iloc[j]['time_to_eruption'] - train.iloc[i]['time_to_eruption']) <= 60001):\n            segments_overlapping_segments.append(train.iloc[j]['segment_id'])\n            j = j+1\n        else:\n            break\n    return segments_overlapping_segments\n\n\nfor i, row in train.iterrows():\n    segments_overlapping_segments = []\n    \n    segments_overlapping_segments = go_backwards(\n                                        i, train, \n                                        segments_overlapping_segments\n                                    )\n    segments_overlapping_segments = go_forwards(\n                                        i, train, \n                                        segments_overlapping_segments\n                                    )\n    ten_closest = train[(~train.segment_id.isin(segments_overlapping_segments)) & (train.segment_id != train.iloc[i].segment_id)].copy()\n    ten_closest['closeness_by_tte'] = np.abs((ten_closest['time_to_eruption'] - train.iloc[i]['time_to_eruption']).to_numpy())\n    segments_ten_closest_non_overlapping = list(ten_closest.sort_values(['closeness_by_tte']).head(10)['segment_id'])\n    \n    overlapping_segments[train.iloc[i].segment_id] = segments_overlapping_segments\n    ten_closest_non_overlapping_segments[train.iloc[i].segment_id] = segments_ten_closest_non_overlapping","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics.pairwise import euclidean_distances\nimport networkx as nx\n\n\ntrain_segment_similarity = pd.DataFrame()\ntrain_segment_similarity['segment_id'] = train.segment_id\ntrain_segment_similarity = train_segment_similarity.set_index('segment_id')\n#train_segment_similarity.columns = train.index\n\nk = 0\ntrain_segment_similarity_graph = nx.Graph()\nfor i, row in train_segment_similarity.iterrows():\n    for e, j in enumerate(ten_closest_non_overlapping_segments[i]):\n        x = train.set_index('segment_id').loc[i].drop(['time_to_eruption', 'start_time']).to_numpy()\n        y = train.set_index('segment_id').loc[j].drop(['time_to_eruption', 'start_time']).to_numpy()\n        distance = euclidean_distances(x.reshape(1, -1),y.reshape(1, -1)).flatten()[0]\n        train_segment_similarity.loc[i, e + 1] = distance\n        \n        # Build Graph\n        if (i not in train_segment_similarity_graph.nodes):\n            train_segment_similarity_graph.add_node(i)\n        if (j not in train_segment_similarity_graph.nodes):\n            train_segment_similarity_graph.add_node(j)\n            \n        train_segment_similarity_graph.add_weighted_edges_from(\n            [(i, j, distance)]\n        )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_segment_similarity.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 15))\n\nnx.draw(train_segment_similarity_graph, with_labels=False) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nx.write_gpickle(train_segment_similarity_graph,\"train_segment_similarity_graph.gpickle\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_segment_similarity_graph","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Final Thoughts\n\n- If there are multiple eruptions the holdout set may very contain an eruption sequence we have not seen in the train data, one were sensor 10 failed for a significant amount of time. This might then also explain the overfitting we are seeing. We overfit on individual eruption sequences that are hidden in our train data. Once we are presented with a new eruption sequence in test our scores go from `2.6e6` -> `6e6`. My guess is that there are about 6-7 individual eruptions in train\n\n- If you can identify the segments that belong to the same eruption sequence you should CV by holding out an eruption sequence or, even better, time series cross-validate over each subsequent eruption sequence. However, the latter might be difficult because you would then have to determine ordering of the identified eruption sequences\n\n- A good CV model should have a local CV close to the LB CV and it should be higher than `2.6e6` but I imagine that the LB scores will not directly match because the test set does have those segments that are quite unusual. It will be interesting to see how people manage to deal with those."}],"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}