{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Import and Explore Sample Training File\n\nWe'll take a look at a single `segment_id` to review a file example. At the same time importing the `train.csv` file."},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nfrom scipy import signal\nimport seaborn as sns\nimport glob","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_labels = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train.csv\")\ntrain_labels.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Total segment files: {}\".format(len(train_labels['segment_id'])))\ntrain_labels.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_example = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/train/1003154738.csv\")\ndf_example.head()\n\ndata_columns = list(df_example.columns)\n\nprint('Index Dataframe Shape: {}'.format(df_example.shape))\nprint('Column Headers:\\n')\nprint(data_columns)\ndf_example.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"10 sensors with 60001 readings. We also see a significant amount of NaN in there as well. Let's plot the time series data."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(nrows=5, ncols=2)\nfig.set_size_inches(20,10)\nfig.subplots_adjust(hspace=0.5)\n\nfor col,ax in zip(data_columns, axs.flatten()):\n    ax.plot(range(len(df_example[col])),df_example[col])\n    ax.set_title(col)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are going to build a feature set for each `segment_id` using the spectral density. This approach evaluates the magnitude of a signal over the range of frequencies. With 100 Hz as the sampling rate, the `signal.welch` function returns 129 features for each sample totaling 1290 features. This is far more manageable than the 60001 features per sensor. This characterizes the event into the plot shown below for each `segment_id`."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,6))\nplt.title('PSD')\nplt.xlabel('Frequency')\nplt.ylabel('Power')\nplt.tight_layout()\n\nexample_PSD = []\n\nfor col in data_columns:\n    col_mean = df_example[col].mean()\n    df_example[col].fillna(col_mean, inplace=True)\n    freq, psd = signal.welch(df_example[col],500)\n    plt.loglog(freq,psd)\n    example_PSD.append(psd)\n    \nplt.legend(data_columns)\n\nprint('Length of PSD: {}'.format(len(psd)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example_PSD = np.transpose(example_PSD)\nprint(\"Transposed PSD array shape: {}\".format(example_PSD.shape))\ndf_PSD = pd.DataFrame(data=example_PSD, columns=data_columns)\ndf_PSD.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate Training DataFrame\n\nWe'll conduct the power spectral density function on every segment id and flatten it out. We should end up with an array of 4431 samples of 1290 features."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_input = []\ni=0\n\nfor segment in train_labels['segment_id']:\n    output_psd = []\n    dataframe = pd.read_csv(f'../input/predict-volcanic-eruptions-ingv-oe/train/{segment}.csv')\n    for col in data_columns:\n        freq, psd = signal.welch(dataframe[col],100)\n        output_psd = np.append(output_psd,psd)\n    \n    train_input = np.append(train_input,output_psd)\n    \n    i=i+1\n    #print('Manipulating segment {}, {} out of {}'.format(segment,i,len(train_labels['segment_id'])))\n\nnum_features = len(output_psd)\ntrain_input = np.reshape(train_input,(i,num_features))\nprint('Finalized input shape: {}'.format(train_input.shape))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Shape of training input: {}\".format(train_input.shape))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate Testing DataFrame\n\nNow we'll do the exact same for all the test data. The files will be read in sequence of the `sample_submission.csv`."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_sample = pd.read_csv(\"../input/predict-volcanic-eruptions-ingv-oe/sample_submission.csv\")\ntest_sample.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_input = []\nk=0\n\nfor segment in test_sample['segment_id']:\n    output_psd = []\n    dataframe = pd.read_csv(f'../input/predict-volcanic-eruptions-ingv-oe/test/{segment}.csv')\n    for col in data_columns:\n        freq, psd = signal.welch(dataframe[col],100)\n        output_psd = np.append(output_psd,psd)\n    \n    test_input = np.append(test_input,output_psd)\n    \n    k=1+k\n    #print('Manipulating segment {}, {} out of {}'.format(segment,k,len(test_sample['segment_id'])))\n\n\nnum_features = len(output_psd)\ntest_input = np.reshape(test_input,(-1,num_features))\nprint('Finalized input shape: {}'.format(test_input.shape))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Shape of test input: {}\".format(test_input.shape))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training and Test Uniformity and Distribution\n\nThere seems to be a discrepency between the distribution of data between the two sets. We'll do a ks_2samp test to illustrate that. We'll take a look at this for each of the 129 features for the 10 sensors and plot it against the 10 sensors."},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import stats\n\nks_value = []\np_value = []\nindex = range(0,1290)\n\nfor i in index:\n    train = train_input[:,i]\n    test = test_input[:,i]\n    statistic,pvalue = stats.ks_2samp(train,test)\n    ks_value = np.append(ks_value,statistic)\n    p_value = np.append(p_value,pvalue)\n    \nsensor_ks = np.reshape(ks_value,(-1,129))\nsensor_p = np.reshape(p_value,(-1,129))\n\nsensor_array = []\nfor i in range(len(data_columns)):\n    sensor_array = np.append(sensor_array,np.full((129),i+1))\n    \nsensor_df = pd.DataFrame({'sensor_id':sensor_array,'ks_value':ks_value,'p_value':p_value})\n\nfig, axes = plt.subplots(1, 2,figsize=(15, 5))\nfig.suptitle(\"K Statistic and P-values for Train and Test\")\nsns.stripplot(ax=axes[0],x='sensor_id',y='ks_value',data=sensor_df)\nsns.stripplot(ax=axes[1],x='sensor_id',y='p_value',data=sensor_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are a few sensors that have decent variation and have non-uniform distribution mong the train and test data. We can also plot the features and see we have not addressed the `NaN` within the tables."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(nrows=5, ncols=2)\nfig.set_size_inches(20,10)\nfig.subplots_adjust(hspace=0.5)\n\nindex_3 = test_input[3,:]\nindex_7 = test_input[7,:]\n\nindex_3 = np.reshape(index_3,(10,129))\nindex_3 = np.transpose(index_3)\nindex_3 = pd.DataFrame(data=index_3,columns=data_columns)\n\nindex_7 = np.reshape(index_7,(10,129))\nindex_7 = np.transpose(index_7)\nindex_7 = pd.DataFrame(data=index_7,columns=data_columns)\n\nfor col,ax in zip(data_columns, axs.flatten()):\n    ax.loglog(freq,index_3[col])\n    ax.loglog(freq,index_7[col])\n    ax.set_title(col)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Standard Scaler\n\nThe best performance I found was simply replacing the `NaN` with zeroes. Below I conduced a Nonlinear Kernal PCA to reduce the features and better compair test/train segments. I had tried a KNN Imputer but it seems there is too much variation between the test/train segments to help in the modell itself. "},{"metadata":{"trusted":true},"cell_type":"code","source":"train_input = np.nan_to_num(train_input)\ntest_input = np.nan_to_num(test_input)\ny_train = train_labels['time_to_eruption']\n\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nscaler.fit(train_input)\n\nX_train = scaler.transform(train_input)\nX_test = scaler.transform(test_input)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Nonlinear Kernal PCA\n\nWith the high number of features, I used a Nonlinear PCA to reduce the number of featurs in order to better visualize the relationship between train and test segments."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.decomposition import KernelPCA\ntransformer = KernelPCA(n_components=100,kernel=\"linear\")\ntrain_transformed = transformer.fit_transform(X_train)\ntest_transformed = transformer.transform(X_test)\n\ncombo_input = np.vstack((train_transformed,test_transformed))\n\ntrain_seg_list = train_labels['segment_id'].to_list()\ntest_seg_list = test_sample['segment_id'].to_list()\ncombo_seg = np.array(train_seg_list)\ncombo_seg = np.append(combo_seg,test_seg_list)\ncombo_seg","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.covariance import EmpiricalCovariance, MinCovDet\n\nrobust_cov = MinCovDet().fit(combo_input[:,:])\n\nm = robust_cov.mahalanobis(combo_input[:,:])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,12))\nplt.title('Train versus Test Scatter of Principal Components')\nplt.xlabel('PCA Component 1')\nplt.ylabel('PCA Component 2')\n\ncm = plt.cm.get_cmap('viridis')\nplt.scatter(train_transformed[:,0],train_transformed[:,1], c=m[:train_input.shape[0]], cmap=cm,s=100)\nplt.scatter(test_transformed[:,0],test_transformed[:,1], c=m[train_input.shape[0]:], cmap=cm,marker=\"P\",s=100)\nplt.colorbar()\n\ntest_x = test_transformed[:,0]\ntest_y = test_transformed[:,1]\n\nfor i,x,y in zip(range(0,len(test_sample['segment_id'])),test_x,test_y):\n    if ((x > 50) | (y > 50)):\n        label = test_sample['segment_id'][i]\n        #plt.annotate(label,(x,y),ha=\"left\")\n        \ntrain_x = train_transformed[:,0]\ntrain_y = train_transformed[:,1]\n        \nfor i,x,y in zip(range(0,len(train_labels['segment_id'])),train_x,train_y):\n    if ((x > 50) | (y > 50)):\n        label = train_labels['segment_id'][i]\n        #plt.annotate(label,(x,y),ha=\"left\")\n        #outlier_list.append(label)\nplt.legend(['train','test'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"outlier_list = []\nplt.figure(figsize=(20,12))\nplt.title('Train versus Test Scatter of Principal Components')\nplt.xlabel('PCA Component 1')\nplt.ylabel('PCA Component 2')\n\ncm = plt.cm.get_cmap('viridis')\nplt.scatter(train_transformed[:,0],train_transformed[:,1],alpha=.6)\nplt.scatter(test_transformed[:,0],test_transformed[:,1], alpha=.4)\n\n\ntest_x = test_transformed[:,0]\ntest_y = test_transformed[:,1]\n\nfor i,x,y in zip(range(0,len(test_sample['segment_id'])),test_x,test_y):\n    if ((x > 100) | (y>100)):\n        label = test_sample['segment_id'][i]\n        plt.annotate(label,(x,y),ha=\"left\")\n        \ntrain_x = train_transformed[:,0]\ntrain_y = train_transformed[:,1]\n        \nfor i,x,y in zip(range(0,len(train_labels['segment_id'])),train_x,train_y):\n    if ((x > 100) | (y>100)):\n        label = train_labels['segment_id'][i]\n        plt.annotate(label,(x,y),ha=\"left\")\n        outlier_list.append(label)\nplt.legend(['train','test'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# List of unclustered segments.\n\noutlier_list","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Random Forest Regression"},{"metadata":{},"cell_type":"markdown","source":"We see there is a wide scattering between number of test and training segments. Most pronounced is there is a high concentration of training semgents along the y axis which have no similar training sets. This may be due to the high number of 'NaN' in the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"#from sklearn.model_selection import cross_val_score\n#from sklearn.model_selection import RepeatedKFold\n#from sklearn.ensemble import RandomForestRegressor\n\n#model = RandomForestRegressor(max_features=700,criterion='mae',random_state=42,\n#                              max_samples=0.8,n_jobs=-1,min_samples_leaf=3)\n\n#cv = RepeatedKFold(n_splits=5, n_repeats=3, random_state=1)\n\n#n_scores = cross_val_score(model,X_train,y_train,scoring='neg_mean_absolute_error',\n#                          cv=cv, n_jobs=-1, error_score='raise',\n#                          verbose=10)\n\n#print('MAE: %.3f (%.3f)' % (np.mean(n_scores), np.std(n_scores)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import joblib\nfrom joblib import dump,load\nmodel = joblib.load('../input/reg-model/original_reg_model.joblib')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.fit(X_train,y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = model.predict(X_test).astype('int64')\ndf_submit = test_sample.copy()\ndf_submit['time_to_eruption'] = abs(predictions)\ndf_submit.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Minimum event time is: {}\".format(df_submit['time_to_eruption'].min()))\nprint(\"Maximum event time is: {}\".format(df_submit['time_to_eruption'].max()))\ndf_submit.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The random forest appears to handle the missing sensors better. I've had scores with CV and LB differencing by 5M but with this model I've at least been able to reduce the gap between the CV & LB by less than 3 million. One option may be to explore other feature characteristics that align the train and test set better."}],"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}