{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Began from the following starter kernel, https://www.kaggle.com/johnnyd113/baseline-with-explanations-how-to-get-started\n#    processed the data a little differently and utilized a different benchmark\n# Added some additional metrics from the following kernel, https://www.kaggle.com/harshel7/earthquake-prediction-ensemble-nn\n#    used a deep NN for a benchmark\n# optimizers = https://towardsdatascience.com/preventing-deep-neural-network-from-overfitting-953458db800a\n# grid search method = https://blogs.oracle.com/meena/simple-neural-network-model-using-keras-and-grid-search-hyperparameterstuning\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"traindata = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv',\n                    dtype={'acoustic_data': np.int16,\n                           'time_to_failure': np.float64}) \nprint(traindata.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Smaller sample set of training data \n\n#traindata = pd.read_csv('../input/test-sample/train_sample.csv',\n#                    dtype={'acoustic_data': np.int16,\n#                           'time_to_failure': np.float64}) \n#print(traindata.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Loop to find all earthquake events and print the index of the beginning of the next data segment\n\n#prev = 0\n#index = 0\n#event_index = []\n#for index in range(traindata.shape[0]):\n#    row = traindata.iloc[index]\n#    if (row[1] - prev) >= 0.1: # Check for jump in the values to know where the earthquake event stops\n#        event_index.append(index)\n#        print('')\n#        print(index)\n#    if (index % 100000) == 0:\n#        print('*', end='') # Each * represents 100,000 data points, print this to know the loop is running\n#    prev = row[1] \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Visualize the data\nimport matplotlib.pyplot as plt\n\n#\n# Figure 1 - plot all data\n##########################\nx = pd.Series(traindata['acoustic_data'].values[::50])\ny = pd.Series(traindata['time_to_failure'].values[::50])\n\n# source for plotting, https://matplotlib.org/gallery/api/two_scales.html\nplt.figure(1)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_xlabel('data (#)')\nax1.set_ylabel('input (acoustic data)', color=color)\nax1.plot(x, color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"All Data\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data\n\n#\n# Figure 2 - zoom in to 4th event\n######################################\nx = pd.Series(traindata['acoustic_data'].values[104677356:139772453:20])\ny = pd.Series(traindata['time_to_failure'].values[104677356:139772453:20])\n\nplt.figure(2)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_xlabel('data (#)')\nax1.set_ylabel('input (acoustic data)', color=color)\nax1.plot(x, color=color)\n#ax1.plot(np.abs(x), color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"Zoom to 4th event\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#store and print out all of the earthquake events from the training data\n#Data splits will be created from these event points\nseg_size = 150000\n\nevent_index=[5656574,\n50085878,\n104677356,\n138772453,\n187641820,\n218652630,\n245829585,\n307838917,\n338276287,\n375377848,\n419368880,\n461811623,\n495800225,\n528777115,\n585568144,\n621985673]\n\nindex = 0\n\n# update the event_index with the earthquake events instead of the beginnings\nfor eq in event_index:\n    event_index[index] = eq - 1\n    index = index + 1\n# append the final earthquake event\nevent_index.append(traindata.shape[0])\n# event_index now contains all earthquake events\nprint(event_index)\n\n# process the data into equal segments starting from each earthquake\nindex = 0\nevent_index_segments = []\nfor eq in event_index:\n    if index == 0:\n        event_index_segments.append(int(eq / seg_size))\n    else :\n        event_index_segments.append(int((eq - event_index[index-1]) / seg_size))\n    index = index + 1\n\nprint(event_index_segments)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Feature for training data\ncols = ['mean','median','std','max','skew',\n        'min','sum','var','ptp','mean_change_abs','max_to_min','max_to_min_diff',\n        'abs_avg','abs_std','abs_max','abs_min','abs_med','abs_skew','abs_sum','abs_var','abs_ptp',\n        '10p','25p','50p','75p','90p','abs_1p','abs_5p','abs_30p','abs_60p','abs_95p','abs_99p',\n#        'mean_first_10000','mean_last_10000','mean_first_50000','mean_last_50000',\n#        'std_first_10000','std_last_10000','std_first_50000','std_last_50000'\n       ]\n\n# function for condensing segments into statistical data\ndef condense_segment(segment_, x_condensed, y_label, loc_):        \n    #Converts to numpy array\n    x = pd.Series(segment_['acoustic_data'].values)\n\n    #Grabs the final 'time_to_failure' value\n    y = segment_['time_to_failure'].values[-1]\n    y_label.loc[loc_, 'time_to_failure'] = y\n\n    #For every 150,000 rows, we make these calculations\n    x_condensed.loc[loc_, 'mean'] = np.mean(x)\n    x_condensed.loc[loc_, 'median'] = np.median(x)\n    x_condensed.loc[loc_, 'std'] = np.std(x)\n    x_condensed.loc[loc_, 'max'] = np.max(x)\n    x_condensed.loc[loc_, 'skew'] = x.skew()\n    x_condensed.loc[loc_, 'min'] = np.min(x)\n    x_condensed.loc[loc_, 'sum'] = np.sum(x)\n    x_condensed.loc[loc_, 'var'] = np.var(x)\n    x_condensed.loc[loc_, 'ptp'] = np.ptp(x) #Peak-to-peak is like range\n    x_condensed.loc[loc_, 'mean_change_abs'] = np.mean(np.diff(x))\n    x_condensed.loc[loc_, 'max_to_min'] = np.max(x) / np.abs(np.min(x))\n    x_condensed.loc[loc_, 'max_to_min_diff'] = np.max(x) - np.abs(np.min(x))\n\n    x_condensed.loc[loc_, 'abs_avg'] = np.abs(x).mean()\n    x_condensed.loc[loc_, 'abs_std'] = np.abs(x).std()\n    x_condensed.loc[loc_, 'abs_max'] = np.abs(x).max()\n    x_condensed.loc[loc_, 'abs_min'] = np.abs(x).min()\n    x_condensed.loc[loc_, 'abs_med'] = np.abs(x).median()\n    x_condensed.loc[loc_, 'abs_skew'] = np.abs(x).skew()\n    x_condensed.loc[loc_, 'abs_sum'] = np.abs(x).sum()\n    x_condensed.loc[loc_, 'abs_var'] = np.abs(x).var()\n    x_condensed.loc[loc_, 'abs_ptp'] = np.abs(x).ptp() #Peak-to-peak is like range\n\n    x_condensed.loc[loc_, '10p'] = np.percentile(x,q=10) \n    x_condensed.loc[loc_, '25p'] = np.percentile(x,q=25) #We can also grab percentiles\n    x_condensed.loc[loc_, '50p'] = np.percentile(x,q=50)\n    x_condensed.loc[loc_, '75p'] = np.percentile(x,q=75)\n    x_condensed.loc[loc_, '90p'] = np.percentile(x,q=90)   \n\n    x_condensed.loc[loc_, 'abs_1p'] = np.percentile(x, np.abs(0.01))\n    x_condensed.loc[loc_, 'abs_5p'] = np.percentile(x, np.abs(0.05))\n    x_condensed.loc[loc_, 'abs_30p'] = np.percentile(x, np.abs(0.30))\n    x_condensed.loc[loc_, 'abs_60p'] = np.percentile(x, np.abs(0.60))\n    x_condensed.loc[loc_, 'abs_95p'] = np.percentile(x, np.abs(0.95))\n    x_condensed.loc[loc_, 'abs_99p'] = np.percentile(x, np.abs(0.99))\n\n#    x_condensed.loc[loc_, 'mean_first_10000'] = x[:10000].mean()\n#    x_condensed.loc[loc_, 'mean_last_10000']  =  x[-10000:].mean()\n#    x_condensed.loc[loc_, 'mean_first_50000'] = x[:50000].mean()\n#    x_condensed.loc[loc_, 'mean_last_50000'] = x[-50000:].mean()\n\n#    x_condensed.loc[loc_, 'std_first_10000'] = x[:10000].std()\n#    x_condensed.loc[loc_, 'std_last_10000']  =  x[-10000:].std()\n#    x_condensed.loc[loc_, 'std_first_50000'] = x[:50000].std()\n#    x_condensed.loc[loc_, 'std_last_50000'] = x[-50000:].std()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Dataframe for condensed training data\ntotal_segments = sum(event_index_segments)\nprint(\"Total Segments = \" + str(total_segments))\nX_train = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=cols) # Feature list\ny_train = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=['time_to_failure']) #Our target variable\n\n# Condense the segments into the training dataframe\nindex = 0\ni_total = 0 # i counter restarts at 0 for each event_index_segment iteration, keep track of total segments\np_done = 0\nprint(\"Data condensing...\")\nfor segments in event_index_segments: # Load total segments available for this event\n    event_end = event_index[index]  # Earthquake index in the training data\n    index = index + 1\n\n    for i in range(segments):\n        # Creating segments backwards from the earthquake events\n        # Start of segment = event_end - ((i+1)*seg_size)\n        # End of segment = event_end - (i*seg_size)        \n        sample_segment = traindata.iloc[event_end - ((i+1)*seg_size):event_end - (i*seg_size)]\n        condense_segment(sample_segment, X_train, y_train, i_total)\n        # Print status %\n        i_total = i_total + 1\n        if (i_total % int(total_segments/10)) == 0:\n            p_done = p_done + 1\n            print(str(p_done) + \"0%\\r\", end='') # Each * represents 400 data segments\n\nprint(i_total)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Shift and create addiitional alternate training sets, subtracting 1 segment from each event because of shifting\ntotal_segments = total_segments - len(event_index_segments)\nprint(\"Alt Total Segments = \" + str(total_segments))\nX_train_alt1 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=cols) # Feature list\ny_train_alt1 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=['time_to_failure']) #Our target variable\nX_train_alt2 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=cols) # Feature list\ny_train_alt2 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=['time_to_failure']) #Our target variable\nX_train_alt3 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=cols) # Feature list\ny_train_alt3 = pd.DataFrame(index=range(total_segments), dtype=np.float64, columns=['time_to_failure']) #Our target variable\n\n#\n# Condense alt 1 dataframe (offset 50,0000)\n##########\nindex = 0\ni_total = 0 # i counter restarts at 0 for each event_index_segment iteration, keep track of total segments\np_done = 0\noffset = 37500\nprint(\"Alt1 condensing...  offset=\" + str(offset))\nfor segments in event_index_segments: # Load total segments available for this event\n    event_end = event_index[index] - offset  # Earthquake index in the training data\n    index = index + 1\n    segments = segments - 1 # decrement 1 since we are shifting into the last segment\n\n    for i in range(segments):   \n        sample_segment = traindata.iloc[event_end - ((i+1)*seg_size):event_end - (i*seg_size)]\n        condense_segment(sample_segment, X_train_alt1, y_train_alt1, i_total)\n        # Print status %\n        i_total = i_total + 1\n        if (i_total % int(total_segments/10)) == 0:\n            p_done = p_done + 1\n            print(str(p_done) + \"0%\\r\", end='') # Each * represents 400 data segments\n\nprint(i_total)\n\n#\n# Condense alt 2 dataframe (offset 100,0000)\n##########\nindex = 0\ni_total = 0 # i counter restarts at 0 for each event_index_segment iteration, keep track of total segments\np_done = 0\noffset = 37500 * 2\nprint(\"Alt2 condensing...  offset=\" + str(offset))\nfor segments in event_index_segments: # Load total segments available for this event\n    event_end = event_index[index] - offset  # Earthquake index in the training data\n    index = index + 1\n    segments = segments - 1\n\n    for i in range(segments):    \n        sample_segment = traindata.iloc[event_end - ((i+1)*seg_size):event_end - (i*seg_size)]\n        condense_segment(sample_segment, X_train_alt2, y_train_alt2, i_total)\n        # Print status %\n        i_total = i_total + 1\n        if (i_total % int(total_segments/10)) == 0:\n            p_done = p_done + 1\n            print(str(p_done) + \"0%\\r\", end='') # Each * represents 400 data segments\n\nprint(i_total)\n\n#\n# Condense alt 3 dataframe (offset 100,0000)\n##########\nindex = 0\ni_total = 0 # i counter restarts at 0 for each event_index_segment iteration, keep track of total segments\np_done = 0\noffset = 37500 * 3\nprint(\"Alt3 condensing...  offset=\" + str(offset))\nfor segments in event_index_segments: # Load total segments available for this event\n    event_end = event_index[index] - offset  # Earthquake index in the training data\n    index = index + 1\n    segments = segments - 1\n\n    for i in range(segments):    \n        sample_segment = traindata.iloc[event_end - ((i+1)*seg_size):event_end - (i*seg_size)]\n        condense_segment(sample_segment, X_train_alt3, y_train_alt3, i_total)\n        # Print status %\n        i_total = i_total + 1\n        if (i_total % int(total_segments/10)) == 0:\n            p_done = p_done + 1\n            print(str(p_done) + \"0%\\r\", end='') # Each * represents 400 data segments\n\nprint(i_total)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_alt1.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_alt2.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_alt3.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Visualize the data\nimport matplotlib.pyplot as plt\n\n#\n# Figure 4 - condensed data visualization\n######################################\n\nplt.figure(4)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_xlabel('data (#)')\nax1.set_ylabel('input (acoustic abs mean)', color=color)\nax1.plot(X_train['abs_avg'].values[0:2000], color=color)\nax1.plot(X_train_alt1['abs_avg'].values[0:2000], color='green')\nax1.plot(X_train_alt2['abs_avg'].values[0:2000], color='yellow')\nax1.plot(X_train_alt3['abs_avg'].values[0:2000], color='orange')\n#ax1.plot(np.abs(x), color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y_train.values[0:2000], color=color)\nax2.plot(y_train_alt1.values[0:2000], color='cyan')\nax2.plot(y_train_alt2.values[0:2000], color='magenta')\nax2.plot(y_train_alt3.values[0:2000], color='black')\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"1ST 2000 of Condensed Data - ABS Mean (all 3 data structures)\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Split the training data for additional tests to tune performance\nfrom sklearn.model_selection import train_test_split\n\nX_train_s, X_test, y_train_s, y_test = train_test_split(X_train, y_train, test_size=0.2, random_state=1)\nX_train_s, X_val, y_train_s, y_val = train_test_split(X_train_s, y_train_s, test_size=0.2, random_state=1)\n\n#Alt1 split\nX_train_alt1_s, X_test_alt1, y_train_alt1_s, y_test_alt1 = train_test_split(X_train_alt1, y_train_alt1, test_size=0.2, random_state=1)\nX_train_alt1_s, X_val_alt1, y_train_alt1_s, y_val_alt1 = train_test_split(X_train_alt1_s, y_train_alt1_s, test_size=0.2, random_state=1)\n#Alt2 split\nX_train_alt2_s, X_test_alt2, y_train_alt2_s, y_test_alt2 = train_test_split(X_train_alt2, y_train_alt2, test_size=0.2, random_state=1)\nX_train_alt2_s, X_val_alt2, y_train_alt2_s, y_val_alt2 = train_test_split(X_train_alt2_s, y_train_alt2_s, test_size=0.2, random_state=1)\n#Alt3 split\nX_train_alt3_s, X_test_alt3, y_train_alt3_s, y_test_alt3 = train_test_split(X_train_alt3, y_train_alt3, test_size=0.2, random_state=1)\nX_train_alt3_s, X_val_alt3, y_train_alt3_s, y_val_alt3 = train_test_split(X_train_alt3_s, y_train_alt3_s, test_size=0.2, random_state=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Example Benchmark from the Starter Kernel ###\n#Fit a random forest\n#from sklearn.ensemble import RandomForestRegressor\n\n#This creates the Randomforest with the given parameters\n#rf = RandomForestRegressor(n_estimators=100, #100 trees (Default of 10 is too small)\n#                          max_features=0.5, #Max number of features each tree can use \n#                          min_samples_leaf=30, #Min amount of samples in each leaf\n#                          random_state=11)\n\n#This trains the random forest on our training data\n#rf.fit(X_train,y_train)\n\n### Score the Example Benchmark\n#from sklearn.metrics import mean_absolute_error\n\n#mean_absolute_error(y_val, rf.predict(X_val))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Begin creation of benchmark Simple Deep NN\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Activation\n\nnn_bm_model = Sequential()\nnn_bm_model.add(Dense(128, activation='relu', input_shape=(X_train_s.shape[1],)))\nnn_bm_model.add(Dropout(.2))\nnn_bm_model.add(Dense(64, activation='relu'))\nnn_bm_model.add(Dropout(.1))\nnn_bm_model.add(Dense(1, activation='linear'))\n\n# TODO: Compile the model using a loss function and an optimizer.\nnn_bm_model.compile(loss = 'mean_absolute_error', optimizer='adam', metrics=['mean_absolute_error'])\nnn_bm_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nn_bm_model.fit(X_train_s, y_train_s, epochs=100, batch_size=20, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"score = nn_bm_model.evaluate(X_test, y_test, verbose=0)\nprint(\"        \" + nn_bm_model.metrics_names[0] + \"                \" + nn_bm_model.metrics_names[1])\nprint(\"Test: \", score)\nscore = nn_bm_model.evaluate(X_val, y_val, verbose=0)\nprint(\"Val : \", score)\nscore = nn_bm_model.evaluate(X_test_alt1, y_test_alt1, verbose=0)\nprint(\"Test_alt1 : \", score)\nscore = nn_bm_model.evaluate(X_val_alt1, y_val_alt1, verbose=0)\nprint(\"Val_alt1  : \", score)\nscore = nn_bm_model.evaluate(X_test_alt2, y_test_alt2, verbose=0)\nprint(\"Test_alt2 : \", score)\nscore = nn_bm_model.evaluate(X_val_alt2, y_val_alt2, verbose=0)\nprint(\"Val_alt2  : \", score)\nscore = nn_bm_model.evaluate(X_test_alt3, y_test_alt3, verbose=0)\nprint(\"Test_alt3 : \", score)\nscore = nn_bm_model.evaluate(X_val_alt3, y_val_alt3, verbose=0)\nprint(\"Val_alt3  : \", score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Utilize Grid Search to attempt to tune paramters...\nfrom keras import regularizers\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Activation, BatchNormalization\nfrom keras.wrappers.scikit_learn import KerasClassifier\nfrom sklearn.model_selection import GridSearchCV\n\n# Create Tunable Model\ndef nn_model_(dropout_rate=0.2,optimizer='adam',activation='relu',regularizer=0.0):\n    #default parameters\n    \n    #kernel_regularizer=regularizers.l1(0.)\n    model = Sequential()\n    model.add(Dense(128, activation=activation, input_shape=(X_train_s.shape[1],)))\n    nn_model.add(BatchNormalization()) # added batch normalization to normalize data between layers\n    model.add(Dropout(dropout_rate))\n    model.add(Dense(64, activation=activation))\n    nn_model.add(BatchNormalization())\n    model.add(Dropout(dropout_rate))\n    model.add(Dense(1, activation='linear'))\n    model.compile(loss='mean_absolute_error', optimizer=optimizer, metrics=['mean_absolute_error'])\n    return model\n\nmodel = KerasClassifier(build_fn=nn_model_,epochs=100, batch_size=20)\n# parameters\noptimizer_=['adam','adamax'] ### best - adam\nactivation_=['relu','tanh'] ### best - relu\nepochs_=[80,100,120]\ndropout_rate_=[0.2,0.3,0.4,0.5,0.6] \nbatch_size_=[20,32,48] \nregularizer=[0.0,0.0001,0.001,0.01,0.1]\nparam_grid1 = dict(batch_size=batch_size_,dropout_rate=dropout_rate_)\nparam_grid2 = dict(activation=activation_,optimizer=optimizer_)\nparam_grid3 = dict(epochs=epochs_)\n#grid = GridSearchCV(estimator=model,param_grid=param_grid3,scoring='neg_mean_absolute_error',n_jobs=-1)\n#grid_result = grid.fit(X_train_s, y_train_s)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Print results of grid search\n#print(\"Best: %f using %s\" % (grid_result.best_score_, grid_result.best_params_))\n\n### 128-64 nodes ###\n# param_grid1 results \n# Best: -5.664138 using {'batch_size': 20, 'dropout_rate': 0.2}\n##\n# param_grid2 results (batch_size, dropout_rate set from previous run)\n# Best: -5.664138 using {'activation': 'relu', 'optimizer': 'adam'}\n##\n# param_grid3 results (L1 regularizer)\n# Best: -5.664138 using {'regularizer': 0.0}\n##\n# param_grid4 results (L2 regularizer)\n# Best: -5.664138 using {'regularizer': 0.0}\n\n### 512-256 nodes ###\n# param_grid1 results \n# Best: -5.635036 using {'batch_size': 20, 'dropout_rate': 0.2}\n##\n# param_grid2 results (batch_size, dropout_rate set from previous run)\n# Best: -5.635036 using {'activation': 'relu', 'optimizer': 'adam'}\n##\n# param_grid3 results (epochs)\n# Best: -5.635036 using {'epochs': 80}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Final network with tuned hyperparameters...\nfrom keras import regularizers\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Dropout, Activation, BatchNormalization\nfrom numpy.random import seed\nseed(7)\nfrom tensorflow import set_random_seed\nset_random_seed(77)\n\n#\n# Model Hyperparamters\n#===================\ndropout_rate = 0.20    # 0.2, 0.15, 0.1\nactivation = 'relu'\noptimizer = 'adam'\nepochs = 80      # 80, 100, 250\nbatch_size = 20  # 20, 15\nl1r = 0.1        # 0.1, 0.15, 0.20\n\n# at default, 2 fits seems to do the trick\n# dropout=0.15, 1 fit gave 2.4, 3rd fit 2.5, looks very good\n# dropout=0.1, 3rd fit dropped to 2.4\n\n# 3-layer Deep NN with batch normalization and dropout\nnn_model = Sequential()\nnn_model.add(Dense(64, activation=activation, input_shape=(X_train_s.shape[1],),kernel_regularizer=regularizers.l2(l1r)))\nnn_model.add(BatchNormalization())\nnn_model.add(Dropout(dropout_rate))\nnn_model.add(Dense(32, activation=activation))\nnn_model.add(BatchNormalization())\nnn_model.add(Dropout(dropout_rate))\nnn_model.add(Dense(1, activation='linear'))\n\n# TODO: Compile the model using a loss function and an optimizer.\nnn_model.compile(loss = 'mean_absolute_error', optimizer=optimizer, metrics=['mean_absolute_error'])\nnn_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nn_model.fit(X_train_s, y_train_s, epochs=epochs, batch_size=batch_size, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualize predictions\ny_pred = nn_model.predict(X_train)\n#y_pred_bm = nn_bm_model.predict(X_train)\n\n#\n# Figure 6 - Condensed Training Data Final Prediction\n######################################\n\nplt.figure(6)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_ylim(0,10)\nax1.set_xlabel('data (#)')\nax1.set_ylabel('pred output (s)', color=color)\nax1.plot(y_pred, color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylim(0,10)\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y_train.values, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"Predicted Values for Final Model\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"              \" + nn_model.metrics_names[0] + \"                \" + nn_model.metrics_names[1])\nscore = nn_model.evaluate(X_test, y_test, verbose=0)\nprint(\"Test_s      : \", score)\nscore = nn_model.evaluate(X_val, y_val, verbose=0)\nprint(\"Val_s       : \", score)\nscore = nn_model.evaluate(X_test_alt1, y_test_alt1, verbose=0)\nprint(\"Test_alt1 : \", score)\nscore = nn_model.evaluate(X_val_alt1, y_val_alt1, verbose=0)\nprint(\"Val_alt1  : \", score)\nscore = nn_model.evaluate(X_test_alt2, y_test_alt2, verbose=0)\nprint(\"Test_alt2 : \", score)\nscore = nn_model.evaluate(X_val_alt2, y_val_alt2, verbose=0)\nprint(\"Val_alt2  : \", score)\nscore = nn_model.evaluate(X_test_alt3, y_test_alt3, verbose=0)\nprint(\"Test_alt3 : \", score)\nscore = nn_model.evaluate(X_val_alt3, y_val_alt3, verbose=0)\nprint(\"Val_alt3  : \", score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nn_model.fit(X_train_alt1_s, y_train_alt1_s, epochs=epochs, batch_size=batch_size, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualize predictions\ny_pred = nn_model.predict(X_train)\n#y_pred_bm = nn_bm_model.predict(X_train)\n\n#\n# Figure 6 - Condensed Training Data Final Prediction\n######################################\n\nplt.figure(6)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_ylim(0,10)\nax1.set_xlabel('data (#)')\nax1.set_ylabel('pred output (s)', color=color)\nax1.plot(y_pred, color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylim(0,10)\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y_train.values, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"Predicted Values for Final Model - 2nd fit\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data#","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"              \" + nn_model.metrics_names[0] + \"                \" + nn_model.metrics_names[1])\nscore = nn_model.evaluate(X_test, y_test, verbose=0)\nprint(\"Test_s      : \", score)\nscore = nn_model.evaluate(X_val, y_val, verbose=0)\nprint(\"Val_s       : \", score)\nscore = nn_model.evaluate(X_test_alt1, y_test_alt1, verbose=0)\nprint(\"Test_alt1 : \", score)\nscore = nn_model.evaluate(X_val_alt1, y_val_alt1, verbose=0)\nprint(\"Val_alt1  : \", score)\nscore = nn_model.evaluate(X_test_alt2, y_test_alt2, verbose=0)\nprint(\"Test_alt2 : \", score)\nscore = nn_model.evaluate(X_val_alt2, y_val_alt2, verbose=0)\nprint(\"Val_alt2  : \", score)\nscore = nn_model.evaluate(X_test_alt3, y_test_alt3, verbose=0)\nprint(\"Test_alt3 : \", score)\nscore = nn_model.evaluate(X_val_alt3, y_val_alt3, verbose=0)\nprint(\"Val_alt3  : \", score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nn_model.fit(X_train_alt2_s, y_train_alt2_s, epochs=epochs, batch_size=batch_size, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualize predictions\ny_pred = nn_model.predict(X_train)\n#y_pred_bm = nn_bm_model.predict(X_train)\n\n#\n# Figure 6 - Condensed Training Data Final Prediction\n######################################\n\nplt.figure(6)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_ylim(0,10)\nax1.set_xlabel('data (#)')\nax1.set_ylabel('pred output (s)', color=color)\nax1.plot(y_pred, color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylim(0,10)\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y_train.values, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"Predicted Values for Final Model - 3rd fit\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"              \" + nn_model.metrics_names[0] + \"                \" + nn_model.metrics_names[1])\nscore = nn_model.evaluate(X_test, y_test, verbose=0)\nprint(\"Test_s      : \", score)\nscore = nn_model.evaluate(X_val, y_val, verbose=0)\nprint(\"Val_s       : \", score)\nscore = nn_model.evaluate(X_test_alt1, y_test_alt1, verbose=0)\nprint(\"Test_alt1 : \", score)\nscore = nn_model.evaluate(X_val_alt1, y_val_alt1, verbose=0)\nprint(\"Val_alt1  : \", score)\nscore = nn_model.evaluate(X_test_alt2, y_test_alt2, verbose=0)\nprint(\"Test_alt2 : \", score)\nscore = nn_model.evaluate(X_val_alt2, y_val_alt2, verbose=0)\nprint(\"Val_alt2  : \", score)\nscore = nn_model.evaluate(X_test_alt3, y_test_alt3, verbose=0)\nprint(\"Test_alt3 : \", score)\nscore = nn_model.evaluate(X_val_alt3, y_val_alt3, verbose=0)\nprint(\"Val_alt3  : \", score)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#nn_model.fit(X_train_alt3_s, y_train_alt3_s, epochs=epochs, batch_size=batch_size, verbose=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualize predictions\ny_pred = nn_model.predict(X_train)\n#y_pred_bm = nn_bm_model.predict(X_train)\n\n#\n# Figure 6 - Condensed Training Data Final Prediction\n######################################\n\nplt.figure(6)\nfig, ax1 = plt.subplots(figsize=(14, 7))\n\ncolor = 'tab:red'\nax1.set_ylim(0,10)\nax1.set_xlabel('data (#)')\nax1.set_ylabel('pred output (s)', color=color)\nax1.plot(y_pred, color=color)\nax1.tick_params(axis='y', labelcolor=color)\n\nax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n\ncolor = 'tab:blue'\nax2.set_ylim(0,10)\nax2.set_ylabel('output time (s)', color=color)  # we already handled the x-label with ax1\nax2.plot(y_train.values, color=color)\nax2.tick_params(axis='y', labelcolor=color)\n\nplt.title(\"Predicted Values for Final Model - 4th fit\")\nfig.tight_layout()  # otherwise the right y-label is slightly clipped\nplt.show() #all data#"},{"metadata":{},"cell_type":"markdown","source":"print(\"              \" + nn_model.metrics_names[0] + \"                \" + nn_model.metrics_names[1])\nscore = nn_model.evaluate(X_test, y_test, verbose=0)\nprint(\"Test_s      : \", score)\nscore = nn_model.evaluate(X_val, y_val, verbose=0)\nprint(\"Val_s       : \", score)\nscore = nn_model.evaluate(X_test_alt1, y_test_alt1, verbose=0)\nprint(\"Test_alt1 : \", score)\nscore = nn_model.evaluate(X_val_alt1, y_val_alt1, verbose=0)\nprint(\"Val_alt1  : \", score)\nscore = nn_model.evaluate(X_test_alt2, y_test_alt2, verbose=0)\nprint(\"Test_alt2 : \", score)\nscore = nn_model.evaluate(X_val_alt2, y_val_alt2, verbose=0)\nprint(\"Val_alt2  : \", score)\nscore = nn_model.evaluate(X_test_alt3, y_test_alt3, verbose=0)\nprint(\"Test_alt3 : \", score)\nscore = nn_model.evaluate(X_val_alt3, y_val_alt3, verbose=0)\nprint(\"Val_alt3  : \", score)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Gather and process the testing data to create the submission results...\n\nsubmission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col = 'seg_id')\nX_test_sub = pd.DataFrame(columns = X_train.columns, dtype = np.float64, index = submission.index)\n\nfor id in X_test_sub.index:\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + id + '.csv')\n    x = pd.Series(seg['acoustic_data'].values)\n\n    #For every segment (150,000), make these calculations\n    X_test_sub.loc[id, 'mean'] = np.mean(x)\n    X_test_sub.loc[id, 'median'] = np.median(x)\n    X_test_sub.loc[id, 'std'] = np.std(x)\n    X_test_sub.loc[id, 'max'] = np.max(x)\n    X_test_sub.loc[id, 'skew'] = x.skew()\n    X_test_sub.loc[id, 'min'] = np.min(x)\n    X_test_sub.loc[id, 'sum'] = np.sum(x)\n    X_test_sub.loc[id, 'var'] = np.var(x)\n    X_test_sub.loc[id, 'ptp'] = np.ptp(x) #Peak-to-peak is like range\n    X_test_sub.loc[id, 'mean_change_abs'] = np.mean(np.diff(x))\n    X_test_sub.loc[id, 'max_to_min'] = np.max(x) / np.abs(np.min(x))\n    X_test_sub.loc[id, 'max_to_min_diff'] = np.max(x) - np.abs(np.min(x))\n\n    X_test_sub.loc[id, 'abs_avg'] = np.abs(x).mean()\n    X_test_sub.loc[id, 'abs_std'] = np.abs(x).std()\n    X_test_sub.loc[id, 'abs_max'] = np.abs(x).max()\n    X_test_sub.loc[id, 'abs_min'] = np.abs(x).min()\n    X_test_sub.loc[id, 'abs_med'] = np.abs(x).median()\n    X_test_sub.loc[id, 'abs_skew'] = np.abs(x).skew()\n    X_test_sub.loc[id, 'abs_sum'] = np.abs(x).sum()\n    X_test_sub.loc[id, 'abs_var'] = np.abs(x).var()\n    X_test_sub.loc[id, 'abs_ptp'] = np.abs(x).ptp() #Peak-to-peak is like range\n\n    X_test_sub.loc[id, '10p'] = np.percentile(x,q=10) \n    X_test_sub.loc[id, '25p'] = np.percentile(x,q=25) #We can also grab percentiles\n    X_test_sub.loc[id, '50p'] = np.percentile(x,q=50)\n    X_test_sub.loc[id, '75p'] = np.percentile(x,q=75)\n    X_test_sub.loc[id, '90p'] = np.percentile(x,q=90)   \n\n    X_test_sub.loc[id, 'abs_1p'] = np.percentile(x, np.abs(0.01))\n    X_test_sub.loc[id, 'abs_5p'] = np.percentile(x, np.abs(0.05))\n    X_test_sub.loc[id, 'abs_30p'] = np.percentile(x, np.abs(0.30))\n    X_test_sub.loc[id, 'abs_60p'] = np.percentile(x, np.abs(0.60))\n    X_test_sub.loc[id, 'abs_95p'] = np.percentile(x, np.abs(0.95))\n    X_test_sub.loc[id, 'abs_99p'] = np.percentile(x, np.abs(0.99))\n\n#    X_test_sub.loc[id, 'mean_first_10000'] = x[:10000].mean()\n#    X_test_sub.loc[id, 'mean_last_10000']  =  x[-10000:].mean()\n#    X_test_sub.loc[id, 'mean_first_50000'] = x[:50000].mean()\n#    X_test_sub.loc[id, 'mean_last_50000'] = x[-50000:].mean()\n\n#    X_test_sub.loc[id, 'std_first_10000'] = x[:10000].std()\n#    X_test_sub.loc[id, 'std_last_10000']  =  x[-10000:].std()\n#    X_test_sub.loc[id, 'std_first_50000'] = x[:50000].std()\n#    X_test_sub.loc[id, 'std_last_50000'] = x[-50000:].std()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Predict the test submission data\nnn_predictions = nn_model.predict(X_test_sub)\nnn_bm_prediction = nn_bm_model.predict(X_test_sub)\n\n# Generate Benchmark Prediction\nsubmission['time_to_failure'] = nn_bm_prediction\nprint(submission.head(5))\nsubmission.to_csv('submission_bm.csv')\n\n# Generate Final Prediction\nsubmission['time_to_failure'] = nn_predictions\nprint(submission.head(5))\nsubmission.to_csv('submission_nn.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# At this point, i've submitted both my benchmark and NN for comparison...\n# Frustratingly so, the benchmark did better...  It obtained a score of 2.659\n# I plan to compute MAE between the NN predicted and the benhcmark to hopefully improve the NN\n# Kaggle limits submissions to 2 per day so hoping this helps as a comparison for checking my work\n\n# pull in previously scored benchmark (keep in mind only 20% is scored so this may not work)\ncomp_sub = pd.read_csv('../input/submission-216/submission_nn_216.csv', index_col = 'seg_id')\ntotal_error = 0\nfor id in comp_sub.index:\n    # calculate MAE\n    total_error = total_error + np.abs(comp_sub.loc[id].values[0] - submission.loc[id].values[0])\n# print MAE between NN and benchmark    \nprint(total_error/comp_sub.shape[0])\n\n# initial MAE was 17 prior to adding L1 regularization, the scoring made me think it was overfitting test data\n# L1 with 0.1 gave MAE of 2.01 so \n# L1 with 0.01 gave MAE of 1.055 and the output comparison chart above looked best fit","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}