{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook summarizes the main steps I followed for this dataset. I mainly compared tree ensembling methods and neural networks.\n\nMachine learning algorithms cannot be straightforwardly applied on the given data. We first have to perform some feature extraction. Three functions are introduced below :\n- ListFiles: this function takes as argument a directory and lists all the csv files present in this directory. It then returns a list containing all csv corresponding to volcanoes, a list containing the id's of the volcanoes and a DataFrame corresponding to the train.csv file.\n- FeatureExtraction: this function takes as argument a DataFrame containing sensor data, and extracts features from it. We aim at summarizing efficiently the information contained in the sensor data. The extracted features are namely the mean, standard deviation, max, quantiles, skewness and kurtosis of each sensor. The positions of the maximum and minimum values are also extracted.\n- ImportData: this function takes as arguments a directory containing the csv files corresponding to each volcano, and the file containing the response variable. It then performs the extraction of features on each csv files corresponding to a volcano. It finally returns a DataFrame containing the extracted features for each volcano, and another DataFrame containing the id's and the remaining time to eruption (the response variable)."},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os\n\ndef ListFiles(path):\n    listfiles = []\n    id = []\n    for root, dir, files in os.walk(path):\n        for x in files:\n            listfiles.append(os.path.join(root, x))\n            id.append(x)\n    id = [x.split(\".csv\")[0] for x in id]\n    ord = np.argsort(np.array(id))\n    files = np.array(listfiles)[ord]\n    id = np.array(id)[ord]\n    return files, id\n\nfrom scipy.stats import skew, kurtosis\ndef FeatureExtraction(X):\n    sumNA = np.apply_along_axis(np.sum,0,(X.isnull()).values)\n    X = (X.fillna(0)).values\n    xmax = np.apply_along_axis(np.max,0,X)\n    xmean = np.apply_along_axis(np.mean,0,X)\n    xq = np.apply_along_axis(lambda x: np.quantile(x,np.array([0.01,0.05,0.1,0.25,0.5,0.75,0.9,0.95,0.99])),0,X).T.reshape(90)\n    xstd = np.apply_along_axis(np.std,0,X)\n    xskew = np.apply_along_axis(skew,0,X)\n    xkurt = np.apply_along_axis(kurtosis, 0, X)\n    xargmax = np.apply_along_axis(np.argmax,0,X)\n    xmin = np.apply_along_axis(np.min, 0, X)\n    xargmin = np.apply_along_axis(np.argmin, 0, X)\n    res = np.concatenate([xmax,xmean,xstd,xskew,xkurt,xmin,xargmax,xargmin,sumNA],axis=0)\n    res = np.concatenate([res,xq],axis=0)\n    return res\n\ndef ImportData(path,Yfile):\n    Y = pd.read_csv(Yfile,sep=\",\")\n    Y[\"segment_id\"] = [str(x) for x in Y[\"segment_id\"].values]\n    Y = Y.iloc[np.argsort(Y[\"segment_id\"]),:]\n    files, id = ListFiles(path)\n    X = np.zeros((len(files), 9*10+9*10))\n    for i in range(len(files)):\n        if i % 500 == 0:\n            print(str(round(i * 100 / len(files))) + \"% \")\n        XX = pd.read_csv(files[i], sep=\",\")\n        X[i, :] = FeatureExtraction(XX)\n    cols1 = [\"max\", \"mean\", \"std\", \"skew\", \"kurt\", \"min\", \"argmax\", \"argmin\", \"sumNA\"]\n    cols2 = [\"Q\" + x for x in [\"1\", \"5\", \"10\", \"25\", \"50\", \"75\", \"90\", \"95\", \"99\"]]\n    cols = cols1 + cols2\n    cols = [[x]*10 for x in cols]\n    cols = [x for l in cols for x in l]\n    S = [\"S\" + str(i) + \"_\" for i in range(1,11)]*int(len(cols)/10)\n    cols = [S[i] + cols[i] for i in range(len(cols))]\n    X = pd.DataFrame(X, columns=cols, index=id)\n    ind_const = np.where(X.std() == 0)[0]\n    X = X.drop(columns=X.columns[ind_const])\n    cols = [i for i in range(X.shape[1]) if \"max\" in X.columns[i] and not \"argmax\" in X.columns[i]]\n    X[\"Maxmax\"] = X[X.columns[cols]].apply(\"max\",axis=1)\n    cols = [i for i in range(X.shape[1]) if \"std\" in X.columns[i]]\n    X[\"Maxstd\"] = X[X.columns[cols]].apply(\"max\", axis=1)\n    cols = [i for i in range(X.shape[1]) if \"skew\" in X.columns[i]]\n    X[\"Maxskew\"] = X[X.columns[cols]].apply(\"max\", axis=1)\n    cols = [i for i in range(X.shape[1]) if \"kurt\" in X.columns[i]]\n    X[\"Maxkurt\"] = X[X.columns[cols]].apply(\"max\", axis=1)\n    cols = [i for i in range(X.shape[1]) if \"mean\" in X.columns[i]]\n    X[\"Maxmean\"] = X[X.columns[cols]].apply(\"max\", axis=1)\n    cols = [i for i in range(X.shape[1]) if \"min\" in X.columns[i] and not \"argmin\" in X.columns[i]]\n    X[\"Maxmin\"] = X[X.columns[cols]].apply(\"max\", axis=1)\n    X[\"Minmin\"] = X[X.columns[cols]].apply(\"min\", axis=1)\n    return X, Y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"[](http://)These functions can be used to build the dataset using the code below. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#path = \"../input/predict-volcanic-eruptions-ingv-oe/train\"\n#file = \"../input/predict-volcanic-eruptions-ingv-oe/train.csv\"\n#X, Y = ImportData(path,file)\n#X.to_csv(\"X.csv\",sep=\";\",columns=X.columns)\n#Y.to_csv(\"Y.csv\",sep=\";\",columns=Y.columns)\nX = pd.read_csv(\"../input/preprocessed-volcano-data/X.csv\",sep=\";\",index_col=0)\nY = pd.read_csv(\"../input/preprocessed-volcano-data/Y.csv\",sep=\";\",index_col=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can now start evaluating algorithms. We first split the dataset in a train and a test sets,  the test set will contain 500 instances of the whole dataset. We then perform linear regression only using the max features, in order to obtain a baseline performance."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nXtrain, Xtest, Ytrain, Ytest = train_test_split(X,Y[\"time_to_eruption\"],test_size=500,random_state=1234)\n\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import mean_absolute_error\n\nreg = LinearRegression()\ncols = [col for col in X.columns if \"max\" in col and not \"argmax\" in col]\nreg.fit(Xtrain[cols].values, Ytrain)\npred = reg.predict(Xtest[cols].values)\nmae = np.mean(np.abs(Ytest-pred))\nprint(mae)\nR2 = np.corrcoef(Ytest,pred)[0,1]**2\nprint(R2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The mean absolute error approximately equals 10^7. According to the R-squared coefficient, the correlation between the predictors and the response is weak. However, there might be a strong non-linear link. It might also be useful to include the remaining features. Let us try using all these features using a RandomForest algorithm."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import GridSearchCV\nfrom sklearn.ensemble import RandomForestRegressor\ngridsearch = GridSearchCV(estimator=RandomForestRegressor(),\n                          param_grid={\"max_features\" : [50, 100, 150],\n                                      \"min_samples_leaf\" : [1, 5],\n                                      \"n_estimators\" : [1000]},\n                          cv=5,n_jobs=5,scoring=\"neg_mean_absolute_error\",verbose=10)\ngridsearch.fit(Xtrain, Ytrain)\nres_RF = pd.DataFrame({\"max_features\" : [list(x.values())[0] for x in gridsearch.cv_results_[\"params\"]],\n                       \"min_samples_leaf\" : [list(x.values())[1] for x in gridsearch.cv_results_[\"params\"]],\n                       \"mean_mae\" : np.round(gridsearch.cv_results_[\"mean_test_score\"],4),\n                       \"sd_mae\" : np.round(gridsearch.cv_results_[\"std_test_score\"],4)})\nprint(res_RF)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The cross-validation scores obtained with the RandomForest algorithm are clearly lower than the score obtained on the test set with linear regression. This demonstrates that the chosen features are useful for predicting the remaining time to eruption. I then tried the ExtraTrees (extremely randomized trees) algorithm."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import ExtraTreesRegressor\ngridsearch = GridSearchCV(estimator=ExtraTreesRegressor(),\n                          param_grid={\"max_features\" : [50, 100, 150],\n                                      \"min_samples_leaf\" : [1, 5],\n                                      \"n_estimators\" : [1000]},\n                          cv=5,n_jobs=5,scoring=\"neg_mean_absolute_error\",verbose=10)\ngridsearch.fit(Xtrain, Ytrain)\nres_ET = pd.DataFrame({\"max_features\" : [list(x.values())[0] for x in gridsearch.cv_results_[\"params\"]],\n                       \"min_samples_leaf\" : [list(x.values())[1] for x in gridsearch.cv_results_[\"params\"]],\n                       \"mean_mae\" : np.round(gridsearch.cv_results_[\"mean_test_score\"],4),\n                       \"sd_mae\" : np.round(gridsearch.cv_results_[\"std_test_score\"],4)})\nprint(res_ET)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The results are better than those obtained with the RandomForest. After fine tuning both algorithms, I obtained the following optimal values for the hyperparameters, ExtraTrees still gave better results than RandomForest. The obtained values of the hyperparameters were max_features = 164 and min_node_size = 1. Let's evaluate it on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import ExtraTreesRegressor\nET = ExtraTreesRegressor(max_features=164,min_samples_leaf=1)\nET.fit(Xtrain, Ytrain)\npred = ET.predict(Xtest)\nmae = np.mean(np.abs(Ytest-pred))\nprint(mae)\nR2 = np.corrcoef(Ytest,pred)[0,1]**2\nprint(R2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These results start looking interesting. At this point, several boosting algorithms (CatBoost, LightGBM, XGBoost) were also compared. All of them had a lower performance than ExtraTrees. \n\nTo further improve the performance, I tried to perform some feature selection. It is very likely that some features are noisy or redundant, which might deter the performance of the regression algorithms. The feature selection procedure is the following :\n- Entries : number of iterations n_iter (30 is enough here), matrix of predictors X and response vector Y\n- For i = 1 to 30 :\n    - Create an ExtraTrees algorithm with max_features = number of features of X and min_node_size = 1\n    - Evaluate the algorithm using 5-fold cross-validation, get the mean, min and max CV scores\n    - Fit the algorithm on the full set to get the feature importances\n    - Remove from X the 5 least important features\n- The result is a dictionary whose elements are :\n    - remaining_cols : a list giving the names of the features remaining at each step\n    - cvmeanMAE, cvminMAE, cvmaxMAE : the mean, min and max CV MAE obtained at each step.\nThis procedure is implemented in the function below. After running it, one can plot the CV scores against the number of feature selection steps. We will then select the remaining features at the step giving the best CV scores.\n\nThe value of max_features at each step was chosen to be the total number of available features since in the previous results, the optimal performance was obtained with max_features very close to the total number of features. At each step, a batch of 5 features is removed because the number of features is quite important, so removing them one by one would take too much time. Moreover, if there truly are many useless features, removing batches of 5 features will not do any harm at the beginning."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold, cross_validate\nfrom sklearn.ensemble import ExtraTreesRegressor\ndef FeatureSelection(X,Y,nfolds=5,random_state=123):\n    XX = X.values\n    n_iter = 30\n    cvmeanMAE = np.zeros(n_iter)\n    cvminMAE = np.zeros(n_iter)\n    cvmaxMAE = np.zeros(n_iter)\n    cols = X.columns\n    remaining_cols = []\n    cvsplit = KFold(n_splits=nfolds, shuffle=True, random_state=random_state)\n    for i in range(n_iter):\n        print(\"Variable selection step \" + str(i + 1) + \" - \" + str(XX.shape[1]) + \" features remaining\")\n        rf = ExtraTreesRegressor(max_features=XX.shape[1], min_samples_leaf=1, n_estimators=1000)\n        cv = cross_validate(estimator=rf, X=XX, y=Y, cv=cvsplit, n_jobs=5,\n                            scoring=\"neg_mean_absolute_error\")\n        cvmeanMAE[i] = np.mean(-cv[\"test_score\"])\n        cvminMAE[i] = np.min(-cv[\"test_score\"])\n        cvmaxMAE[i] = np.max(-cv[\"test_score\"])\n        rf.fit(XX, Y)\n        varimp = rf.feature_importances_\n        ind = np.argsort(varimp)[5:]\n        XX = XX[:, ind]\n        cols = cols[ind]\n        remaining_cols.append(cols)\n    res = {\"remaining_cols\" : remaining_cols,\n           \"cvmeanMAE\" : cvmeanMAE,\n           \"cvminMAE\" : cvminMAE,\n           \"cvmaxMAE\" : cvmaxMAE}\n    return res\n\ntmp = FeatureSelection(Xtrain, Ytrain)\nfrom matplotlib import pyplot\npyplot.bar(np.arange(len(tmp[\"cvmeanMAE\"])),height=tmp[\"cvmeanMAE\"],\n           yerr=np.concatenate([(tmp[\"cvmeanMAE\"]-tmp[\"cvminMAE\"]).reshape(-1,1),\n                                (tmp[\"cvmaxMAE\"]-tmp[\"cvmeanMAE\"]).reshape(-1,1)],axis=1).T)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It can be seen that the MAE slowly decreases as the number of feature selection steps increases, until a number of steps around 25, after which it starts to increase again. Below are given the remaining columns when I ran it for the first time."},{"metadata":{"trusted":true},"cell_type":"code","source":"sel_cols = ['S1_max', 'S10_Q75', 'S3_Q50', 'S1_sumNA', 'S3_Q1', 'S3_max', 'S9_std',\n       'S6_Q5', 'S7_Q1', 'S10_Q95', 'S5_Q90', 'S2_Q5', 'S5_Q5', 'S6_Q10',\n       'S7_Q50', 'S1_Q95', 'S3_Q75', 'S5_Q10', 'S8_Q99', 'S4_Q95', 'S3_sumNA',\n       'S7_Q90', 'S8_Q50', 'S6_Q90', 'S10_Q90', 'S2_sumNA', 'S4_Q99', 'S2_Q75',\n       'S9_Q90', 'S9_Q25', 'S1_Q99', 'S6_Q1', 'S1_Q5', 'S7_Q99', 'S3_Q5',\n       'S3_Q95', 'S5_Q99', 'S4_Q1', 'S9_Q50', 'S2_std', 'S1_Q75', 'S2_Q50',\n       'S7_Q5', 'S10_Q25', 'S5_sumNA', 'S8_Q5', 'S10_Q1']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can now select these features and evaluate the former algorithms again. As previously, the best results were obtained using ExtraTrees (here with max_features = 47 and min_node_size = 1)."},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain = Xtrain[sel_cols]\nXtest = Xtest[sel_cols]\nET = ExtraTreesRegressor(max_features=47,min_samples_leaf=1,n_estimators=1000)\nET.fit(Xtrain, Ytrain)\npred = ET.predict(Xtest)\nmae = np.mean(np.abs(Ytest-pred))\nprint(mae)\nR2 = np.corrcoef(pred,Ytest)[0,1]**2\nprint(R2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The feature selection procedure led to a futher improvement of the results. Using only the selected features, the previous algorithms (RandomForest and boosting algorithms) were compared to ExtraTrees again, but ExtraTrees still gave the best results (this might be biased at this point since ExtraTrees was used to select the features).\n\nI then switched to neural networks. I first tried out on the full set of variables, but could not get a satisfying performance. Trying again on the selected features gave much better results. A class implementing all the required steps for the training of the network and for predicting the response is given below. Those steps namely include proper scaling of the predictors and response and two callbacks:\n- ModelCheckpoint, which enables us to save a fitted model at a given epoch when its score on the validation step beats the previously obtained optimal validation score\n- ReduceLROnPlateau, which automatically multiplies by 0.5 the current value of the learning rate if no improvement on the validation score is observed during 20 successive epochs.\n\nThe arguments of the Network class are the following :\n- file : name of the file to which will be saved the network by the ModelCheckpoint callback, each time an improvement is observed on the validation set\n- nval : size of the validation set\n- epochs : number of epochs\n- batch_size : batch size used for training the network."},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow import keras\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import train_test_split\nclass Network:\n    def __init__(self,file,nval=500,nepochs=300,batch_size=64):\n        self.nval = nval\n        self.nepochs = nepochs\n        self.batch_size = batch_size\n        self.file = file\n    def __build_model(self,input_shape):\n        net = keras.models.Sequential([\n            keras.layers.Input(shape=input_shape),\n            keras.layers.Dense(256, activation=\"relu\"),\n            keras.layers.Dense(256, activation=\"relu\"),\n            keras.layers.Dense(256, activation=\"relu\"),\n            keras.layers.Dense(256, activation=\"relu\"),\n            keras.layers.Dense(512, activation=\"relu\"),\n            keras.layers.Dropout(0.5),\n            keras.layers.Dense(1, activation=\"linear\")\n        ])\n        return net\n    def fit(self,X,Y):\n        Xtrain, Xval, Ytrain, Yval = train_test_split(X, Y, test_size=self.nval)\n        scaler = StandardScaler()\n        scaler.fit(Xtrain)\n        Xtrain_sc = scaler.transform(Xtrain)\n        Xval_sc = scaler.transform(Xval)\n        my = np.mean(Ytrain)\n        sy = np.std(Ytrain)\n        Ytrain_sc = (Ytrain - my) / sy\n        Yval_sc = (Yval - my) / sy\n        net = self.__build_model(X.shape[1])\n        checkpoint = keras.callbacks.ModelCheckpoint(self.file,monitor=\"val_loss\",\n                                                     save_best_only=True,mode=\"min\",\n                                                     verbose=1)\n        lr_callback = keras.callbacks.ReduceLROnPlateau(factor=0.5, patience=20)\n        net.compile(loss=\"mean_absolute_error\", optimizer=keras.optimizers.Nadam(learning_rate=0.001))\n        fitted = net.fit(x=Xtrain_sc, y=Ytrain_sc, batch_size=self.batch_size, \n                         epochs=self.nepochs, verbose=1,\n                         callbacks=[lr_callback, checkpoint], \n                         validation_data=(Xval_sc, Yval_sc))\n        net = keras.models.load_model(self.file)\n        fitted = {\"net\" : net, \"scaler\" : scaler, \"my\" : my, \"sy\" : sy}\n        self.fitted = fitted\n    def predict(self,Xtest):\n        fitted = self.fitted\n        net = fitted[\"net\"]\n        scaler = fitted[\"scaler\"]\n        my = fitted[\"my\"]\n        sy = fitted[\"sy\"]\n        pred = net.predict(scaler.transform(Xtest))[:,0]*sy+my\n        return pred","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The network can now be trained on the training set and evaluated on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"net = Network(file=\"fitted_network.h5\",nval=500,nepochs=300,batch_size=64)\nnet.fit(Xtrain,Ytrain)\npred = net.predict(Xtest)\nmae = np.mean(np.abs(Ytest-pred))\nprint(mae)\nR2 = np.corrcoef(Ytest,pred)[0,1]**2\nprint(R2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So far, the best performance was obtained with this neural network. We can retrain it on the full dataset and predict the values for the test dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"X = pd.read_csv(\"../input/preprocessed-volcano-data/X.csv\",sep=\";\",index_col=0)\nY = pd.read_csv(\"../input/preprocessed-volcano-data/Y.csv\",sep=\";\",index_col=0)\nXtest = pd.read_csv(\"../input/preprocessed-volcano-data/Xtest.csv\",sep=\";\",index_col=0)\nsegment_id = Xtest.index\nX = X[sel_cols]\nXtest = Xtest[sel_cols]\nnet = Network(file=\"fitted_network.h5\",nval=500,nepochs=300,batch_size=64)\nnet.fit(X,Y[\"time_to_eruption\"])\npred = net.predict(Xtest)\nres = pd.DataFrame({\"segment_id\" : segment_id, \"time_to_eruption\" : pred})\nres.to_csv(\"prediction_test.csv\",sep=\",\",columns=res.columns,index=False)","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}