{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport gc\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import mean_squared_error\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nfrom sklearn.linear_model import SGDRegressor\nfrom sklearn.linear_model import Lasso\nfrom sklearn.linear_model import HuberRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.preprocessing import StandardScaler\nimport xgboost as xgb\n\nimport librosa\nimport librosa.display\n\nimport os\n#for dirname, _, filenames in os.walk('/kaggle/input'):\n    #for filename in filenames:\n        #print(os.path.join(dirname, filename))\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n        \nfrom sklearn.decomposition import PCA\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nfrom scipy.ndimage import maximum_filter1d\nfrom scipy.ndimage import minimum_filter1d\n\nfrom datetime import datetime\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def reduce_mem(df):\n    for col in df.columns:\n        if df[col].dtype in ['float64','float32']:\n            df[col] = df[col].astype(np.float16)\n            \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"V_PATH = '/kaggle/input/predict-volcanic-eruptions-ingv-oe/'\nTRAIN_PATH = V_PATH + 'train/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"SENSOR_COLS = ['sensor_1', 'sensor_2', 'sensor_3', 'sensor_4', 'sensor_5', 'sensor_6',\n       'sensor_7', 'sensor_8', 'sensor_9', 'sensor_10']\n\n\nSENSOR_RMEANS = [x+'_rmin' for x in SENSOR_COLS] \nSENSOR_RSTDS = [x+'_rstd' for x in SENSOR_COLS] \nSENSOR_RMINS = [x+'_rmin' for x in SENSOR_COLS] \nSENSOR_RMAXES = [x+'_rmax' for x in SENSOR_COLS]\nSENSOR_RGRADMEAN = [x+'_grad_rmean' for x in SENSOR_COLS]\nSENSOR_RGRADSTD = [x+'_grad_rstd' for x in SENSOR_COLS]\n\nSENSOR_RSTATS = [SENSOR_RMEANS, SENSOR_RSTDS, SENSOR_RMINS, SENSOR_RMAXES,\n               SENSOR_RGRADMEAN,  SENSOR_RGRADSTD]\n\nROLL_DESCR = ['rmin', 'rstd', 'rmin', 'rmax', 'grad_rmean','grad_rstd']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"dtypes_dict = {'sensor_1': 'float32',\n 'sensor_2': 'float32',\n 'sensor_3': 'float32',\n 'sensor_4': 'float32',\n 'sensor_5': 'float32',\n 'sensor_6': 'float32',\n 'sensor_7': 'float32',\n 'sensor_8': 'float32',\n 'sensor_9': 'float32',\n 'sensor_10': 'float32'}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I don't have tons and tons of time, so this is not particularly neat or tidy, just trying to get some ideas of the data as quickly as possible"},{"metadata":{},"cell_type":"markdown","source":"# 1. Inspect train and test CSV\nThere are 4431 train data files, each with a row in train"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(V_PATH+'train.csv')\nprint(train.shape)\nprint(train.columns)\ntrain.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Individual data files are around 60000 rows each"},{"metadata":{"trusted":true},"cell_type":"code","source":"total_rows_estimate = 60001 * len(train) / 1000000\nprint('estimate of total TRAIN rows (millions)',total_rows_estimate)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv(V_PATH+'sample_submission.csv')\nprint(sample_submission.shape)\nprint(sample_submission.columns)\n\ntotal_rows_estimate = 60001 * len(sample_submission) / 1000000\nprint('estimate of total TEST rows (millions)',total_rows_estimate)\n\nsample_submission.head(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The rows of the train file show a pretty even distribution of times until next eruption"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(figsize=(10,5))\naxes.set_title('Distribution of Train time to eruption (scaled 0 to 1 = max)', size=16)\nsns.kdeplot(train['time_to_eruption'] / train['time_to_eruption'].max(), color='Red')\naxes.set_xlabel('Time To Eruption (Scaled 0 - 1)', size=12)\naxes.set_ylabel('Density of Train Data', size=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Minimum time to eruption', np.round(train['time_to_eruption'].min(),0)) \nprint('Max time to eruption', np.round(train['time_to_eruption'].max(),0))\nprint('Mean time to eruption', np.round(train['time_to_eruption'].mean(),0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#rescaling time to eruption to simplify looking at some of the numbers later\nsz = train['time_to_eruption'].size-1\ntrain['PCNT_TIME'] = np.round(train['time_to_eruption'].rank(method='max').apply(lambda x: 100.0*(x-1)/sz),0)\ntrain.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def get_rolling(df, cols, window=50):\n    for col in cols:\n        df[col+'_grad'] = np.gradient(df[col])\n        df[col+'_grad'] = df[col+'_grad'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_grad_abs'] = np.gradient(np.abs(df[col]))\n        df[col+'_grad'] = df[col+'_grad'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_rmin'] = minimum_filter1d(df[col].values, size=window)\n        df[col+'_rmax'] = maximum_filter1d(df[col].values, size=window)\n        \n        df[col+'_rmin'] = df[col+'_rmin'].fillna(method='bfill').fillna(method='ffill')\n        df[col+'_rmax'] = df[col+'_rmax'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_rmean'] = df[col].rolling(window=window, center=True).mean().fillna(method='bfill').fillna(method='ffill')\n        df[col+'_rstd'] = df[col].rolling(window=window, center=True).std().fillna(method='bfill').fillna(method='ffill')\n        \n        #add also for gradients\n        df[col+'_grad_rmin'] = minimum_filter1d(df[col+'_grad_abs'].values, size=window)\n        df[col+'_grad_rmax'] = maximum_filter1d(df[col+'_grad_abs'].values, size=window)\n        \n        df[col+'_grad_rmin'] = df[col+'_grad_rmin'].fillna(method='bfill').fillna(method='ffill')\n        df[col+'_grad_rmax'] = df[col+'_grad_rmax'].fillna(method='bfill').fillna(method='ffill')\n        \n        df[col+'_grad_rmean'] = df[col+'_grad_abs'].rolling(window=window, center=True).mean().fillna(method='bfill').fillna(method='ffill')\n        df[col+'_grad_rstd'] = df[col+'_grad_abs'].rolling(window=window, center=True).std().fillna(method='bfill').fillna(method='ffill')\n        \n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2. Explore some sample data\nSelect some example train data frames to run some initial analysis, without needing to import all of the data - use some quantiles to ensure a mix of durations to next eruption (simplified - only for EDA)\n\nCreate some rolling data analysis as they are imported (get rolling mins, maxes, stds for each sensor)\n\nNote: there are quite a few NAs / missing values in the data. For now these are just filled with zeros."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"select_quantiles = np.arange(0.01, 0.99, 0.02)\nprint('Quantiles list',select_quantiles)\n\nloaded_dfs = pd.DataFrame()\n\nfor count,q in enumerate(select_quantiles):\n    \n    s_ID = train['segment_id'][train['PCNT_TIME']==int(q*100)].values[0]\n    #print(s_ID)\n    \n    temp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv', dtype=dtypes_dict)\n                    #     na_values=0)\n    #print('count NAs Percent, ', temp_df.isna().sum().sum() / len(temp_df.values.flatten()))\n    #for col in temp_df.columns:\n    temp_df = temp_df.fillna(method='ffill').fillna(method='bfill')\n    temp_df = temp_df.fillna(value=0)\n    \n    temp_df = get_rolling(temp_df, SENSOR_COLS)\n    \n    \n    #print('count NAs Percent, ', temp_df.isna().sum().sum() / len(temp_df.values.flatten()))\n    temp_df['time_to_eruption'] = q\n    temp_df['segment'] = count\n    loaded_dfs = pd.concat([loaded_dfs, temp_df])\n    #print('    ')\nloaded_dfs = loaded_dfs.reset_index(drop=True) \nprint(' Total DF Size ',loaded_dfs.shape)\nloaded_dfs.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2a - Ranges of Sensor Data\nI realised that my earlier versions neglected any basic analysis of the range of values covered by each sensor in these example extracts\nLets see if the sensors report similar ranges in this sample data"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Sensor Max values (absolute)')\nnp.abs(loaded_dfs[SENSOR_COLS]).max()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Sensor Mean values (absolute)')\nnp.abs(loaded_dfs[SENSOR_COLS]).mean()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Sensor Median values (absolute)')\nnp.abs(loaded_dfs[SENSOR_COLS]).median()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see there are some significant differences in the max and median sensor readings.\nNext question - how do the distributions look if we plot them? I'm going to ignore zero values as these are 'missings'."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes = plt.subplots(figsize=(20,5))\n\nfor sc in SENSOR_COLS:\n    sns.kdeplot(np.abs(loaded_dfs[sc])[loaded_dfs[sc]!=0], ax=axes)\n    \naxes.set_xlim(0, 2500)\naxes.set_title('Absolute Sensor Reading Distributions (zero values dropped)', size=18)\n\naxes.set_xlabel('Absolute Sensor Reading (graph clipped at max = 2500)', size=12)\naxes.set_ylabel('Density of Sample Train Data', size=12)\naxes.legend(fontsize=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK, so we definitely have some significant differences in the spread of readings. \n\n# 2b. Correlations between sensor readings\n\nLet's check how well correlated the sensors are with one another."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sensor_corr = np.abs(loaded_dfs[SENSOR_COLS]).corr()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(figsize=(8,8))\n\nsns.heatmap(sensor_corr, annot=True, cbar=False, cmap='seismic_r',\n           vmin=0,vmax=1,ax=axes)\n\naxes.set_title('Sensor Signal (Absolute) Correlations', size=20)\naxes.yaxis.set_tick_params(labelsize=14)\naxes.xaxis.set_tick_params(labelsize=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Correlations between the signals seem to be modest. Are the correlations stronger closer to an eruption?"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(ncols=2,figsize=(20,8))\n\nsensor_corr = np.abs(loaded_dfs[SENSOR_COLS][loaded_dfs['time_to_eruption']<0.5]).corr()\n\nsns.heatmap(sensor_corr, annot=True, cbar=False, cmap='seismic_r',\n           vmin=0,vmax=1,ax=axes[0])\n\naxes[0].set_title('Sensor Signal (Absolute) Correlations Closer to Eruption', size=20)\n\nsensor_corr = np.abs(loaded_dfs[SENSOR_COLS][loaded_dfs['time_to_eruption']>=0.5]).corr()\n\nsns.heatmap(sensor_corr, annot=True, cbar=False, cmap='seismic_r',\n           vmin=0,vmax=1,ax=axes[1])\n\naxes[1].set_title('Sensor Signal (Absolute) Correlations Further From Eruption', size=20)\n\n\naxes[0].yaxis.set_tick_params(labelsize=12)\naxes[0].xaxis.set_tick_params(labelsize=12)\n\naxes[1].yaxis.set_tick_params(labelsize=12)\naxes[1].xaxis.set_tick_params(labelsize=12)\n\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Interesting. It looks like the correlation between signals is extremely poor when we are far away from an eruption, at least for this sample data.\nIt also looks to me like there are some plausible 'groupings' of sensors"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(ncols=1,figsize=(8,8))\n\nsensor_corr = np.abs(loaded_dfs[SENSOR_COLS][loaded_dfs['time_to_eruption']<0.5]).corr()\n\ngrouped_list = ['sensor_1','sensor_2','sensor_3','sensor_4','sensor_9',\n               'sensor_5','sensor_10',\n             'sensor_6','sensor_8','sensor_7'  ]\n\nsensor_corr = sensor_corr[grouped_list].T[grouped_list].T\n\nsns.heatmap(sensor_corr, annot=True, cbar=False, cmap='seismic_r',\n           vmin=0,vmax=1,ax=axes)\n\naxes.set_title('Sensor Signal (Absolute) Correlations - Reordered', size=20)\naxes.yaxis.set_tick_params(labelsize=14)\naxes.xaxis.set_tick_params(labelsize=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This kind of looks like 'blocks' to me. \n* Sensors 1-2-3-4-9\n* Sensors 5-10\n* Sensors 6-8-7\n\nThe last 3 only seem to share the characteristic that they are not really very well correlated with any of the others.\nWe'll have to see whether this can provide any useful outcome with predictions"},{"metadata":{},"cell_type":"markdown","source":"Finally in this section, lets check how often each sensor has the largest absolute value."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('Percentage of Rows where the sensor is the highest absolute value')\n\nloaded_dfs['strongest_signal']  = np.argmax(np.abs(loaded_dfs[SENSOR_COLS].values), axis=1) + 1\nloaded_dfs['strongest_signal'].value_counts() / loaded_dfs['strongest_signal'].value_counts().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sensors 10, 6 and 2 have the highest proportion of rows where the are the largest absolute value recorded. \nHow about any correlation with eruptions?"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Average Time to Eruption grouped by strongest signal at point in time')\nprint('Note: time to eruption is scaled 0-1.0 for easier viewing')\nloaded_dfs.groupby(['strongest_signal'])['time_to_eruption'].mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The best indicator from this data seems to be when signal 2 is the strongest. But overall it doesn't look all that helpful - many of the values are around the 0.5 mark"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"So we know there is some correlation between sensor readings, I guess this makes sense as we could guess they would be picking up the same events, to some extent. Let's compare how they record during some data which shows what look like some strong signals."},{"metadata":{},"cell_type":"markdown","source":"Note: this is not an attempt to start trying to identify actual real sensor locations etc, I'm assuming that comparing and analysing readings from multiple sensors is a generally applicable approach in this field."},{"metadata":{"trusted":true},"cell_type":"code","source":"s_ID=1424510231 \n\n#print('Segment ID' ,s_ID)\n\nprint(train[train['segment_id']==s_ID])\n\nfig,axes=plt.subplots(nrows=1,ncols=1,figsize=(20,6))\ntemp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv')\n\nfor s in SENSOR_COLS:\n    axes.plot(temp_df[s], linewidth=1)\n    \naxes.set_title('Readings from each Sensor for segment ID '+str(s_ID), size=18)\naxes.legend(SENSOR_COLS, fontsize=12)\n\n\naxes.set_xlabel('Time (total duration = one 10 minute segment)', size=14)\naxes.set_ylabel('Absolute Sensor Reading', size=14)\n#axes.legend(fontsize=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What does the correlation look like for this segment of data?"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(figsize=(8,8))\n\nsns.heatmap(np.abs(temp_df[SENSOR_COLS]).corr()[grouped_list].T[grouped_list], annot=True, cbar=False, cmap='seismic_r',\n           vmin=0,vmax=1,ax=axes)\n\naxes.set_title('Sensor Signal (Absolute) Correlations for Segment '+str(s_ID), size=20)\n\naxes.yaxis.set_tick_params(labelsize=14)\naxes.xaxis.set_tick_params(labelsize=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets look at 3 comparisons: \n* Sensors 1 and 3 show high correlation\n* Sensors 4 and 6 show moderate correlation\n* Sensors 7 and 8 are fairly well correlated between themselves, but less so with other sensors."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=1,ncols=1,figsize=(20,6))\ntemp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv')\n\nplot_sensors=['sensor_1', 'sensor_3']\n\nfor s in plot_sensors:\n    axes.plot(temp_df[s], linewidth=1)\n    \naxes.set_title('Selected Sensor Readings for segment ID '+str(s_ID), size=18)\naxes.legend(plot_sensors, fontsize=12)\n\n\naxes.set_xlabel('Time (total duration = one 10 minute segment)', size=14)\naxes.set_ylabel('Absolute Sensor Reading', size=14)\n#axes.legend(fontsize=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=1,ncols=1,figsize=(20,6))\ntemp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv')\n\nplot_sensors=['sensor_4', 'sensor_6']\n\nfor s in plot_sensors:\n    axes.plot(temp_df[s], linewidth=1)\n    \naxes.set_title('Selected Sensor Readings for segment ID '+str(s_ID), size=18)\naxes.legend(plot_sensors, fontsize=12)\n\naxes.set_xlabel('Time (total duration = one 10 minute segment)', size=14)\naxes.set_ylabel('Absolute Sensor Reading', size=14)\n#axes.legend(fontsize=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=1,ncols=1,figsize=(20,6))\ntemp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv')\n\nplot_sensors=['sensor_7', 'sensor_8', 'sensor_1']\n\nfor s in plot_sensors:\n    axes.plot(temp_df[s], linewidth=1)\n    \naxes.set_title('Selected Sensor Readings for segment ID '+str(s_ID), size=18)\naxes.legend(plot_sensors, fontsize=12)\n\naxes.set_xlabel('Time (total duration = one 10 minute segment)', size=14)\naxes.set_ylabel('Absolute Sensor Reading', size=14)\n#axes.legend(fontsize=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Distributions of data / gradients by time to next eruption"},{"metadata":{},"cell_type":"markdown","source":"\nRed = Closer to Eruption, Green = further away"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"SEGMENTS = pd.Series(loaded_dfs['time_to_eruption'])\nnr = len(SENSOR_COLS)\nnc = 4\ntg = nr * nc\n\nfig,axes=plt.subplots(nrows=nr, ncols=nc,figsize=(20,nr*4))\n\n\ncolor_grades = ['Red', 'darkorange','wheat','yellowgreen','seagreen']\n\nfor count, SC in enumerate(SENSOR_COLS):\n    for Q in np.flip(select_quantiles):\n        sns.kdeplot(loaded_dfs[SC][loaded_dfs['time_to_eruption']==Q], ax=axes[count, 0],\n                   alpha=0.7, color=color_grades[int(0.5+Q*4)])\n        \n        sns.kdeplot(loaded_dfs[SC+'_grad_rmax'][loaded_dfs['time_to_eruption']==Q], ax=axes[count, 1],\n                   alpha=0.7, color=color_grades[int(0.5+Q*4)], bw=0.1)\n        \n        sns.kdeplot(loaded_dfs[SC+'_grad_rstd'][loaded_dfs['time_to_eruption']==Q], ax=axes[count, 2],\n                   alpha=0.7, color=color_grades[int(0.5+Q*4)], bw=0.1)\n        \n        sns.kdeplot(loaded_dfs[SC+'_grad_rmin'][loaded_dfs['time_to_eruption']==Q], ax=axes[count, 3],\n                   alpha=0.7, color=color_grades[int(0.5+Q*4)], bw=0.1)       \n        \n        \n        \n    xulim = loaded_dfs[SC].quantile(0.95)\n    xllim = loaded_dfs[SC].quantile(0.05)\n    #axes[count, 0].legend(None)\n    axes[count, 0].set_xlim(xllim, xulim)\n    axes[count, 0].get_legend().remove()    \n    axes[count, 0].set_title(SC+' Signal Values')\n    \n    xulim = loaded_dfs[SC+'_grad_rmax'].quantile(0.95)\n    xllim = loaded_dfs[SC+'_grad_rmax'].quantile(0.05)\n    axes[count, 1].set_xlim(xllim, xulim)\n    axes[count, 1].get_legend().remove()    \n    axes[count, 1].set_title(SC+' Grads RollMax')\n    \n    xulim = loaded_dfs[SC+'_grad_rstd'].quantile(0.95)\n    xllim = loaded_dfs[SC+'_grad_rstd'].quantile(0.05)\n    axes[count, 2].set_xlim(xllim, xulim)\n    axes[count, 2].get_legend().remove()    \n    axes[count, 2].set_title(SC+' Grads RollSTD')\n    \n    xulim = loaded_dfs[SC+'_grad_rmin'].quantile(0.95)\n    xllim = loaded_dfs[SC+'_grad_rmin'].quantile(0.05)\n    axes[count, 3].set_xlim(xllim, xulim)\n    axes[count, 3].get_legend().remove()    \n    axes[count, 3].set_title(SC+' Grads RollMin')\n    \nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The data looks quite scattered, and the sensors don't look that consistent. Sensor 5 readings look a bit odd in these data samples. It looks like there are some specific high concentrations of readings in a narrow range."},{"metadata":{},"cell_type":"markdown","source":"# 4. Check on some sample data for NA / missing values\n\nThere are quite a lot of missing values. Let's check how many / where these are. Just looking at the first 1000 data files for train and test to save some time and RAM"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def get_nas(path, file_id):\n    temp_df = pd.read_csv(path+str(file_id)+'.csv', dtype=dtypes_dict)\n    return temp_df.isna().sum().values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\nSENSOR_COLS_MISSING = [x+'_missing' for x in SENSOR_COLS]\nfor SC in SENSOR_COLS_MISSING:\n    train[SC] = 0.0\n    sample_submission[SC] = 0.0\n\nfor count,i in enumerate(train.index[0:1000]):\n    #print(count/len(train))\n    s_ID = train.loc[i,'segment_id']    \n    #temp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv')\n    train.loc[i,SENSOR_COLS_MISSING] = get_nas(TRAIN_PATH,s_ID)\n    #del temp_df\n    #gc.collect()\nprint('finished Train load')\n    \nTEST_PATH = V_PATH + 'test/'  \nfor count,i in enumerate(sample_submission.index[0:1000]):\n    #print(count/len(sample_submission))\n    s_ID = sample_submission.loc[i,'segment_id']    \n    sample_submission.loc[i,SENSOR_COLS_MISSING] = get_nas(TEST_PATH,s_ID)\n    \ngc.collect()\nprint('finished Test load')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sample_submission[SENSOR_COLS_MISSING] = sample_submission[SENSOR_COLS_MISSING] / sample_submission[SENSOR_COLS_MISSING].max().max()\ntrain[SENSOR_COLS_MISSING] = train[SENSOR_COLS_MISSING] / train[SENSOR_COLS_MISSING].max().max()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"train_missings = pd.DataFrame(columns=['mean','pc dfs with missing'],\n                             index=SENSOR_COLS_MISSING)\ntrain_missings['mean'] = train[SENSOR_COLS_MISSING].mean().values\nfor sc in SENSOR_COLS_MISSING:\n    train_missings.loc[sc, 'pc dfs with missing'] = len(train[train[sc]>0]) / len(train)\n#train_missings['count with missing'] = train[SENSOR_COLS_MISSING].count().va","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"test_missings = pd.DataFrame(columns=['mean','pc dfs with missing'],\n                             index=SENSOR_COLS_MISSING)\ntest_missings['mean'] = sample_submission[SENSOR_COLS_MISSING].mean().values\nfor sc in SENSOR_COLS_MISSING:\n    test_missings.loc[sc, 'pc dfs with missing'] = len(sample_submission[sample_submission[sc]>0]) / len(sample_submission)\n#train_missings['count with missing'] = train[SENSOR_COLS_MISSING].count().values","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false},"cell_type":"markdown","source":"It looks like distributions of missing values are different in Test data. we have a lot more missing values for sensors 9 and 10 in particular, it seems"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig,axes=plt.subplots(ncols=2,figsize=(15,6), sharey=True)\n\nsns.heatmap(train_missings.fillna(value=0), annot=True, cbar=False, cmap='seismic',fmt='.1%', annot_kws={\"fontsize\":12},\n           vmin=0,vmax=0.05,ax=axes[0])\n\nsns.heatmap(test_missings.fillna(value=0), annot=True, cbar=False, cmap='seismic',fmt='.1%', annot_kws={\"fontsize\":12},\n           vmin=0,vmax=0.05, ax=axes[1])\n\naxes[0].set_title('Sample Train DFs - Missings', size=18)\naxes[1].set_title('Sample Test DFs - Missings', size=18)\n\naxes[0].yaxis.set_tick_params(labelsize=14)\naxes[0].xaxis.set_tick_params(labelsize=14)\n\naxes[1].yaxis.set_tick_params(labelsize=14)\naxes[1].xaxis.set_tick_params(labelsize=14)\n\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"for col in SENSOR_COLS_MISSING:\n    train = train.drop(col, axis=1)\n    sample_submission = sample_submission.drop(col, axis=1)\n    \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# 5. Add in some analysis across sensors\nLet's take some mins, maxes, means etc across rows (of original sensor data, and of the groups of rolling columns e.g. the max, min sensor rolling gradients) to see what information combining sensor readings can give us."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def get_stats(df, sensor_cols=SENSOR_COLS, rolling_cols=SENSOR_RSTATS):\n    #we create the min max etc of original sensor columns\n    df['max'] = df[sensor_cols].max(axis=1).astype(np.float16)\n    df['min'] = df[sensor_cols].min(axis=1).astype(np.float16)\n    df['std'] = df[sensor_cols].std(axis=1).astype(np.float16)\n    \n    #and with absolute values\n    df['max_abs'] = np.abs(df[sensor_cols]).max(axis=1).astype(np.float16)\n    df['min_abs'] = np.abs(df[sensor_cols]).min(axis=1).astype(np.float16)\n    df['std_abs'] = np.abs(df[sensor_cols]).std(axis=1).astype(np.float16)\n    \n    #we take mins and maxes of groups of rolling columns\n    for count,rc in enumerate(rolling_cols): #this takes a SINGLE mean, max across each GROUP of rolling\n        #columns - e.g. the max of all rolling mins\n        df[ROLL_DESCR[count]+'_max'] = df[rolling_cols[count]].max(axis=1).astype(np.float16)\n        df[ROLL_DESCR[count]+'_min'] = df[rolling_cols[count]].min(axis=1).astype(np.float16)\n        df[ROLL_DESCR[count]+'_std'] = df[rolling_cols[count]].std(axis=1).astype(np.float16)\n        df[ROLL_DESCR[count]+'_mean'] = df[rolling_cols[count]].mean(axis=1)  .astype(np.float16)  \n   \n    gc.collect()\n    return df\n\nloaded_dfs = get_stats(loaded_dfs, sensor_cols=SENSOR_COLS, rolling_cols=SENSOR_RSTATS)\nloaded_dfs = reduce_mem(loaded_dfs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ldf_corr = loaded_dfs.corr()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot some correlations with our target column"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ldf_corr = ldf_corr.sort_values('time_to_eruption', ascending=False)\n\nldf_corr = ldf_corr[['time_to_eruption'] + [ col for col in ldf_corr.columns if col != 'time_to_eruption' ]]\n\nfilt = np.abs(ldf_corr['time_to_eruption'])>0.05\n\nfig,axes=plt.subplots(figsize=(20,20))\nsns.heatmap(ldf_corr[filt], annot=False, cmap='seismic', vmin=-1, vmax=1, cbar=False)\n\ngc.collect()\n\naxes.set_title('Correlation with Time to Eruption (all features with >5% correlation)', size=18)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Save some of the features which are >5% correlated for later use"},{"metadata":{"trusted":true},"cell_type":"code","source":"TOP_FTS = ldf_corr[filt].index\nprint(len(TOP_FTS))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Rolling features definitely seem useful (no surprise - as important variations in the data would likely last longer than 1 row, but less than ten minutes)"},{"metadata":{},"cell_type":"markdown","source":"As a side note:\n\nhttps://www.nature.com/articles/s41467-020-17375-2\n\n\"We processed 9 years of data, from 1 January 2011 to 1 January 2020 (Fig. 1b), to obtain four time series that capture different parts of the tremor signal (RSAM, DSAR, medium and high frequency bands—MF and HF—sampled every 10 min; see “Methods”)\"\n\nInteresting that this also looked at 10 minutes of data as being a good time frame for detection.    "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ldf_corr = ldf_corr.sort_values('time_to_eruption', ascending=False)\n\nfilt = np.abs(ldf_corr['time_to_eruption'])>0.3\n\nfig,axes=plt.subplots(figsize=(10,7))\nsns.heatmap(ldf_corr[['time_to_eruption', 'segment']][filt], \n            annot=True, cmap='seismic', vmin=-1, vmax=1, cbar=False)\n\ngc.collect()\n\naxes.set_title('Correlation with Time to Eruption (all features with >30% correlation)')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems that Gradient Rolling Standard Deviation min() has a negative correlation with time to eruption - i.e. when it is higher, eruptions are closer.\n# 6. Look at some segments (10 min batches)"},{"metadata":{},"cell_type":"markdown","source":"Let's look at how the gradient rolling std min is distributed across some sub segment (1000 rows)"},{"metadata":{"trusted":true},"cell_type":"code","source":"#create subsegments of 1000 rows\nloaded_dfs['sub_segment'] = loaded_dfs.index // 1000\nsegment_means = loaded_dfs.groupby(['sub_segment'])['grad_rstd_min'].mean().sort_values(ascending=False)\n\n#lets ignore any where this value is zero (possibly a filled NA or other error)\nsegment_means = segment_means[segment_means>0]\nfig,axes=plt.subplots(figsize=(15,5))\nsns.kdeplot(segment_means.values,ax=axes, color='Red')\naxes.set_title('Sub Segment Gradient Rolling STD Mins', size=18)\n\naxes.set_xlabel('Gradient Rolling STD Mins', size=12)\naxes.set_ylabel('Density of Sample Train Data', size=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Segments with higher Gradients')\nprint(segment_means.head(10))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Segments with medium Gradients')\nprint(segment_means[len(segment_means)//2-5:len(segment_means)//2+5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Segments with lower Gradients')\nprint(segment_means.tail(10))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Lets check out sensor readings from 3 contrasting segments"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"SEGMENT_HIGH = 179\nSEGMENT_MEDIUM = 445\nSEGMENT_LOW = 1609\n\nnr = len(SENSOR_COLS)\nnc = 2\ntg = nr * nc\n\nTITLES = ['Close to Eruption','Mid Point', 'Far from Eruption']\ncolors= ['Red','Orange','Green']\n\nfig,axes=plt.subplots(nrows=nr, ncols=nc,figsize=(20,nr*4))\n\nfor count1, SEG in enumerate([SEGMENT_HIGH, SEGMENT_MEDIUM,SEGMENT_LOW]):\n    \n    seg_filt = loaded_dfs['sub_segment']==SEG\n    \n    for count, SC in enumerate(SENSOR_COLS):\n\n        sns.lineplot(x=range(sum(seg_filt)),y=loaded_dfs[SC][seg_filt], ax=axes[count, 0],\n                        color=colors[count1])\n\n        sns.lineplot(x=range(sum(seg_filt)), y=loaded_dfs[SC+'_grad_rstd'][seg_filt], ax=axes[count, 1],\n                     color=colors[count1])\n\n        axes[count, 0].legend(TITLES, fontsize=12)\n       \n        axes[count, 0].set_title(SC + ' Readings', size=18)\n\n        axes[count, 1].set_title(SC+' Gradients Rolling Max', size=18)\n        axes[count, 1].legend(TITLES, fontsize=12)\n    \nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**And the minimum across sensors\nNote: the green line probably includes an error (missing sensor data) - will need to tackle NA values**"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=1, ncols=1,figsize=(12,7))\nfor count1, SEG in enumerate([SEGMENT_HIGH, SEGMENT_MEDIUM,SEGMENT_LOW]):\n    \n    seg_filt = loaded_dfs['sub_segment']==SEG\n    sns.lineplot(x=range(sum(seg_filt)), y=loaded_dfs['grad_rstd_min'][seg_filt], ax=axes,\n                     color=colors[count1])\n    \n    axes.legend(TITLES, fontsize=12)\n    \naxes.set_title('Gradient Rolling STD Mins Comparison', size=18)\naxes.set_xlabel('Time (one subsegment, 1000 rows)', size=14)\naxes.set_ylabel('Gradient Rolling STD Mins', size=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**The mean, to mitigate the issue with missing sensor data**"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=1, ncols=1,figsize=(12,7))\n\nfor count1, SEG in enumerate([SEGMENT_HIGH, SEGMENT_MEDIUM,SEGMENT_LOW]):\n    \n    seg_filt = loaded_dfs['sub_segment']==SEG\n    sns.lineplot(x=range(sum(seg_filt)), y=loaded_dfs['grad_rstd_mean'][seg_filt], ax=axes,\n                     color=colors[count1])\n    \n    axes.legend(TITLES, fontsize=12)\n    \naxes.set_title('Gradient Rolling STD Means Comparison', size=18)\naxes.set_xlabel('Time (one subsegment, 1000 rows)', size=14)\naxes.set_ylabel('Gradient Rolling STD Means', size=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 7. Spectrograms of sensor data"},{"metadata":{},"cell_type":"markdown","source":"As I wasn't really sure the best way to generate the spectrogram, I've borrowed some code from here:\nhttps://www.kaggle.com/michael422/spectrogram-convolution/notebook"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from scipy.signal import spectrogram\n\ndef return_spectrogram(sig_in, dsamp):\n    nperseg = 256 # default 256\n    noverlap = nperseg // 4 # default: nperseg // 8\n    fs = 4000000 // dsamp # raw signal sample rate is 4MHz\n    window = 'triang'\n    scaling = 'density' # {'density', 'spectrum'}\n    detrend = 'linear' # {'linear', 'constant', False}\n    eps = 1e-11\n    f, t, Sxx = spectrogram(sig_in, nperseg=nperseg, noverlap=noverlap,\n                                   fs=fs, window=window,\n                                   scaling=scaling, detrend=detrend)\n    return f, t, np.log(Sxx + eps)\n\ndef plot_segment(segment):\n    fig,axes=plt.subplots(ncols=2,nrows=len(SENSOR_COLS), figsize=(20,2.5*len(SENSOR_COLS)))\n    time_to_eruption = loaded_dfs['time_to_eruption'][loaded_dfs['segment']==segment].mean() * train['time_to_eruption'].max()\n    print('Spectrograms at Time to Eruption: ', time_to_eruption)\n    for count,s in enumerate(SENSOR_COLS):\n        sensor_data = loaded_dfs[s][loaded_dfs['segment']==segment].copy()\n        #print(sensor_data.shape)\n        axes[count,1].imshow(return_spectrogram(sensor_data.values, 100)[2],\n                   cmap='seismic')\n        axes[count,1].set_title(s + ' Spectrogram', size=16)\n        \n        axes[count,0].plot(sensor_data, color='Red')\n        axes[count,0].set_title(s + ' Signal', size=16)\n        \n    plt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_segment(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_segment(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_segment(40)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# 8. Check out some possible features\nAs it's going to be more challenging to load all train data (millions of rows), lets see what statistcs from each 10-minute section of data are most useful first"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#credit - stack overflow\n\ndef percentile(n):\n    def percentile_(x):\n        return np.percentile(x, n)\n    percentile_.__name__ = 'percentile_%s' % n\n    return percentile_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gc.collect()\n\nsegment_summaries = loaded_dfs.groupby(['segment'])[TOP_FTS].agg(['max','skew','std','min','mean',\n                                    percentile(0.1),percentile(0.25),\n                                        percentile(0.75),percentile(0.9)])\n\nsegment_summaries.columns = [a+b for a,b in segment_summaries.columns]\nsegment_summaries_corr = segment_summaries.corr()\nsegment_summaries_corr = segment_summaries_corr.sort_values('time_to_eruptionmean')\n\ndrop_fts = ['eruption', 'segment']\nFTS=segment_summaries_corr.columns\nfor d in drop_fts:\n    FTS = [x for x in FTS if d not in x]\n\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets look at some of the correlations with our target time to eruption. We don't have all that much train data loaded for the exploration, so these correlations are probably quite dependent on the specific sections of data we happen to be looking at. But hopefully they will give some ideas."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"filt = (np.abs(segment_summaries_corr['time_to_eruptionmean'])>0.3) & \\\nsegment_summaries_corr.index.isin(FTS)\n\nfig,axes=plt.subplots(figsize=(20,20))\nsns.heatmap(segment_summaries_corr[filt], annot=False, cmap='seismic', vmin=-1, vmax=1, cbar=False)\n\ngc.collect()\n\naxes.set_title('Segment Level Correlation with Time to Eruption (all features with >30% correlation)', size=20)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Trying to reduce our feature list a bit further for next step to reduce calculation time"},{"metadata":{},"cell_type":"markdown","source":"I think we're seeing based on this set of data that quite a lot of the most positive correlated data is rolling mins and gradients, and a lot of the negatively correlated (closer to eruption) features are percentiles. This is a small subset of the data though, so we'll have to explore further.\n"},{"metadata":{},"cell_type":"markdown","source":"Lets try looking at subsegments of the data. This could maybe also help to think about what time frames / windows are useful to apply.\nTo do this, I'm going to try to group by varying 'window' lengths, and look at some average / max correlations of the groupby summary data. I'm just trying some large window sizes, because much smaller sizes take too long (maybe I can think of ways to speed up...). Based on this analysis, with this quite simple approach, it doesn't seem to help to group summary stats into a window of 1000 rather than 10000 or 30000. Doubtless more sophisticated approaches looking at specific events in the sensor data would do much better on shorter batches.\n\nNote: again think we have to keep in mind that this is not a very perfect analysis, with limited train data sets."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"windows_test = [1000, 5000, 10000, 30000]\nwindows_test_results = pd.Series(index=windows_test)\nfig,axes=plt.subplots(figsize=(20,6))\nfor wt in windows_test:\n    print('checking window size', wt)\n    loaded_dfs['sub_segment'] = loaded_dfs.index // wt\n\n    segment_summaries = loaded_dfs.groupby(['sub_segment'])[TOP_FTS].agg(['mean','max','skew','std','min',\n                                    percentile(0.1),percentile(0.25),\n                                        percentile(0.75),percentile(0.9)])\n    \n    segment_summaries.columns = [a+b for a,b in segment_summaries.columns]\n    segment_summaries_corr = segment_summaries.corr()\n    segment_summaries_corr = segment_summaries_corr.fillna(value=0)\n    segment_summaries_corr = segment_summaries_corr.sort_values('time_to_eruptionmean')\n    \n    print('Mean Correlation',np.abs(segment_summaries_corr.loc[FTS,'time_to_eruptionmean']).mean())\n    windows_test_results[wt] = np.abs(segment_summaries_corr.loc[FTS,'time_to_eruptionmean']).mean()\n    \n    sns.kdeplot(np.abs(segment_summaries_corr.loc[FTS,'time_to_eruptionmean']))\n    \naxes.legend(windows_test, fontsize=16)\naxes.set_title('Distribution of feature correlations by Window size (note: negative correlations plotted as absolutes)', size=18)\n\n\naxes.set_xlabel('Feature Correlation with proximity of eruption', size=12)\naxes.set_ylabel('Density of Features', size=12)\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 9. Analyse Gradient feature(s)\nNote: the sections of the data are roughly ordered with the closest to eruption being on the left - but the data is not consecutive..so this is indicative, rather than a thorough analysis."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"loaded_dfs['sub_segment'] = loaded_dfs.index // 100\n\nsegment_summaries = loaded_dfs.groupby(['sub_segment'])['grad_rstd_min',].agg(['mean','max','skew','std','min',\n                                    percentile(0.1),percentile(0.25),\n                                        percentile(0.75),percentile(0.9)])\n\nsegment_summaries.columns = [b for a,b in segment_summaries.columns]\n\nfig,axes=plt.subplots(nrows=segment_summaries.shape[1],\n                     figsize=(20,segment_summaries.shape[1]*7))\n\nfor count,c in enumerate(segment_summaries.columns):\n    axes[count].plot(segment_summaries.loc[:,c].values, color='Red')\n    axes[count].set_title('Grad. Rolling Std ' + c+' Left Side = Closer To Eruption', size=18)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have some quite high peaks on the far right (far from next eruption).\nI dont have any domain knowledge, but could these be tremors shortly after an eruption? Googling the subject seems to suggest that aftershocks could be possible. This could be a source of confusion for our modelling.\nMore in line with what might be expected, there are some increasingly large peaks on the left (close to an eruption)"},{"metadata":{},"cell_type":"markdown","source":"Lets see if its possible to cluster some of the raw data based on some top features. This does seem to reveal some smaller clusters closer to eruption - probably some of the 'events' visible in the graphs above."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"CLUSTER_FTS = ['min','sensor_2_grad_rmin',\n              'sensor_5_grad_rmin',\n              'sensor_3_grad_rmin',\n              'grad_rstd_min',\n              'rstd_min','rstd_mean','min_abs']\ngc.collect()\nkm = KMeans(n_clusters=8, random_state=42)\nStSc = StandardScaler()\nkm.fit(StSc.fit_transform(loaded_dfs[CLUSTER_FTS]))\ngc.collect()\nloaded_dfs['km_labels'] = km.labels_.astype(str)\n\n#print()\n#print(loaded_dfs['km_labels'].value_counts())\n\nloaded_dfs.groupby(['km_labels'])['time_to_eruption'].agg(['mean','count'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(figsize=(12,8))\nkm_summ = loaded_dfs.groupby(['km_labels',\n        'time_to_eruption'])['segment'].count().unstack('km_labels').fillna(value=0)\n\norder = loaded_dfs.groupby(['km_labels'])['time_to_eruption'].agg(['mean','count']).sort_values('mean').index\nkm_summ.index = np.round(km_summ.index,2)\nsns.heatmap(km_summ[order],\n           cmap='seismic', annot=False,vmin=0,cbar=False)\n\naxes.set_title('Clustered & Sorted by mean label proximity to eruption', size=18)\n\naxes.set_xlabel('KMeans Label', size=14)\naxes.set_ylabel('Proximity to Eruption (lower = closer)', size=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at the range of sensor readings for each of these 'clusters' - we'd probably guess that the clusters close to eruption will show a wider spread\nDarker lines indicate closer to eruption"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes=plt.subplots(nrows=len(SENSOR_COLS),\n                     figsize=(15,5*len(SENSOR_COLS)),\n                     sharex=True)\n\nlabel_prox = loaded_dfs.groupby(['km_labels'])['time_to_eruption'].mean()\n\nfor count,s in enumerate(SENSOR_COLS):\n    \n    for ksl in loaded_dfs['km_labels'].unique():\n        #print(ksl)\n        sns.kdeplot(loaded_dfs[s][loaded_dfs['km_labels']==ksl].values,ax=axes[count],\n                   color=color_grades[int(0.5+label_prox[ksl]*4)], alpha=1)\n    axes[count].set_xlim(-5000,5000)\n    axes[count].set_title(s + ' Distribution of Readings by Cluster Label - Red = Closer To Eruption', size=16)\n    axes[count].legend(order)\n        \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Some further analysis of sub segments to see what correlates well with approach eruptions"},{"metadata":{"trusted":true},"cell_type":"code","source":"segment_stats = loaded_dfs.groupby(['sub_segment'])[loaded_dfs.columns].mean()\n\ncorrelations=pd.Series(index=segment_stats.columns,\n                      data=0.0)\n\nfor i in correlations.index:\n    #print(i)\n    correlations[i] = np.corrcoef(segment_stats['time_to_eruption'],\n                                 segment_stats[i])[0,1]\n    \ncorrelations.sort_values(ascending=False).head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"correlations.sort_values(ascending=False).tail(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 10. Train vs Test Sensor Ranges\nThe predictions so far have shown a different shape for train and test. Lets load some randomised data"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"TEST_PATH = V_PATH + 'test/'\ntest_examples=pd.Series(sample_submission.index).sample(n=100,random_state=42).values\nprint('Examples list ',test_examples)\n\nloaded_test_dfs = pd.DataFrame()\n\nfor count,q in enumerate(test_examples):\n    \n    s_ID = sample_submission.loc[q,'segment_id']\n        \n    temp_df = pd.read_csv(TEST_PATH+str(s_ID)+'.csv', dtype=dtypes_dict)\n    \n    temp_df = temp_df.fillna(method='ffill').fillna(method='bfill')\n    temp_df = temp_df.fillna(value=0)\n    \n    temp_df['time_to_eruption'] = q\n    temp_df['segment'] = count\n    loaded_test_dfs = pd.concat([loaded_test_dfs, temp_df])\n    \n    gc.collect()\n    \nloaded_test_dfs = loaded_test_dfs.reset_index(drop=True) \nprint(' Total DF Size ',loaded_test_dfs.shape)\nloaded_test_dfs.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"train_examples=pd.Series(train.index).sample(n=100,random_state=42).values\nprint('Examples list',test_examples)\n\nloaded_train_dfs = pd.DataFrame()\n\nfor count,q in enumerate(train_examples):\n    \n    s_ID = train.loc[q,'segment_id']\n        \n    temp_df = pd.read_csv(TRAIN_PATH+str(s_ID)+'.csv', dtype=dtypes_dict)\n    \n    temp_df = temp_df.fillna(method='ffill').fillna(method='bfill')\n    temp_df = temp_df.fillna(value=0)\n    \n    temp_df['time_to_eruption'] = q\n    temp_df['segment'] = count\n    loaded_train_dfs = pd.concat([loaded_train_dfs, temp_df])\n    \n    gc.collect()\n    \nloaded_train_dfs = loaded_train_dfs.reset_index(drop=True) \nprint(' Total DF Size ',loaded_train_dfs.shape)\nloaded_train_dfs.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Examine some sample distributions of train vs test sensor data\nHiding readings of 0 or readings which are very extreme, to better view the central distributions"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"nc=2\nfig,axes=plt.subplots(nrows=10//nc,ncols=nc,figsize=(20,18))\n#fil\nfor count,S in enumerate(SENSOR_COLS):\n    \n    sns.kdeplot(loaded_train_dfs[S][(loaded_train_dfs[S]!=0) & (np.abs(loaded_train_dfs[S])<6000)],color='Green',\n               ax=axes[count//nc,count%nc])\n    \n    sns.kdeplot(loaded_test_dfs[S][(loaded_test_dfs[S]!=0) & (np.abs(loaded_test_dfs[S])<6000)],color='Red',\n               ax=axes[count//nc,count%nc])    \n    \n    axes[count//nc,count%nc].set_xlim(-5000,5000)\n    axes[count//nc,count%nc].legend(['Sample Train data', 'Sample Test data'], fontsize=12)\n    axes[count//nc,count%nc].set_title(S, size=18)\n    axes[count//nc,count%nc].set_xlabel('Sensor Reading', size=12)\n    \nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Within these randomly selected examples, it looks like the distributions are pretty similar, though the test data does look more clustered around the centre for 8 of the 10 sensors. This could imply that the average test time to the next eruption is a bit higher, if it is true of the whole Test data set.\nThe NAs were filtered out though, and the analysis further up highlighted that some sensors could have a lot more missing values in Test data."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}