{"cells":[{"metadata":{},"cell_type":"markdown","source":"# INGV - Volcanic Eruption Prediction\n\nIn this notebook I will describe the solution based on a simple workflow. This notebook will have the next sections:\n\n1. Feature engineering\n2. Experimentation\n3. Feature selection\n4. Hyperparameter tuning\n5. Features importances\n6. Submission"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":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 pickle\nimport time\nimport os\nimport datetime\nfrom tqdm import tqdm\nfrom pathlib import Path\n\nfrom tsfresh import extract_features\nfrom tsfresh.feature_extraction.settings import MinimalFCParameters\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.ensemble import StackingRegressor\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom sklearn.model_selection import RepeatedKFold\nfrom sklearn.tree import DecisionTreeRegressor\nfrom sklearn.ensemble import AdaBoostRegressor, GradientBoostingRegressor, RandomForestRegressor\nfrom sklearn.feature_selection import SelectFromModel\n\nimport lightgbm as lgb\nfrom catboost import Pool, CatBoostRegressor\n\nfrom skopt import BayesSearchCV\nfrom skopt.space import Real, Categorical, Integer\nimport eli5","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature engineering\n\nFor feature engineering I used the library [tsfresh](https://tsfresh.readthedocs.io/en/latest/). I customized the calculated features to be as efficient as it can. Note that more features could improve the results."},{"metadata":{"trusted":true},"cell_type":"code","source":"def feature_engineering(path, params, df=None, is_train=True):\n    features_df = None\n    \n    if is_train:\n        segments = df.iterrows()\n    else:\n        segments = enumerate(os.listdir(path))\n    \n    for idx, row in segments:\n        \n        if is_train:\n            segment, time_to_eruption = row\n            segment_timeseries_path = path/f'{segment}.csv'\n        else:\n            segment = row\n            segment_timeseries_path = path/f'{segment}'\n            \n        \n        segment_timeseries = pd.read_csv(segment_timeseries_path)\n        segment_timeseries=segment_timeseries.fillna(0).reset_index() # Drop nan columns and adds a \"time\" identifier loc[:,~segment_timeseries.isna().any().values]\n        segment_timeseries['id'] = idx # Every segment is from the same timeseries so lets asignate the same identifier for them\n        extracted_features = extract_features(segment_timeseries, column_id=\"id\", column_sort=\"index\", default_fc_parameters = params) # Extract the features ~ 3 minutes per segment\n        extracted_features['segment'] = segment\n        \n        if is_train:\n            extracted_features['time_to_eruption'] = time_to_eruption # Set the target\n\n\n        # Finally update the new dataset\n        if features_df is None:\n            features_df = extracted_features\n        else:\n            features_df = pd.concat([features_df,extracted_features], axis=0, ignore_index=True, sort=True)\n        \n    return features_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"volcanic_path = Path('/kaggle/input/predict-volcanic-eruptions-ingv-oe')\ntrain_path = volcanic_path/'train.csv'\nsubmission_path = volcanic_path/'sample_submission.csv'\ntrain_df = pd.read_csv(train_path)\nsubmission = pd.read_csv(submission_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params_tsfresh = MinimalFCParameters()\ndel params_tsfresh['length']\nparams_tsfresh['skewness'] = None\nparams_tsfresh['kurtosis'] = None\nparams_tsfresh['last_location_of_maximum'] = None\nparams_tsfresh['first_location_of_maximum'] = None\nparams_tsfresh['last_location_of_minimum'] = None\nparams_tsfresh['first_location_of_minimum'] = None\nparams_tsfresh['first_location_of_minimum'] = None\nparams_tsfresh['benford_correlation'] = None\nparams_tsfresh['number_peaks'] =  [{'n': 1}, {'n': 3}, {'n': 5}, {'n': 10}, {'n': 50}]\nparams_tsfresh['binned_entropy']  = [{'max_bins': 10}]\nparams_tsfresh['fft_aggregated']  = [{'aggtype': 'centroid'},\n  {'aggtype': 'variance'},\n  {'aggtype': 'skew'},\n  {'aggtype': 'kurtosis'}]\n\nparams_tsfresh['percentage_of_reoccurring_values_to_all_values'] = None\nparams_tsfresh['percentage_of_reoccurring_datapoints_to_all_datapoints'] = None\nparams_tsfresh[ 'autocorrelation'] =[{'lag': 0},\n  {'lag': 1},\n  {'lag': 2},\n  {'lag': 3},\n  {'lag': 4},\n  {'lag': 5},\n  {'lag': 6},\n  {'lag': 7},\n  {'lag': 8},\n  {'lag': 9}]\nparams_tsfresh['agg_autocorrelation'] = [{'f_agg': 'mean', 'maxlag': 40},\n  {'f_agg': 'median', 'maxlag': 40},\n  {'f_agg': 'var', 'maxlag': 40}]\nparams_tsfresh['friedrich_coefficients'] = [{'coeff': 0, 'm': 3, 'r': 30},\n  {'coeff': 1, 'm': 3, 'r': 30},\n  {'coeff': 2, 'm': 3, 'r': 30},\n  {'coeff': 3, 'm': 3, 'r': 30}]\nparams_tsfresh['count_above'] = [{'t': 0}]\nparams_tsfresh['count_below'] = [{'t': 0}]\nparams_tsfresh","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The calculation will take more than 2 hours so I will comment it. The results is in *features-volcanic* folder."},{"metadata":{"_kg_hide-output":false,"trusted":true},"cell_type":"code","source":"# features_train = feature_engineering(volcanic_path/'train', params_tsfresh, df)\n# features_test = feature_engineering(volcanic_path/'test', params_tsfresh, is_train=False) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Experimentation\n\nIn this section I will experiment with several models and select the best one, the tested models are:\n\n* KNN\n* AdaBoost\n* Gradient Boosting\n* Decission tree\n* Random Forest\n* LGBM\n* CatBoost"},{"metadata":{},"cell_type":"markdown","source":"First lets define our metrics."},{"metadata":{"trusted":true},"cell_type":"code","source":"EPSILON = 1e-10\n\ndef _error(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\" Simple error \"\"\"\n    return actual - predicted\n\n\ndef _percentage_error(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\"\n    Percentage error\n\n    Note: result is NOT multiplied by 100\n    \"\"\"\n    return _error(actual, predicted) / (actual + EPSILON)\n\ndef mse(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\" Mean Squared Error \"\"\"\n    return np.mean(np.square(_error(actual, predicted)))\n\ndef mae(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\" Mean Absolute Error \"\"\"\n    return np.mean(np.abs(_error(actual, predicted)))\n\ndef mape(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\"\n    Mean Absolute Percentage Error\n\n    Properties:\n        + Easy to interpret\n        + Scale independent\n        - Biased, not symmetric\n        - Undefined when actual[t] == 0\n\n    Note: result is NOT multiplied by 100\n    \"\"\"\n    return np.mean(np.abs(_percentage_error(actual, predicted)))\n\ndef rmse(actual: np.ndarray, predicted: np.ndarray):\n    \"\"\" Root Mean Squared Error \"\"\"\n    return np.sqrt(mse(actual, predicted))\n\ndef evaluate(actual: np.ndarray, predicted: np.ndarray):\n    results = {}\n    metrics = [('MAE', mae), ('MAPE', mape), ('RMSE', rmse)]\n    for name, metric in metrics:\n        try:\n            results[name] = metric(actual, predicted)\n        except Exception as err:\n            results[name] = np.nan\n            print('Unable to compute metric {0}: {1}'.format(name, err))\n    return results\n\ndef calculate_metrics(y_test, y_pred, target, model_name):\n\n    metrics = pd.DataFrame(data=evaluate(y_test, y_pred), index=[model_name])\n    \n    return metrics","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_scale_generator(df, target_name, n_splits=10, selector=None):\n    kf = KFold(n_splits=n_splits,random_state=123, shuffle=True)\n    y = df.loc[:, target_name]\n    X = df.drop(target_name, axis=1).values\n        \n    if selector is not None:\n        X = selector.transform(X)\n        \n    for train_index, test_index in kf.split(df.values):\n        scaler = StandardScaler()\n        \n        X_train, X_test = X[train_index], X[test_index]\n        y_train, y_test = y[train_index], y[test_index]\n        \n        X_train = scaler.fit_transform(X_train)\n        X_test = scaler.transform(X_test)\n        \n        yield X_train, y_train, X_test, y_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Train function for the model"},{"metadata":{"trusted":true},"cell_type":"code","source":"def train(name, model, target, X_train, y_train, X_test, y_test):\n    start_time = time.time()\n    print(\"Training \"+name+\"....\")\n\n    model.fit(X_train, y_train)\n\n    print(f\"Trained on {round((time.time()-start_time)/60, 3)} minutes....\")\n    \n    test_pred = model.predict(X_test)\n\n    test_metrics = calculate_metrics(y_test, test_pred, target, name)\n    \n    return test_metrics","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This function will train and evaluate all model with cross validation."},{"metadata":{"trusted":true},"cell_type":"code","source":"def execute_baseline(df, models, target,n_splits=10, selector=None):\n    valid_metrics_global = None\n\n    for idx, (X_train, y_train, X_test, y_test) in enumerate(split_scale_generator(df, target, n_splits,selector)):\n        print(f'******************** iteration {idx} ********************')\n        for name, model in models.items():\n            valid_metrics = train(name, model, target, X_train, y_train, X_test, y_test)\n            if valid_metrics_global is None:\n                valid_metrics_global = valid_metrics\n            else:\n                valid_metrics_global = valid_metrics_global.append(valid_metrics)\n        \n    return valid_metrics_global\n            ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#KNN\nknn_regr = KNeighborsRegressor(n_neighbors=3, weights='distance', n_jobs=-1)\n\n# AdaBoost with Decission tree\n#distributions_ada = {'estimator__base_estimator__max_depth':list(range(2,16))}\nada_regrt_model = AdaBoostRegressor(DecisionTreeRegressor(),\n                          n_estimators=300,  random_state=123)\n\n\n# Gradient boosting\ngrboost_model = GradientBoostingRegressor(n_estimators=300,loss='ls', learning_rate=0.1, random_state=123)\n\n#distributions_dt = dict(max_depth=list(range(2,16)))\nregrt_model = DecisionTreeRegressor(random_state=123)\n\n#Random forest\nrand_forest_model = RandomForestRegressor(n_estimators=300, n_jobs=-1, random_state=123)\n\n\n# Light GBM\nlgb_params = {\n    \n    'min_data_in_leaf': 2, \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\": 'mae',\n    \"lambda_l1\": 0.1,\n    \"verbosity\": -1,\n    \"nthread\": -1,\n    \"random_state\": 123\n}\nlgbm_regr = lgb.LGBMRegressor(**lgb_params, n_estimators = 300, n_jobs = -1)\n\n# CatBoost                    \ncatboost_regr = CatBoostRegressor(learning_rate=1, random_seed=123, n_estimators=600, logging_level='Silent')\n\n\nmodels = {\n    \"knn_regr\": knn_regr,\n    \"adaboost_tree_regressor\": ada_regrt_model,\n    \"gradient_boost\": grboost_model,\n    \"decision_tree_regressor\": regrt_model,\n    \"random_forest\": rand_forest_model,\n    \"lgbm_regressor\": lgbm_regr,\n    \"catboost_regressor\": catboost_regr,\n    \n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfeatures_path = Path('/kaggle/input/features-volcanic')\n\ntrain_features = pd.read_csv(features_path/'features_train.csv')\nmetrics = execute_baseline(train_features, models, 'time_to_eruption', n_splits=10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets see the metrics"},{"metadata":{"trusted":true},"cell_type":"code","source":"metrics.groupby(metrics.index).mean().sort_values('MAE')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Adaboost solution is the best solution, if I submit this model the MAE is 4749935.\nLets improve it with feature selection and hyperparameter tuning."},{"metadata":{},"cell_type":"markdown","source":"## Feature selection"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(features_path/'features_train.csv')\ntest = pd.read_csv(features_path/'features_test.csv')\n\nX_test = test.drop(['segment', 'time_to_eruption'], axis=1).values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nmodel = AdaBoostRegressor(DecisionTreeRegressor(), n_estimators=300)\n\nscaler = StandardScaler()\ny = train.loc[:, 'time_to_eruption']\nX = scaler.fit_transform(train.drop('time_to_eruption', axis=1).values)\nX_test = scaler.transform(X_test)\n\nselector = SelectFromModel(estimator=model).fit(X, y)\nX_test_selected = selector.transform(X_test)\nX_train_selected = selector.transform(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"selected_features = train.drop('time_to_eruption', axis=1).columns[selector.get_support()].tolist()\nselected_features","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With feature selection the MAE improved to 4596505."},{"metadata":{},"cell_type":"markdown","source":"## Hyperparameter tuning"},{"metadata":{},"cell_type":"markdown","source":"For hyperparameter tuning I used Bayesian Optimization from [skopt](https://scikit-optimize.github.io/stable/)."},{"metadata":{"trusted":true},"cell_type":"code","source":"search_space = {\"loss\": Categorical(['linear', 'square', 'exponential']), # values for boostrap can be either True or False\n        \"base_estimator__max_depth\": Integer(3, 200), # values of max_depth are integers from 6 to 20\n        \"base_estimator__splitter\": Categorical(['best', 'random']), \n        \"n_estimators\": Integer(200, 5000),\n        \"learning_rate\": Real(1e-12, 1)\n    }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = AdaBoostRegressor(DecisionTreeRegressor())\n\nclf = BayesSearchCV(model, search_space, cv=10, n_iter=15, n_jobs=-1, random_state=123, scoring='neg_mean_absolute_error')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This will take a lot of time so lets comment this cell."},{"metadata":{"trusted":true},"cell_type":"code","source":"#%%time\n# search = clf.fit(X_train_selected, y)\n\n# search.best_params_, search.best_score_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The best params I found are:\n* base_estimator__max_depth: 43\n* base_estimator__splitter: random\n* learning_rate: 0.09210493994507518\n* loss: linear\n* n_estimators: 943"},{"metadata":{},"cell_type":"markdown","source":"## Features importances"},{"metadata":{"trusted":true},"cell_type":"code","source":"model = AdaBoostRegressor(DecisionTreeRegressor(max_depth=68, splitter='random'), learning_rate=.162163403819552, n_estimators=971, loss='exponential')\n\nperm = eli5.sklearn.PermutationImportance(model, cv=10).fit(X_train_selected, y)\neli5.show_weights(perm, feature_names=selected_features)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission"},{"metadata":{},"cell_type":"markdown","source":"Finally lets see the MAE with hyperparameter tuning and selection."},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"#model = AdaBoostRegressor(DecisionTreeRegressor(max_depth=search.best_params_['base_estimator__max_depth'], splitter=search.best_params_['base_estimator__splitter']), learning_rate=search.best_params_['learning_rate'], n_estimators=search.best_params_['n_estimators'], loss=search.best_params_['loss'])\nmodel = AdaBoostRegressor(DecisionTreeRegressor(max_depth=68, splitter='random'), learning_rate=.162163403819552, n_estimators=971, loss='exponential')\n\nmodel.fit(X_train_selected, y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = model.predict(X_test_selected)\nsubmission = pd.read_csv(submission_path)\nsubmission['time_to_eruption'] = predictions\nsubmission.to_csv('submission.csv', index=None)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The final MAE in private leaderboard was: 3928087."}],"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}