{"cells":[{"metadata":{},"cell_type":"markdown","source":"## INGV Challenge - Volcanic Eruption Prediction\nhttps://www.kaggle.com/c/predict-volcanic-eruptions-ingv-oe  \n\n## Summary\nIn part 1, we will read data and do some basic EDA.  \nIn the part 2 of this notebook, we will in generate almost 900 hundreds features in less than 35 minutes using multiprocessing.  \nFeatures will be saved in a Kaggle dataset.\n\nIn part 3, we will fit a LGBM model, like in several notebooks. And we will see how to boost our CV MAE Score by removing many features.\n\nBut we will see that CV MAE Score is much better than the Public LB Score : test sets and train sets are really different as we will see in part 4. We have a big problem of overfitting.\n\nThanks to  \nhttps://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft  \nhttps://www.kaggle.com/jesperdramsch/introduction-to-volcanology-seismograms-and-lgbm   \nhttps://www.kaggle.com/kyakovlev/m5-three-shades-of-dark-darker-magic  \nhttps://www.kaggle.com/gpreda/lanl-earthquake-eda-and-prediction  \nhttps://www.kaggle.com/tunguz/ms-malware-adversarial-validation  \nand many more..."},{"metadata":{},"cell_type":"markdown","source":"## Constants and Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"__seed = 666\n\n__create_dataframes = True  # To create Train and Test datasets. If False, data will we read in \"../input/ingvchallengefeatures/\"\n__dataframe_size = 5000     # to debug : use small value like 10 or 100, or use 5000 to create complete dataframes\n\n__path_to_my_data = \"../input/ingvchallengefeatures/\"  # Path to read Train and Test dataframes, if we didn't create them\n\n__n_folds = 5","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np \nnp.random.seed(__seed)\n\nimport scipy\nimport scipy.signal\n\nimport pandas as pd\npd.set_option('max_columns', 100)\npd.set_option('max_rows', 200)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport gc, pickle, os, itertools, datetime\nfrom functools import partial\n\nimport psutil\n__n_cores = psutil.cpu_count()     # Available CPU cores\nfrom multiprocessing import Pool   # Multiprocess Runs\n\nfrom sklearn.model_selection import KFold, train_test_split\nfrom sklearn.metrics import mean_absolute_error, r2_score\n\nimport lightgbm as lgbm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 1 : Train and Sample_submission files"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"y_train = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train.csv\"\n                      , index_col=\"segment_id\"\n                      , dtype={\"time_to_eruption\":np.int32}\n                     , squeeze=True)   # to return a Series and not a DataFrame\nprint(\"Length of time_to_eruption serie : {}\\n\\nThe five first lines :\".format(len(y_train)))\nprint(y_train.head().map('{:,.0f}'.format))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Time to Eruption in Train\\n-------------------------\\nMin : {:,}\\nMedian : {:,}\\nMax : {:,}\".format(\n    y_train.min(), y_train.median(), y_train.max()))\n\n# Convert 'time_to_eruption'to hours:minutes:seconds (Just for reference)\n# Thanks to Amanooo in\n#    https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft\n_temp = y_train.apply(lambda x:datetime.timedelta(milliseconds = x))\nprint(\"\\nTime to Eruption in Train\\n-------------------------\\nMin : {}\\nMedian : {}\\nMax : {}\".format(\n    _temp.min(), _temp.median(), _temp.max()))\n\ndel _temp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# matplotlib histogram\nplt.hist(y_train, color = 'blue', edgecolor = 'black', bins=150, alpha=0.5)\n\n# Add labels\nplt.title('Distribution of time to eruption in Train')\nplt.xlabel('Time')\nplt.ylabel('Frequency');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Number of sequences\\n-------------------\\nIn Train : {} files\\nIn Test : {} files\".format(\n    len(os.listdir('../input/predict-volcanic-eruptions-ingv-oe/train/'))\n    , len(os.listdir('../input/predict-volcanic-eruptions-ingv-oe/test/'))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv\")\ny_pred = sample_submission.set_index(\"segment_id\")[\"time_to_eruption\"].astype(np.int32)\ny_pred.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Can we do something with Segment ID ? No, it's just a index."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Quantile of segment Id in Train : \\n{}\\n\".format(y_train.index.to_series().quantile([0, .1, .5, .9, 1]).map('{:,.0f}'.format)))\nprint(\"Quantile of segment Id in Test : \\n{}\\n\".format(y_pred.index.to_series().quantile([0, .1, .5, .9, 1]).map('{:,.0f}'.format)))\n\nplt.boxplot([y_train.index.to_series(), y_pred.index.to_series()])\nplt.ylabel(\"Segment ID\")\nplt.xticks([1, 2], ['Train segment ID', 'Test segment ID']);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_file_submission(pred, filename = \"submission.csv\", verbose = False):\n    \n    sample_submission[\"time_to_eruption\"] = pred\n    sample_submission.loc[sample_submission[\"time_to_eruption\"]<0, \"time_to_eruption\"]=0\n    \n    sample_submission.to_csv(filename, index=False)\n    \n    if verbose:\n        print(sample_submission.head())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Show some sequences"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Thanks to jesperdramsch in\n#    https://www.kaggle.com/jesperdramsch/introduction-to-volcanology-seismograms-and-lgbm\ndef plot_sequence(segment_id=\"1000015382\"):\n    sequence = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/\"+segment_id+\".csv\", dtype=\"Int16\")\n    sequence.fillna(0).plot(subplots=True, figsize=(25, 10))\n    plt.tight_layout()\n    plt.show()\n    \nprint(y_train.sort_values().head())\n\n# Some sequences with low value of time to eruption\nplot_sequence(\"1658693785\")\n#plot_sequence(\"1957235969\")\n#plot_sequence(\"442994108\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Some sequences with high value of time to eruption\nprint(y_train.sort_values().tail())\nplot_sequence(\"1162128945\")\n#plot_sequence(\"1131527270\")\n#plot_sequence(\"356854390\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 2 : Features engineering\nAs fast as possible with Pandas, pandas.agg and multiprocessing Pool.\n\nAbout multiprocessing, look at :\n* https://www.kaggle.com/kyakovlev/m5-three-shades-of-dark-darker-magic  \n* https://sites.google.com/site/python3tutorial/multiprocessing_map/multiprocessing_partial_function_multiple_arguments \n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"## Multiprocess Runs\n# Thanks to Konstantin Yakovlev in :\n#    https://www.kaggle.com/kyakovlev/m5-three-shades-of-dark-darker-magic\ndef df_parallelize_run(func, t_split):\n    \n    num_cores = np.min([__n_cores, len(t_split)])\n    pool = Pool(num_cores)\n    df = pd.concat(pool.map(func, t_split), axis=0)\n    pool.close()\n    pool.join()\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Put all statistics for many columns on one unique row\n## and put name to columns\ndef add_stats_on_the_row(df_stats, df_result = None):\n    \n    # Name of columns\n    input_list = list(itertools.product(*[list(df_stats.columns), list(df_stats.index.values)]))\n    columns = [f'{e[0]}_{str(e[1])}' for e in input_list]\n    \n    # Create dataframe with statistics, on a unique row\n    res = pd.DataFrame(df_stats.values.reshape(1, df_stats.shape[0] * df_stats.shape[1], order=\"F\"), columns = columns)\n    \n    # Add those statistics to a pre-existing data frame (or not)\n    if df_result is None:\n        return res\n    else:\n        return pd.concat([df_result, res], axis = 1)\n\n# An example :\nprint(\"Little example, how I put statistics on a unique row :\")\ndf_ex = pd.DataFrame({'a':[100, 1, 2], 'b':[10, 20, 3], 'c':[20, 60, 7]})\n\n# Use the function and pandas.agg\ndf_stats = add_stats_on_the_row(df_ex[[\"a\", \"c\"]].agg([\"min\", \"mean\", \"max\"]))\n\nprint(\"Let's create a litte data frame \\n {}\\n\\nAnd use this function to compute statistics : \\n{}\".\n      format(df_ex, df_stats))\n\ndel df_ex, df_stats     ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Create lists of features names\ndef lfeat_with_suff(suffix=\"\"):\n    return [f\"sensor_{i+1}{suffix}\" for i in range (10)]\n\nlfeat = lfeat_with_suff()             # = [\"sensor_1\", \"sensor_2\", ..., \"sensor_10\"]\nlfeat_abs = lfeat_with_suff(\"_abs\")   # = [\"sensor_1_abs\", \"sensor_2_abs\", ..., \"sensor_10_abs\"]\nlfeat_fft_real = lfeat_with_suff(\"_fft_real\")   \nlfeat_fft_imag = lfeat_with_suff(\"_fft_imag\")   \n\ninput_list = [range(10), [\"min\", \"max\"]]\nlist_result = list(itertools.product(*input_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/amanooo/ingv-volcanic-basic-solution-stft\n# Thanks to Amanooo\n\n# STFT Specifications\nfs = 100                # sampling frequency \n#N = len(segment_df)     # data size\nN = 60001               # data size\nn = 256                 # FFT segment size\nmax_f = 20              # ～20Hz\n\ndelta_f = fs / n        # 0.39Hz\ndelta_t = n / fs / 2    # 1.28s\n\n\ndef STFT_Features(segment_df, segment_id):\n    \n    segment = [segment_id]\n    \n    for sensor in segment_df.columns:\n        x = segment_df[sensor][:N]\n        if x.isna().sum() > 1000:     ##########\n            segment += ([np.NaN] * 10)\n            continue\n        f, t, Z = scipy.signal.stft(x.fillna(0), fs = fs, window = 'hann', nperseg = n)\n        f = f[:round(max_f/delta_f)+1]\n        Z = np.abs(Z[:round(max_f/delta_f)+1]).T    # ～max_f, row:time,col:freq\n\n        th = Z.mean() * 1     ##########\n        Z_pow = Z.copy()\n        Z_pow[Z < th] = 0\n        Z_num = Z_pow.copy()\n        Z_num[Z >= th] = 1\n\n        Z_pow_sum = Z_pow.sum(axis = 0)\n        Z_num_sum = Z_num.sum(axis = 0)\n\n        A_pow = Z_pow_sum[round(10/delta_f):].sum()\n        A_num = Z_num_sum[round(10/delta_f):].sum()\n        BH_pow = Z_pow_sum[round(5/delta_f):round(8/delta_f)].sum()\n        BH_num = Z_num_sum[round(5/delta_f):round(8/delta_f)].sum()\n        BL_pow = Z_pow_sum[round(1.5/delta_f):round(2.5/delta_f)].sum()\n        BL_num = Z_num_sum[round(1.5/delta_f):round(2.5/delta_f)].sum()\n        C_pow = Z_pow_sum[round(0.6/delta_f):round(1.2/delta_f)].sum()\n        C_num = Z_num_sum[round(0.6/delta_f):round(1.2/delta_f)].sum()\n        D_pow = Z_pow_sum[round(2/delta_f):round(4/delta_f)].sum()\n        D_num = Z_num_sum[round(2/delta_f):round(4/delta_f)].sum()\n        segment += [A_pow, A_num, BH_pow, BH_num, BL_pow, BL_num, C_pow, C_num, D_pow, D_num]\n\n    features = [f\"sensor_{i+1}_{f}\" for i in range(10) for f in [\"A_pow\", \"A_num\", \"BH_pow\", \"BH_num\", \"BL_pow\", \"BL_num\"\n                                                    , \"C_pow\", \"C_num\", \"D_pow\", \"D_num\"]]\n    df = pd.DataFrame(np.array([segment]), columns = [\"segment_id\"] + features)\n    df[\"segment_id\"] = df[\"segment_id\"].astype(np.int32)\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Thanks to\n# https://www.kaggle.com/isaienkov/ingv-volcanic-eruption-prediction-eda-modeling\n# https://www.kaggle.com/ajcostarino/ingv-volcanic-eruption-prediction-lgbm-baseline\n# https://www.kaggle.com/gpreda/lanl-earthquake-eda-and-prediction\ndef features_engineering_one_file(onefile, source):\n    \n    # Read a file\n    sequence = pd.read_csv(f'/kaggle/input/predict-volcanic-eruptions-ingv-oe/{source}/{onefile}.csv')\n    \n    res = STFT_Features(sequence[lfeat], onefile)\n    \n    # Count NAN values for each sensor, and replace them by 0\n    res = pd.concat([res, pd.DataFrame(sequence[lfeat].isna().sum().values.reshape(1, 10, order=\"F\")\n                   , columns=[f\"sensor_{i+1}_nan\" for i in range(10)])], axis=1)\n    sequence.fillna(0, inplace = True)\n    \n    # Absolute value of each sensor\n    sequence[lfeat_abs] = sequence[lfeat].abs()\n    \n    # Fast Fourier Transformation for all sensors\n    ft = np.fft.fft(sequence[lfeat], axis=0)\n    sequence[lfeat_fft_real] = np.real(ft)\n    sequence[lfeat_fft_imag] = np.imag(ft)\n    \n    # Basic statistic on each sensor\n    lfun = [\"mean\", \"std\", \"min\", \"max\", \"mad\", \"skew\", \"kurtosis\"]\n    lquantiles = [.01, .05, .1, .25, .5, .75, .9, .95, .99]\n    res = add_stats_on_the_row(sequence[lfeat + lfeat_abs + lfeat_fft_real + lfeat_fft_imag].agg(lfun), res)\n    res = add_stats_on_the_row(sequence[lfeat + lfeat_abs + lfeat_fft_real + lfeat_fft_imag].quantile(lquantiles), res)\n            \n    # End\n    \n    return res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndef create_data_for_model(df, source=\"train\"):\n    \n    # Create a new function which need only one parameter : the file\n    func_partial_fe = partial(features_engineering_one_file, source = source)\n    \n    # Read all sequences in parallel for feature engineering\n    df_set = list()\n    df_set.append(df_parallelize_run(func_partial_fe, df.index[:__dataframe_size]))\n            \n    # Transform to Pandas DataFrame\n    df_set = pd.concat(df_set)\n    df_set.reset_index(inplace = True, drop=True)\n    df_set.set_index(\"segment_id\", inplace=True)\n    \n    return df_set\n    \n# Create dataframes or read them in my Data\nif __create_dataframes: \n    \n    print(\"Train dataset creation...\")\n    train_set = create_data_for_model(df = y_train)\n    \n    if __dataframe_size > 100:\n    \n        print(\"Test dataset creation...\")\n        test_set = create_data_for_model(df = y_pred, source = \"test\")\n    \nelse:\n    \n    train_set = pickle.load(open(__path_to_my_data + \"ingv_train_set.pkl\", \"rb\"))\n    test_set = pickle.load(open(__path_to_my_data + \"ingv_test_set.pkl\", \"rb\"))\n    \n# Write train set and test set\npickle.dump(train_set, open( \"ingv_train_set.pkl\", \"wb\" ) ) \nif __dataframe_size > 100:\n    pickle.dump(test_set, open( \"ingv_test_set.pkl\", \"wb\" ) ) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a=train_set.sum(axis=0, skipna=False)\nprint(\"Features with NAN values : \\n{}\".format(list(a[a.isna()].index)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Some other feature engineering\ndef other_fe(df_e):\n    \n    df = df_e.copy()\n    \n    # Inter quantile ratio\n    lfeat = ['abs_0.5', 'abs_0.75', 'abs_0.95', 'abs_0.99', 'abs_max']\n    lfeat_new1 = [\"q75_o_med\", \"q95_o_q75\", \"q99_o_q95\", \"max_o_q99\"]\n    lfeat_new2 = [\"q95_o_med\", \"q99_o_med\", \"max_o_med\"]\n    \n    for i in range(10):\n        for feat1, num, denom in zip(lfeat_new1, lfeat[1:], lfeat[:-1]):\n            df[f'sensor_{i+1}_{feat1}'] = df[f'sensor_{i+1}_{num}'] / df[f'sensor_{i+1}_{denom}']\n            df.loc[ df[f'sensor_{i+1}_{denom}'] == 0, f'sensor_{i+1}_{feat1}'] = 0\n\n        for feat2, num in zip(lfeat_new2, lfeat[2:]):\n            df[f'sensor_{i+1}_{feat2}'] = df[f'sensor_{i+1}_{num}'] / df[f'sensor_{i+1}_abs_0.5']\n            df.loc[ df[f'sensor_{i+1}_abs_0.5'] == 0, f'sensor_{i+1}_{feat2}'] = 0\n            \n        df[f'sensor_{i+1}_q3_q1_abs'] = df[f'sensor_{i+1}_abs_0.75'] - df[f'sensor_{i+1}_abs_0.25']\n        df[f'sensor_{i+1}_q3_q1'] = df[f'sensor_{i+1}_0.75'] - df[f'sensor_{i+1}_0.25']\n        df[f'sensor_{i+1}_d9_d1_abs'] = df[f'sensor_{i+1}_abs_0.9'] - df[f'sensor_{i+1}_abs_0.1']\n        df[f'sensor_{i+1}_d9_d1'] = df[f'sensor_{i+1}_0.9'] - df[f'sensor_{i+1}_0.1']\n        \n    df[[f\"sensor_{i+1}_empty\" for i in range(10)]]=0\n    for i in range(10):\n        df.loc[df[f\"sensor_{i+1}_abs_max\"]==0, f\"sensor_{i+1}_empty\"]=1\n        \n    df.fillna(0, inplace=True)\n        \n    return df\n\n\n# Some others features\ntrain_set = other_fe(train_set)\nif __dataframe_size > 100:\n    test_set = other_fe(test_set)\ntrain_set[\"time_to_eruption\"] = y_train\ntrain_set.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Do have some features with only one value ? They are useless."},{"metadata":{"trusted":true},"cell_type":"code","source":"a = train_set.nunique()\nto_del = list(a[a<2].index)\nprint(\"Features deleted : \\n{}\".format(to_del))\ntrain_set.drop(to_del, axis=1, inplace=True)\ntest_set.drop(to_del, axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 3 : LGBM"},{"metadata":{"trusted":true},"cell_type":"code","source":"def fit_and_predict_with_lgbm(df_train, params\n            , seed=__seed, features=None, X_test=None, check_feat_importance = True, n_folds=__n_folds\n            , early_stopping_rounds = 100):\n    \n    folds = KFold(n_splits = n_folds, shuffle = True, random_state = seed)\n    list_mae, list_r2 = [], []\n    predictions = None\n    feature_importance_df = pd.DataFrame()\n    feat_imp={}\n\n    params[\"random_state\"] = seed + 1\n    params[\"bagging_seed\"] = seed - 1\n    \n    # Keep only a list of features\n    if not features is None:\n        lfeat = features\n    else:\n        lfeat = [f for f in df_train.columns if f not in [\"time_to_eruption\"]]\n    print(\"Train DF shape : {} rows and {} features\".format(df_train.shape[0], len(lfeat)))\n    \n    if not X_test is None:\n        X_test = X_test[lfeat]\n        print(\"Test DF shape  : {} rows and {} features\".format(X_test.shape[0], X_test.shape[1]))\n        predictions = np.zeros(len(X_test))\n\n    for fold_, (trn_idx, val_idx) in enumerate(folds.split(df_train)):\n    \n        print(\"Fold n°{}\".format(fold_+1))\n        \n        X_tr, X_val = df_train[lfeat].iloc[trn_idx], df_train[lfeat].iloc[val_idx]\n        y_tr, y_val = df_train['time_to_eruption'].iloc[trn_idx], df_train['time_to_eruption'].iloc[val_idx]\n    \n        model = lgbm.LGBMRegressor(**params, n_estimators = 3000, n_jobs = -1)\n        model.fit(X_tr, y_tr, \n              eval_set=[(X_tr, y_tr), (X_val, y_val)], eval_metric='mae',\n              verbose=1000, early_stopping_rounds = early_stopping_rounds)\n    \n        y_val_pred = model.predict(X_val, num_iteration=model.best_iteration_)\n        _mae = mean_absolute_error(y_val_pred, y_val) ; list_mae.append(_mae)\n        _r2 = r2_score(y_val_pred, y_val) ; list_r2.append(_r2)\n        print(\"Scores on valid set for fold {} : MAE {:,.0f} & R2 {:.3f}\".format(fold_+1, _mae, _r2))\n        \n        # Predictions\n        if not X_test is None:\n            predictions += model.predict(X_test[lfeat], num_iteration=model.best_iteration_) / folds.n_splits\n    \n        # Feature importance\n        if check_feat_importance:\n            \n            # to do a nice and useless plot\n            fold_importance_df = pd.DataFrame()\n            fold_importance_df[\"Feature\"] = X_tr.columns\n            fold_importance_df[\"importance\"] = model.feature_importances_[:len(X_tr.columns)]\n            fold_importance_df[\"fold\"] = fold_ + 1\n            feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)\n        \n            # Evaluate each feature importance by random permutation of the feature values\n            # More usefull\n            # https://eli5.readthedocs.io/en/latest/blackbox/permutation_importance.html\n            # https://www.kaggle.com/kyakovlev/m5-three-shades-of-dark-darker-magic\n            for feat in X_tr.columns:\n                if fold_ == 0:\n                    feat_imp[feat] = []\n                temp_df = X_val.copy()\n                temp_df[feat] = np.random.permutation(temp_df[feat])\n                y_temp = model.predict(temp_df, num_iteration=model.best_iteration_)\n                feat_imp[feat].append(mean_absolute_error(y_temp, y_val) - _mae)\n        \n    \n    # Mean MAE\n    _mae = np.array(list_mae).mean()\n    _r2 = np.array(list_r2).mean()\n    print(\"Mean scores on valid sets : MAE {:,.0f} & R2 {:.3f}\".format(_mae, _r2))\n    \n    # Feature importance\n    #   https://eli5.readthedocs.io/en/latest/blackbox/permutation_importance.html :\n    #   \"eli5 provides a way to compute feature importances for any black-box estimator by measuring how score \n    #   decreases when a feature is not available; the method is also known as “permutation importance” \n    #   or “Mean Decrease Accuracy (MDA)”.\n    # Here I count folds where feature permutation produce a score decrease\n    # When a score in a fold decreases for a feature, we consider this feature can be removed.\n    feat_imp_df = None\n    if check_feat_importance:\n        \n        feat_imp_df = pd.DataFrame(feat_imp).transpose()\n\n        lfeat = [i for i in range(__n_folds)]\n        lfeat_sign = [f\"sign_{i}\" for i in range(n_folds)]\n\n        feat_imp_df[lfeat_sign] = 0\n        for mae_loss, mae_loss_sign in zip(lfeat, lfeat_sign):\n            feat_imp_df.loc[feat_imp_df[mae_loss]>0, mae_loss_sign] = 1\n    \n        feat_imp_df[\"sum_sign\"] = feat_imp_df[lfeat_sign].sum(axis=1)\n    \n    return {\"pred\":predictions, \"feat_imp\":feat_imp_df, \"list_mae\":list_mae, \"mean_mae_valid\":_mae\n           , \"list_r2\":list_r2, \"mean_r2_valid\":_r2, \"plot_feat_imp\":feature_importance_df}\n\n\nparams = {\n    'num_leaves': 28,          # small value to avoid overfitting\n    'min_data_in_leaf': 10, \n#    'max_depth': 5,\n    'learning_rate': 0.15,\n    'max_bins': 50,            # small value to avoid overfitting\n    \"feature_fraction\": 0.5,\n    \"bagging_freq\": 1,\n    \"bagging_fraction\": 0.8,\n    \"lambda_l1\": 0.1,\n    \"boosting\": \"gbdt\",\n    'objective':'regression',\n    \"metric\": 'mae',\n    \"verbosity\": -1,\n    \"nthread\": -1,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nresult1 = fit_and_predict_with_lgbm(train_set, params)\n\n# Now we will keep features usefull in all folds\nfeat_imp_df = result1[\"feat_imp\"]\nfeats = list(feat_imp_df.loc[feat_imp_df[\"sum_sign\"] == __n_folds].index)\nprint(\"\\nThere are {} features with no loss of MAE on the {} folds.\\n\".format(len(feats), __n_folds))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nresult2 = fit_and_predict_with_lgbm(train_set, params, seed=__seed+10, features=feats, X_test = test_set)\n\nprint(\"\\n\\nResults summary :\")\nprint(\"CV n°1 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result1[\"mean_mae_valid\"], result1[\"mean_r2_valid\"]))\nprint(\"CV n°2 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result2[\"mean_mae_valid\"], result2[\"mean_r2_valid\"]))\n\nfeat_imp_df = result2[\"feat_imp\"]\nfeats = list(feat_imp_df.loc[feat_imp_df[\"sum_sign\"] == __n_folds].index)\nprint(\"\\nThere are {} features with no loss of MAE on the {} folds.\\n\".format(len(feats), __n_folds))\n\nmake_file_submission(result2[\"pred\"], \"submission2.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nresult3 = fit_and_predict_with_lgbm(train_set, params, seed=__seed+69, features=feats, X_test = test_set)\n\nprint(\"\\n\\nResults summary :\")\nprint(\"CV n°1 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result1[\"mean_mae_valid\"], result1[\"mean_r2_valid\"]))\nprint(\"CV n°2 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result2[\"mean_mae_valid\"], result2[\"mean_r2_valid\"]))\nprint(\"CV n°3 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result3[\"mean_mae_valid\"], result3[\"mean_r2_valid\"]))\n\nfeat_imp_df = result3[\"feat_imp\"]\nfeats = list(feat_imp_df.loc[feat_imp_df[\"sum_sign\"] == __n_folds].index)\nprint(\"\\nThere are {} features with no loss of MAE on the {} folds.\\n\".format(len(feats), __n_folds))\n\nmake_file_submission(result3[\"pred\"], \"submission3.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nresult4 = fit_and_predict_with_lgbm(train_set, params, seed=__seed-41, features=feats, X_test = test_set)\n\nprint(\"\\n\\nResults summary :\")\nprint(\"CV n°1 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result1[\"mean_mae_valid\"], result1[\"mean_r2_valid\"]))\nprint(\"CV n°2 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result2[\"mean_mae_valid\"], result2[\"mean_r2_valid\"]))\nprint(\"CV n°3 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result3[\"mean_mae_valid\"], result3[\"mean_r2_valid\"]))\nprint(\"CV n°4 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result4[\"mean_mae_valid\"], result4[\"mean_r2_valid\"]))\n\nfeat_imp_df = result4[\"feat_imp\"]\nfeats = list(feat_imp_df.loc[feat_imp_df[\"sum_sign\"] == __n_folds].index)\nprint(\"\\nThere are {} features with no loss of MAE on the {} folds.\\n\".format(len(feats), __n_folds))\n\nmake_file_submission(result4[\"pred\"], \"submission4.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nresult5 = fit_and_predict_with_lgbm(train_set, params, seed=__seed-41, features=feats, X_test = test_set)\n\nprint(\"\\n\\nResults summary :\")\nprint(\"CV n°1 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result1[\"mean_mae_valid\"], result1[\"mean_r2_valid\"]))\nprint(\"CV n°2 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result2[\"mean_mae_valid\"], result2[\"mean_r2_valid\"]))\nprint(\"CV n°3 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result3[\"mean_mae_valid\"], result3[\"mean_r2_valid\"]))\nprint(\"CV n°4 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result4[\"mean_mae_valid\"], result4[\"mean_r2_valid\"]))\nprint(\"CV n°4 results  - Mean MAE : {:,.0f}. Mean R2 : {:.3f}\".format(result5[\"mean_mae_valid\"], result5[\"mean_r2_valid\"]))\n\nfeat_imp_df = result5[\"feat_imp\"]\nfeats = list(feat_imp_df.loc[feat_imp_df[\"sum_sign\"] == __n_folds].index)\nprint(\"\\nThere are {} features with no loss of MAE on the {} folds.\\n\".format(len(feats), __n_folds))\n\nmake_file_submission(result5[\"pred\"], \"submission5.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So I submit n°3."},{"metadata":{},"cell_type":"markdown","source":"## About features importances "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Features which could be removed...\nfeat_imp_df[feat_imp_df[\"sum_sign\"] != __n_folds]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importance_df = result5[\"plot_feat_imp\"]\n\ncols = (feature_importance_df[[\"Feature\", \"importance\"]]\n        .groupby(\"Feature\")\n        .mean()\n        .sort_values(by=\"importance\", ascending=False)[:100].index)\n\nbest_features = feature_importance_df.loc[feature_importance_df.Feature.isin(cols)]\n\nplt.figure(figsize=(14,20))\nsns.barplot(x=\"importance\",\n            y=\"Feature\",\n            data=best_features.sort_values(by=\"importance\",\n                                           ascending=False))\nplt.title('LightGBM Features (avg over folds) - The 100th most important features')\nplt.tight_layout();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Part 4 : Adversarial Validation\nTrain and test are really different.  \nThnaks to Bojan Tunguz  \nhttps://www.kaggle.com/tunguz/ms-malware-adversarial-validation  \nor in https://www.kaggle.com/tunguz/ms-malware-adversarial-validation\nand in so many others challenges...\n\nThanks to https://www.kaggle.com/ajcostarino/ignv-adversarial-validation-cv-lb-differences\n\nI think my MAE will grow on the Private LB..."},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndf = pd.concat([train_set.drop(\"time_to_eruption\", axis=1), test_set], axis=0)\nts = train_set.shape\nli = df.iloc[:ts[0]].index.values\n\ndf[\"is_test\"] = 1\ndf.loc[li, \"is_test\"] = 0\nprint(df[\"is_test\"].value_counts(), \"\\n\")\n\n# I tried many features combinations of features, I had always the same results : test AUC > 0.65 !\nfeat=[f\"sensor_{i+1}_d9_d1\" for i in range(10)] # Inter deciles (1st and 9th) differences for each sensor.\n\ntr, te, tr_y, te_y  = train_test_split(df[feat], df[\"is_test\"], test_size=0.33, random_state=__seed, shuffle=True)\n\nparam = {'num_leaves': 30,\n         'min_data_in_leaf': 30, \n         'max_depth': -1,\n         'learning_rate': 0.1,\n         \"min_child_samples\": 20,\n         \"feature_fraction\": 0.9,\n         \"bagging_freq\": 1,\n         \"bagging_fraction\": 0.9 ,\n         \"bagging_seed\": __seed,\n         \"boosting\": \"gbdt\",\n         'objective':'binary',\n#         \"metric\": 'binary_logloss',\n         \"metric\": 'auc',\n         \"verbosity\": -1}\n\nclf = lgbm.LGBMClassifier(**param, n_estimators = 1000, n_jobs = -1)\nclf.fit(tr, tr_y, eval_set = [(tr, tr_y), (te, te_y)], verbose=50, early_stopping_rounds = 100)\n\n\nfeature_importances = pd.DataFrame(clf.feature_importances_\n            , index = tr.columns, columns=['importance']).sort_values('importance', ascending=False)\n\nfeature_importances = feature_importances.reset_index()\nfeature_importances.columns = ['feature', 'importance']\n\nfig, ax = plt.subplots(figsize = (18, 8))\nsns.set()\nplt.subplot(1, 1, 1);\nsns.barplot(x=\"importance\", y=\"feature\", orient='h', data=feature_importances.head(50));\nplt.title('Feature Importance to detect Train or Test obs');","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}