{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Earth Quake Prediction\n\n## Overview\nEvery year many lives are lost and infastructure worh billions of dollars is destroyed in earh quakes. What if we could predict when an earh quake will occur so that this demage can be avoided or reduced.\n\nIn this project we predict the time remaining before laboratory earthquakes occurs from real-time seismic data.\n\nFor details about the kaggle challene please visit the below link;\n\nhttps://www.kaggle.com/c/LANL-Earthquake-Prediction"},{"metadata":{},"cell_type":"markdown","source":"## Data Description\n\nWe will use the (acoustic_data) input signal to predict the time remaining before the next laboratory earthquake (time_to_failure). Please note that all this data has been generated in lab. If you want to know how a laboratory earth quake happens then please check out the below youtube video.\n\nhttps://www.youtube.com/watch?v=m_dBwwDJ4uo\n\nThe training data is a single, continuous segment of experimental data. The test data consists of a folder containing many small segments. The data within each test file is continuous, but the test files do not represent a continuous segment of the experiment; thus, the predictions cannot be assumed to follow the same regular pattern seen in the training file.\n\nFor each seg_id in the test folder, we have to predict a single time_to_failure corresponding to the time between the last row of the segment and the next laboratory earthquake.\n\n**Data Fields**\n\n* acoustic_data - the seismic signal [int16]\n* time_to_failure - the time (in seconds) until the next laboratory earthquake [float64]\n* seg_id - the test segment ids for which predictions should be made (one prediction per segment)"},{"metadata":{},"cell_type":"markdown","source":"### Package Imports and EDA"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#making the imports / ignore warnings\nimport pandas as pd\nimport numpy as np\nimport gc\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom IPython.display import Image\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# although not related I will show some images to get an idea about where the most earthquakes occur\nglobal_earth_quakes = Image('../input/earth-quake-images/global_earth_quakes.jpg', width = 1000)\nglobal_earth_quakes","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see in this heat map that most of the earth quakes happen in Asia Pacific and South America Region."},{"metadata":{"trusted":true},"cell_type":"code","source":"#image showing nuclear plant locations and earth quake hot zones\nnuclear_plants_locations = Image('../input/earth-quake-images/earth_quakes_nuclear_p_locations.jpg')\nnuclear_plants_locations","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From this image we can tell that nuclear plants located in Japan are at risk (if they get some demage in earth quake). Just as a precaution nuclear facilites should not be developed in high risk areas."},{"metadata":{},"cell_type":"markdown","source":"### Reading the training file"},{"metadata":{"trusted":true},"cell_type":"code","source":"#reading the training file (warning: huge size) specify data types to save memory\n#I will be using garbage collection frequently to clear the memory\n\ndata_type = {'acoustic_data': np.int16, 'time_to_failure': np.float32}\ntrain = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', dtype=data_type)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot to see the relation between given variable and target variable\n\ntrain_ad_sample_df = train['acoustic_data'].values[::1000]\ntrain_ttf_sample_df = train['time_to_failure'].values[::1000]\n\n#function for plotting based on both features\ndef plot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% sampled data\"):\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(title)\n    plt.plot(train_ad_sample_df, color='r')\n    ax1.set_ylabel('acoustic data', color='r')\n    plt.legend(['acoustic data'], loc=(0.01, 0.95))\n    ax2 = ax1.twinx()\n    plt.plot(train_ttf_sample_df, color='g')\n    ax2.set_ylabel('time to failure', color='g')\n    plt.legend(['time to failure'], loc=(0.01, 0.9))\n    plt.grid(True)\n\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#delete the old frame\ndel train_ad_sample_df\ndel train_ttf_sample_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot to show zoomed in view\n\ntrain_ad_sample_df = train['acoustic_data'].values[:6291455]\ntrain_ttf_sample_df = train['time_to_failure'].values[:6291455]\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% of data\")\ndel train_ad_sample_df\ndel train_ttf_sample_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that there are series of little jumps in the acoustic data and then there is a big spike which is followed by failure event. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature Enineering"},{"metadata":{"trusted":true},"cell_type":"code","source":"#lets create a function to generate some statistical features based on the training data\n# this is necessay as only one variable [acoustic data] is given to us in training set\n\ndef generate_features(X):\n    strain = []\n    strain.append(X.mean())\n    strain.append(X.std())\n    strain.append(X.min())\n    strain.append(X.max())\n    strain.append(X.kurtosis())\n    strain.append(X.skew())\n    strain.append(np.quantile(X,0.01))\n    strain.append(np.quantile(X,0.05))\n    strain.append(np.quantile(X,0.95))\n    strain.append(np.quantile(X,0.99))\n    strain.append(np.abs(X).max())\n    strain.append(np.abs(X).mean())\n    strain.append(np.abs(X).std())\n    return pd.Series(strain)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# check the head\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets apply feature generation function\n# also we will read the training file in chunks. chunk size specifies the number of rows which pandas will \n# read in one chunk\n\nc_s = 10 ** 6\ntrain = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', iterator=True, chunksize= c_s, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n\nX_train = pd.DataFrame()\ny_train = pd.Series()\nfor df in train:\n    ch = generate_features(df['acoustic_data'])\n    X_train = X_train.append(ch, ignore_index=True)\n    y_train = y_train.append(pd.Series(df['time_to_failure'].values[-1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#describe the data\n\nX_train.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Catboost "},{"metadata":{},"cell_type":"markdown","source":"CatBoost is an algorithm for gradient boosting on decision trees. It is developed by Yandex researchers and engineers (Russian), and is used at Yandex for search, recommendation systems, personal assistant, self-driving cars, weather prediction and many other tasks. It is in open-source and can be used by anyone now.\n\nMore details @\n\nhttps://catboost.ai/docs/concepts/about.html\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# just a base line for cat boost\n# get the best score without hyper parameter tuning\n\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV\n\n\ntrain_pool = Pool(X_train, y_train)\nm = CatBoostRegressor(iterations=10000, loss_function='MAE', boosting_type='Ordered')\nm.fit(X_train, y_train, silent=True)\nm.best_score_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Support Vector Machines"},{"metadata":{},"cell_type":"markdown","source":"The objective of the support vector machine algorithm is to find a hyperplane in an N-dimensional space(N — the number of features) that distinctly classifies the data points. It seems line it is for classification tasks but \nit can be used for regression also. \n\nMore details @\n\nhttps://towardsdatascience.com/support-vector-machine-introduction-to-machine-learning-algorithms-934a444fca47\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# now lets try SVM with rbf kernel + grid search for hyper paramter tuning\n\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.model_selection import KFold\n\nscaler = StandardScaler()\nscaler.fit(X_train)\nX_train_scaled = scaler.transform(X_train)\n\nfolds = KFold(n_splits= 5, shuffle= True, random_state= 101)\n\nparameters = [{'gamma': [0.001, 0.005, 0.01, 0.02, 0.05, 0.1],\n               'C': [0.1, 0.2, 0.25, 0.5, 1, 1.5, 2]}]\n               \n\nreg1 = GridSearchCV(SVR(kernel='rbf', tol=0.01), parameters, cv= folds, scoring='neg_mean_absolute_error')\nreg1.fit(X_train_scaled, y_train.values.flatten())\ny_pred1 = reg1.predict(X_train_scaled)\n\nprint(\"Best CV score: {:.4f}\".format(reg1.best_score_))\nprint(reg1.best_params_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Neural Nets"},{"metadata":{},"cell_type":"markdown","source":"Go through the below link to get an idead about the NNs. \n\nhttps://towardsdatascience.com/machine-learning-for-beginners-an-introduction-to-neural-networks-d49f22d238f9"},{"metadata":{"trusted":true},"cell_type":"code","source":"#making the imports\n#TQDM is a progress bar library with good support for nested loops and Jupyter/IPython notebooks.\n\nfrom sklearn.preprocessing import StandardScaler\nfrom keras.models import Sequential\nfrom keras.layers import Dense\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#reading the training file with data types int16 and float32\n\ndata_type = {'acoustic_data': np.int16, 'time_to_failure': np.float32}\ntrain_data = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', dtype=data_type)\ntrain_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# making the training file ready to be fed into a NN\n\nrows = 150000\nsegments = int(np.floor(train_data.shape[0] / rows))\n\nX_train = pd.DataFrame(index = range(segments),dtype = np.float32,columns = ['mean','std','99quat','50quat','25quat','1quat'])\ny_train = pd.DataFrame(index = range(segments),dtype = np.float32,columns = ['time_to_failure'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# generating the features like mean/std/quantiles\n\nfor segment in tqdm(range(segments)):\n    x = train_data.iloc[segment*rows:segment*rows+rows]\n    y = x['time_to_failure'].values[-1]\n    x = x['acoustic_data'].values\n    X_train.loc[segment,'mean'] = np.mean(x)\n    X_train.loc[segment,'std']  = np.std(x)\n    X_train.loc[segment,'99quat'] = np.quantile(x,0.99)\n    X_train.loc[segment,'50quat'] = np.quantile(x,0.5)\n    X_train.loc[segment,'25quat'] = np.quantile(x,0.25)\n    X_train.loc[segment,'1quat'] =  np.quantile(x,0.01)\n    y_train.loc[segment,'time_to_failure'] = y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#using standard scaler to scale the data\n\nscaler = StandardScaler()\nX_scaler = scaler.fit_transform(X_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collection\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#compiling the sequential model. Simple model with input shape 6 and activation function rectified linear\n# as it is a regression task so use Mean Absolute Error as measuring matrix\n# will use default optimizer adam\n\nmodel = Sequential()\nmodel.add(Dense(32,input_shape = (6,),activation = 'relu'))\nmodel.add(Dense(32,activation = 'relu'))\nmodel.add(Dense(32,activation = 'relu'))\nmodel.add(Dense(1))\nmodel.compile(loss = 'mae',optimizer = 'adam')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#train the model (30 epochs)\n#feed in the scaled training data\n\nmodel.fit(X_scaler,y_train.values.flatten(),epochs = 30)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#reading the submission file from input directory\n\nsub_data = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv',index_col = 'seg_id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#building the test data frame using same columns as X_train\n\nX_test = pd.DataFrame(columns = X_train.columns,dtype = np.float32,index = sub_data.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#feature generation for test data\n\nfor seq in tqdm(X_test.index):\n    test_data = pd.read_csv('../input/LANL-Earthquake-Prediction/test/'+seq+'.csv')\n    x = test_data['acoustic_data'].values\n    X_test.loc[seq,'mean'] = np.mean(x)\n    X_test.loc[seq,'std']  = np.std(x)\n    X_test.loc[seq,'99quat'] = np.quantile(x,0.99)\n    X_test.loc[seq,'50quat'] = np.quantile(x,0.5)\n    X_test.loc[seq,'25quat'] = np.quantile(x,0.25)\n    X_test.loc[seq,'1quat'] =  np.quantile(x,0.01)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#garbage collect\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#scale the test data using pre-defined scaler\n\nX_test_scaler = scaler.transform(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#making the predictions\npred = model.predict(X_test_scaler)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Will not submit this file. Lets wait for Xgboost with hyper parameter tuning to give us the best predictions. "},{"metadata":{},"cell_type":"markdown","source":"## Xgboost"},{"metadata":{},"cell_type":"markdown","source":"XGBoost is well known to provide better solutions than other machine learning algorithms. In fact, since its inception, it has become the \"state-of-the-art” machine learning algorithm to deal with structured data.\n\nFor hyper parameter tuning of xgboost go through the below link:\n\nhttps://towardsdatascience.com/fine-tuning-xgboost-in-python-like-a-boss-b4543ed8b1e"},{"metadata":{"trusted":true},"cell_type":"code","source":"#import xgboost (we will use regressor)\n\nimport xgboost as xgb","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#use the same already scaled X and y from NN part\n\nxgb_model = xgb.XGBRegressor()\n\nxgb_model.fit(X_scaler,y_train.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#predictions without hyper paramter tuning\n\npred = xgb_model.predict(X_test_scaler)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# hyperparameter tuning with XGBoost (will take some time to run)\n\n# creating a KFold object with 3 splits\n\nfolds = KFold(n_splits= 3, shuffle= True, random_state= 101)\n\n# specify range of hyperparameters\nparam_grid = {'learning_rate': [0.01, 0.1, 0.2, 0.3], \n             'subsample': [0.3, 0.6, 0.9, 1],\n              'n_estimators' : [5, 10, 15, 20],\n              'max_depth' :[2,4,6,8]          \n             }          \n\n\n# specify model\nxgb_model = xgb.XGBRegressor()\n\n# set up GridSearchCV()\nmodel_cv = GridSearchCV(estimator = xgb_model, \n                        param_grid = param_grid, \n                        scoring='neg_mean_absolute_error', \n                        cv = folds, \n                        verbose = 1,\n                        return_train_score=True, \n                        n_jobs= -1)      ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#train the model\nmodel_cv.fit(X_scaler,y_train.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# printing the optimal accuracy score and hyperparameters\nprint('We can get neg_mean_absolute_error:',model_cv.best_score_,'using',model_cv.best_params_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# define model with best paramters and train plus make predictions\n\nxgb_model = xgb.XGBRegressor(learning_rate= 0.2, max_depth= 4, n_estimators= 10, subsample= 0.9)\n\nxgb_model.fit(X_scaler,y_train.values)\n\npred = xgb_model.predict(X_test_scaler)\n\n# read the submission file and populate it with predictions\n\nsample_submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsample_submission['time_to_failure'] = pred\n\nprint(sample_submission.shape)\n\nprint('\\n')\n\nsample_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#write to csv file\n\nsample_submission.to_csv('Final_EQ_sub.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}