{"cells":[{"metadata":{},"cell_type":"markdown","source":"

INGV - Volcanic Eruption Prediction"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nfrom IPython.display import clear_output\nimport gc\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import r2_score\nfrom sklearn.feature_selection import RFE\nfrom sklearn.pipeline import Pipeline\nfrom scipy.stats import moment","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files = os.listdir(\"../input/predict-volcanic-eruptions-ingv-oe/train\")\ntest_files = os.listdir(\"../input/predict-volcanic-eruptions-ingv-oe/test\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(train_files)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(test_files)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Checking if the shape of data is uniform across all the train files"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = []\nrows = []\nfor i,fname in enumerate(train_files):\n train = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/train\",fname))\n cols.append(train.shape[0])\n rows.append(train.shape[1])\n print(f'{i+1} / {len(train_files)}')\n clear_output(wait=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"Rows of all train files: {pd.Series(rows).unique()}\\nColumns of all train files: {pd.Series(cols).unique()}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Checking if the shape of data is uniform across all the test files"},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = []\nrows = []\nfor i,fname in enumerate(test_files):\n test = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/test\",fname))\n cols.append(test.shape[0])\n rows.append(test.shape[1])\n print(f'{i+1} / {len(test_files)}')\n clear_output(wait=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"Rows of all test files: {pd.Series(rows).unique()}\\nColumns of all test files: {pd.Series(cols).unique()}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

All the train and test files have same shape. Hence, we'll inspect the first few rows of first 10 train files to get an idea of data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/train\",train_files[2]))\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Whoa!! NaNs"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.count()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Not only NaNs there are also few 0s"},{"metadata":{"trusted":true},"cell_type":"code","source":"(train==0).sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

There are few sensors with no data at all in some files. Let us also examine if there are sensors with intermittent missing values (NaN) across files"},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_tracker = pd.DataFrame()\nfor i,fname in enumerate(train_files):\n train = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/train\",fname))\n missing_tracker = missing_tracker.append(pd.DataFrame(train.count()).T)\n print(f'{i+1} / {len(train_files)}')\n clear_output(wait=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_tracker","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

There are sensors with intermittent missing values"},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_tracker.nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for col in missing_tracker.columns:\n print(f\"{col}\\n{sorted(missing_tracker[col].unique())}\\n\\n\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    \n
  1. sensor_1: Has moderate no. of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    2. \n
  2. sensor_2: Has lots of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    3. \n
  3. sensor_3: Has moderate no. of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    4. \n
  4. Sensor_4: Looks quite stable. It's having very less number of intermittent NaNs and there are no cases where it hasn't recorded any readings.
    5. \n
  5. sensor_5: Has lots of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    6. \n
  6. sensor_6: Has moderate no. of intermittent missing readings (few NaNs) and no cases with no readings (all NaNs)
    7. \n
  7. sensor_7: Has less no. of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    8. \n
  8. sensor_8: Has very less no. of intermittent missing readings (few NaNs) and cases with no readings (all NaNs)
    9. \n
  9. Sensor_9: Seems to be highly unstable with varying intermittent missing readings / NaNs
    10. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,20))\nsns.heatmap(missing_tracker);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    The heatmap shows sensors 2, 3, 5, 8, and 9 have no readings in many files. Sensors 4 and 6 have atleast 1 reading in every file."},{"metadata":{},"cell_type":"markdown","source":"

    Let's see the distribution of sensors with complete readings in the files"},{"metadata":{"trusted":true},"cell_type":"code","source":"mi = ((missing_tracker==60001).sum()/4431)*100\nprint(f'% of files with complete data per sensor\\n\\n{mi}')\nmi.plot(kind=\"bar\")\nplt.axhline(y=100,color=\"red\")\nplt.title(\"Sensors with complete readings\")\nplt.ylabel(\"% of files with complete data\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    The above bar chart shows sensors 1,4,6, and 7 have almost all the readings in all the files

    \n\n

    Let's see the distribution of sensors with no readings at all across the files"},{"metadata":{"trusted":true},"cell_type":"code","source":"mi = ((missing_tracker==0).sum()/4431)*100\nprint(f'% of files with no data at all per sensor\\n\\n{mi}')\nmi.plot(kind=\"bar\")\nplt.title(\"Sensors with no readings at all\")\nplt.ylabel(\"% of files with no data at all\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    As seen earlier in the heatmap, sensors 2, 3, 5, 8, and 9 have no readings in many files."},{"metadata":{},"cell_type":"markdown","source":"

    Let's look at the data in a single training file"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_1 = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/train\",train_files[0]))\ntrain_1.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_1.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Let's plot the readings of each sensor"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,ax=plt.subplots(5,2,figsize=(12,17))\nr=0\nc=0\nfor i,col in enumerate(train_1.columns):\n ax[r,c].plot(train_1[col])\n ax[r,c].set_title(col)\n c+=1\n if (i+1)%2==0:\n r+=1\n c=0\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    All the sensors don't have the readings in the same range"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_1.plot(figsize=(15,8))\nplt.legend(loc=\"upper left\", ncol=5);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    In the above sample, sensors 1 and 2 have many peaks and troughs in the readings, especially sensor 2"},{"metadata":{},"cell_type":"markdown","source":"

    Are the sensor readings correlated. How's the distribution of sensor readings?

    "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10,7))\nsns.heatmap(train_1.corr(),annot=True,fmt=\".2f\",cbar=False);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    There is no correlation between the sensor readings"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(train_1,diag_kind=\"kde\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Reducing the dimensionality of train and test sets by computing important stats (using describe) on monthly data

    \n\n

    There are 4,000+ train files/locations with 60,000 readings of 10 sensors. Each loacation/file has a single target variable i.e. 60,000 rows of data has only 1 target/label. Hence, we need to reshape the data into single row per location. Doing this for huge data can be tedious. The readings are recorded for every 10 minutes. Hence, I've applied the describe function on monthly data and reshaped it into a single row."},{"metadata":{"trusted":true},"cell_type":"code","source":"def find_stats(x):\n '''\n Function to apply describe on monthly data\n '''\n return x.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Selecting 70% of train files due to memory constraints"},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\nrandom.seed(11)\nrand_sample = random.sample(range(len(train_files)),int(np.floor(len(train_files)*0.7)))\nprint(len(rand_sample))\ntrain_files = list(pd.Series(train_files)[rand_sample])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.DataFrame()\nfor n,i in enumerate(train_files):\n df = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/train\",i))\n df.index=pd.date_range(\"2000-01-01\",periods=len(df),freq=\"10T\") #adding date range with 10 minute interval\n df = df.resample(rule=\"M\").apply(find_stats).values.reshape(1,-1) #resampling the time series data and calculating stats on monthly data\n df = pd.DataFrame(df)\n df[\"segment_id\"] = int(i.replace(\".csv\",\"\"))\n train = train.append(df)\n print(f'{n}')\n gc.collect()\n clear_output(wait=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_labels = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = train.merge(train_labels,on=\"segment_id\",how=\"left\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train=train.drop(columns=\"segment_id\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.fillna(0,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train.iloc[:,:-1].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = train.iloc[:,-1].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#from sklearn.model_selection import GridSearchCV\n#from sklearn.decomposition import PCA\nfrom lightgbm import LGBMRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.model_selection import KFold,cross_val_score","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Though PCA helped me reduce the dimensionality with 99.99 explained variance ratio. But the performance of the model was better without PCA being applied since it ignores the target variable"},{"metadata":{"trusted":true},"cell_type":"code","source":"#pca = PCA(n_components=10)\n#X_pca = pca.fit_transform(X)\n#np.cumsum(pca.explained_variance_ratio_)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Since the dataset is huge, RAM utilization reached to its maximum, for better memory management I've manually removed unwanted variables and performed garbage collection frequently"},{"metadata":{"trusted":true},"cell_type":"code","source":"del df\ndel train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_pca = X.copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Using Random Forest for feature importance"},{"metadata":{"trusted":true},"cell_type":"code","source":"rf = RandomForestRegressor(random_state=11)\nrf.fit(X_pca,y)\nprint(f\"No of Features with Importance > 0: {sum(rf.feature_importances_>0)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Selecting threshold for feature importance. Importance above 75 percentile"},{"metadata":{"trusted":true},"cell_type":"code","source":"feat_imp = pd.Series(rf.feature_importances_).reset_index(drop=True)\nthreshold = feat_imp.quantile(0.75)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"final_feat = list(feat_imp[feat_imp>threshold].index)\nlen(final_feat)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_pca = X_pca.iloc[:,final_feat].copy()\nX_pca.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del rf\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Manually tuning the hyperparameters, since GridSearchCV is taking too long"},{"metadata":{"trusted":true},"cell_type":"code","source":"MAX_DEPTH = 10\nN_ESTIMATORS = 1000\nMIN_SAMPLES_LEAF = 300\nL1 = 50000\n\n\ngscv = LGBMRegressor(n_estimators=N_ESTIMATORS,\n max_depth=MAX_DEPTH,\n num_leaves=2**MAX_DEPTH,\n min_data_in_leaf=MIN_SAMPLES_LEAF,\n lambda_l1 = L1,\n random_state=11,\n n_jobs=-1)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cv = KFold(n_splits=3,random_state=11,shuffle=True)\ncv_score = cross_val_score(gscv,X_pca,y,cv=cv,scoring=\"neg_mean_absolute_error\",n_jobs=-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    CV scores (MAE)"},{"metadata":{"trusted":true},"cell_type":"code","source":"-1 * (cv_score.astype(\"int\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.mean(-1*(cv_score.astype(\"int\")))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gscv.fit(X_pca,y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del X\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    MAE on train set"},{"metadata":{"trusted":true},"cell_type":"code","source":"mean_absolute_error(gscv.predict(X_pca),y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    R2 on Training set"},{"metadata":{"trusted":true},"cell_type":"code","source":"r2_score(gscv.predict(X_pca),y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

    Model is overfit to the training set. Need further hyperparameter tuning."},{"metadata":{"trusted":true},"cell_type":"code","source":"del cv, cv_score\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del X_pca\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.DataFrame()\nfor n,i in enumerate(test_files):\n df = pd.read_csv(os.path.join(\"../input/predict-volcanic-eruptions-ingv-oe/test\",i))\n df.index=pd.date_range(\"2000-01-01\",periods=len(df),freq=\"10T\")\n df = df.resample(rule=\"M\").apply(find_stats).values.reshape(1,-1)\n df = pd.DataFrame(df)\n df = df.iloc[:,final_feat].copy()\n df[\"segment_id\"] = int(i.replace(\".csv\",\"\"))\n test = test.append(df)\n print(f'{n}')\n gc.collect()\n clear_output(wait=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_segment_ids = test[\"segment_id\"]\ntest = test.drop(columns=\"segment_id\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.fillna(0,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#X_test_pca = pca.transform(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test_pca = test.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#X_test_pca = X_test_pca.iloc[:,final_feat].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = gscv.predict(X_test_pca)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.DataFrame({\"segment_id\":test_segment_ids,\"time_to_eruption\":pred})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv(\"ingv_submission.csv\",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}