{"cells":[{"metadata":{},"cell_type":"markdown","source":"**Basic Information**\n\nIn this notebook, I build a Neural Network architecture to predict TTF.\n\nThe goal of this competition is to use seismic signals to predict the timing of laboratory earthquakes. The data comes from a well-known experimental set-up used to study earthquake physics. The acoustic_data input signal is used to predict the time remaining before the next laboratory earthquake (time_to_failure).\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, you should predict a single time_to_failure corresponding to the time between the last row of the segment and the next laboratory earthquake.\n\n* train length: 629,145,480\n* Max time_to_failure = 16.1074\n* Min time_to_failure = 9.5503965e-05\n* test length: 2624 * 150,000 = 393,600,000\n<br>\n"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport feather\nimport os\nimport gc\nimport multiprocessing\nfrom tqdm import tqdm\nfrom numba import jit\nfrom keras import Model, Sequential\nfrom keras.layers import Dense, Flatten, BatchNormalization, Dropout, Activation\nfrom keras.layers import Conv1D, SeparableConv1D, MaxPooling1D, GlobalAveragePooling1D\nfrom keras.layers import Input,Concatenate,Reshape,CuDNNLSTM,CuDNNGRU,GlobalMaxPooling1D\nfrom keras.layers import PReLU, LeakyReLU\nfrom keras.optimizers import adam, rmsprop\nfrom keras.regularizers import l1,l2, l1_l2\nfrom keras.callbacks import ModelCheckpoint\nfrom keras.models import load_model\nfrom scipy.stats import *\nfrom sklearn.metrics import mean_absolute_error\n\n#from scipy import rfft\nfrom numpy.fft import *\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom numpy.random import seed\nseed(1337)\nfrom tensorflow import set_random_seed\nset_random_seed(1337)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir('../input/LANL-Earthquake-Prediction/'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain = feather.read_dataframe('../input/lanl-ft/train.ft')\n# zero center \ntrain['acoustic_data'] = train['acoustic_data'].values - 4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# to plot training/validation history object\ndef plt_dynamic(x, vy, ty, ax, colors=['b'], title=''):\n    ax.plot(x, vy, 'b', label='Validation Loss')\n    ax.plot(x, ty, 'r', label='Train Loss')\n    plt.legend()\n    plt.grid()\n    plt.title(title)\n    fig.canvas.draw()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the earthquake indices. this is where the experiement resets\ndiff = np.diff(train['time_to_failure'].values)\nend = np.nonzero(diff>0)[0]\nstart = end + 1\nstart = np.insert(start, 0, 0)\ndel diff\ngc.collect()\nstart","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def gen_batches_front(col, interval=150_000):\n    high = []\n    low = []\n    splits =[]\n    high_ttf = list(range(9))\n    \n    for i, beg in enumerate(start):\n        counter = 0\n        if beg != 621_985_673:\n            last = start[i+1]\n        else:\n            last = len(train)\n        last = (last-beg)//150_000 * 150_000 + beg\n        \n        for x in range(beg, last, interval):\n            if col == 'acoustic_data':\n                if i in high_ttf:\n                    high.append(train[col].iloc[x:150_000+x].values)\n                else:\n                    low.append(train[col].iloc[x:150_000+x].values)\n            else:\n                if i in high_ttf:\n                    high.append(train[col].iloc[x:150_000+x].values[-1])\n                else:\n                    low.append(train[col].iloc[x:150_000+x].values[-1])\n                    \n        # oversample the end points\n        sample = 15000\n        seg = 150000\n        for z, y in enumerate(range(0, sample*11, sample)):\n            if col == 'acoustic_data':\n                if i in high_ttf:\n                    high.append(train[col].iloc[last-seg*(z+1):last-seg*z].values)\n                else:\n                    low.append(train[col].iloc[last-seg*(z+1):last-seg*z].values)\n            else:\n                if i in high_ttf:\n                    high.append(train[col].iloc[last-seg*(z+1):last-seg*z].values[-1])\n                else:\n                    low.append(train[col].iloc[last-seg*(z+1):last-seg*z].values[-1])\n    return np.asarray(high), np.asarray(low)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# can modify the interval to a factor of 150000 for increased sampling\ndef preprocess_front():\n    xtrain, xtest = gen_batches_front('acoustic_data', interval=150000)\n    ytrain, ytest = gen_batches_front('time_to_failure', interval=150000)\n    xtrain = xtrain.reshape(-1, 150000, 1)\n    xtest = xtest.reshape(-1, 150000, 1)\n    print(xtrain.shape)\n    print(xtest.shape)\n    print(ytrain.shape)\n    print(ytest.shape)\n    return xtrain, xtest, ytrain, ytest","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xtrain, xtest, ytrain, ytest = preprocess_front()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()\ndel train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"checkpoint1 = ModelCheckpoint('best1.hdf5', verbose=0, save_best_only=True, mode='min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nepochs=20\nmdl = Sequential()\nmdl.add(SeparableConv1D(32, 8, activation='relu', input_shape=(xtrain.shape[1],1)))\nmdl.add(CuDNNGRU(32, return_sequences=True))\nmdl.add(GlobalAveragePooling1D())\nmdl.add(Dense(32, activation='relu'))\nmdl.add(Dense(1))\nmdl.compile(loss='mae', optimizer=adam(lr=0.001))\nmdl.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#visualize network architecture\nfrom IPython.display import SVG\nfrom keras.utils.vis_utils import model_to_dot\nSVG(model_to_dot(mdl).create(prog='dot', format='svg'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = mdl.fit(xtrain, ytrain, epochs=epochs, batch_size=64, verbose=2, \n                  validation_data=[xtest,ytest], callbacks=[checkpoint1])\nfig, ax = plt.subplots(1,1)\nvy = history.history['val_loss']\nty = history.history['loss']\nax.set_xlabel('Epoch')\nx = list(range(1,epochs+1))\nax.set_ylabel('Mean Absolute Error')\nplt_dynamic(x,vy,ty,ax, title='history')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x_train = np.concatenate([xtrain, xtest], axis=0)\ny_train = np.concatenate([ytrain, ytest], axis=0)\noof_pred = np.zeros(len(x_train))\nbest_oof = np.zeros(len(x_train))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain_pred = mdl.predict(xtrain)\nvalid_pred = mdl.predict(xtest)\noof_pred[xtrain.shape[0]:] += valid_pred[:,0]\nprint(f'Train MAE: {mean_absolute_error(ytrain, train_pred):.4f}')\nprint(f'Valid MAE: {mean_absolute_error(ytest, valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted train')\nsns.distplot(train_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted test')\nsns.distplot(valid_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#compute absolute error\ntrain_error = np.abs(np.subtract(train_pred.reshape(1,-1)[0], ytrain))\nvalid_error = np.abs(np.subtract(valid_pred.reshape(1,-1)[0], ytest))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the train error distribution\nplt.title('train error')\nplt.xlabel('absolute error')\nsns.distplot(train_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the validation error distribution\nplt.title('validation error')\nplt.xlabel('absolute error')\nsns.distplot(valid_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'max ytrain: {np.max(ytrain):.4f}')\nprint(f'max ytest: {np.max(ytest):.4f}')\nprint(f'max p_train: {np.max(train_pred):.4f}')\nprint(f'max p_test: {np.max(valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot scatter\nplt.figure(figsize=(9,9))\nplt.scatter(ytrain, train_pred)\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot scatter\nplt.figure(figsize=(9,9))\nplt.scatter(ytest, valid_pred)\nplt.title('oof predictions')\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nbest1 = load_model('best1.hdf5')\ntrain_pred = best1.predict(xtrain)\nvalid_pred = best1.predict(xtest)\nbest_oof[xtrain.shape[0]:] += valid_pred[:,0]\nprint(f'Train MAE: {mean_absolute_error(ytrain, train_pred):.4f}')\nprint(f'Valid MAE: {mean_absolute_error(ytest, valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted train')\nsns.distplot(train_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted test')\nsns.distplot(valid_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#compute absolute error\ntrain_error = np.abs(np.subtract(train_pred.reshape(1,-1)[0], ytrain))\nvalid_error = np.abs(np.subtract(valid_pred.reshape(1,-1)[0], ytest))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the train error distribution\nplt.title('train error')\nplt.xlabel('absolute error')\nsns.distplot(train_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the validation error distribution\nplt.title('validation error')\nplt.xlabel('absolute error')\nsns.distplot(valid_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'max ytrain: {np.max(ytrain):.4f}')\nprint(f'max ytest: {np.max(ytest):.4f}')\nprint(f'max p_train: {np.max(train_pred):.4f}')\nprint(f'max p_test: {np.max(valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_pred, valid_pred, train_error, valid_error\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"checkpoint2 = ModelCheckpoint('best2.hdf5', verbose=0, save_best_only=True, mode='min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nepochs=20\nmdl2 = Sequential()\nmdl2.add(SeparableConv1D(32, 8, activation='relu', input_shape=(xtrain.shape[1],1)))\nmdl2.add(CuDNNGRU(32, return_sequences=True))\nmdl2.add(GlobalAveragePooling1D())\nmdl2.add(Dense(32, activation='relu'))\nmdl2.add(Dense(1))\nmdl2.compile(loss='mae', optimizer=adam(lr=0.001))\nhistory = mdl2.fit(xtest, ytest, epochs=epochs, batch_size=64, verbose=2, \n                  validation_data=[xtrain, ytrain], callbacks=[checkpoint2])\nfig, ax = plt.subplots(1,1)\nvy = history.history['val_loss'][4:]\nty = history.history['loss'][4:]\nax.set_xlabel('Epoch')\nx = list(range(5,epochs+1))\nax.set_ylabel('Mean Absolute Error')\nplt_dynamic(x,vy,ty,ax, title='history')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain_pred = mdl2.predict(xtrain)\nvalid_pred = mdl2.predict(xtest)\noof_pred[:xtrain.shape[0]] += train_pred[:,0]\nprint(f'Train MAE: {mean_absolute_error(ytrain, train_pred):.4f}')\nprint(f'Valid MAE: {mean_absolute_error(ytest, valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the train distibution\nplt.title('predicted train')\nsns.distplot(train_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the test distribution\nplt.title('predicted test')\nsns.distplot(valid_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#compute absolute error\ntrain_error = np.abs(np.subtract(train_pred.reshape(1,-1)[0], ytrain))\nvalid_error = np.abs(np.subtract(valid_pred.reshape(1,-1)[0], ytest))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the train error distribution\nplt.title('train error')\nplt.xlabel('absolute error')\nsns.distplot(train_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the test error distribution\nplt.title('validation error')\nplt.xlabel('absolute error')\nsns.distplot(valid_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'max ytrain: {np.max(ytrain):.4f}')\nprint(f'max ytest: {np.max(ytest):.4f}')\nprint(f'max p_train: {np.max(train_pred):.4f}')\nprint(f'max p_test: {np.max(valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nbest2 = load_model('best2.hdf5')\ntrain_pred = best2.predict(xtrain)\nvalid_pred = best2.predict(xtest)\nbest_oof[:xtrain.shape[0]] += train_pred[:,0]\nprint(f'Train MAE: {mean_absolute_error(ytrain, train_pred):.4f}')\nprint(f'Valid MAE: {mean_absolute_error(ytest, valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted train')\nsns.distplot(train_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('predicted test')\nsns.distplot(valid_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot scatter - visualize actual vs prediction\nplt.figure(figsize=(9,9))\nplt.scatter(ytrain, train_pred)\nplt.title('oof predictions')\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot scatter - visualize actual vs prediction\nplt.figure(figsize=(9,9))\nplt.scatter(ytest, valid_pred)\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#compute absolute error\ntrain_error = np.abs(np.subtract(train_pred.reshape(1,-1)[0], ytrain))\nvalid_error = np.abs(np.subtract(valid_pred.reshape(1,-1)[0], ytest))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the train error distribution\nplt.title('train error')\nplt.xlabel('absolute error')\nsns.distplot(train_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot the validation error distribution\nplt.title('validation error')\nplt.xlabel('absolute error')\nsns.distplot(valid_error)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'max ytrain: {np.max(ytrain):.4f}')\nprint(f'max ytest: {np.max(ytest):.4f}')\nprint(f'max p_train: {np.max(train_pred):.4f}')\nprint(f'max p_test: {np.max(valid_pred):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_pred, valid_pred, train_error, valid_error\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# OOF CV Results"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('overtrain mae: {:.4f}'.format(mean_absolute_error(y_train, oof_pred)))\nprint('best mae: {:.4f}'.format(mean_absolute_error(y_train, best_oof)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('oof_pred')\nsns.distplot(oof_pred)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('best_oof')\nsns.distplot(best_oof)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(9,9))\nplt.scatter(y_train, oof_pred)\nplt.title('oof_pred')\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(9,9))\nplt.scatter(y_train, best_oof)\nplt.title('best_oof')\nplt.xlabel('ytrue')\nplt.ylabel('ypred')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Submission File\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsub = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', \n                  dtype={'seg_id': 'category', 'time_to_failure':np.float32})\n\ntest_data = []\nfor fname in sub['seg_id'].values:\n    test_data.append(pd.read_csv('../input/LANL-Earthquake-Prediction/test/'+fname+'.csv', \n                                 dtype={'acoustic_data':np.int16})['acoustic_data'].values)\n# zero center & reshape\ntest_data = np.asarray(test_data) - 4\ntest_data = test_data.reshape(-1, 150_000, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#%time \npred1 = best1.predict(test_data)\nsub['time_to_failure'] = pred1\nsub.to_csv('submission1.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('first half')\nsns.distplot(sub['time_to_failure'].values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time pred2 = best2.predict(test_data)\nsub['time_to_failure'] = pred2\nsub.to_csv('submission2.csv', index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('second half')\nsns.distplot(sub['time_to_failure'].values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# blended first + second half\nfirst = pd.read_csv('submission1.csv')\nsecond = pd.read_csv('submission2.csv')\nblend = first.copy()\nblend['time_to_failure'] = (blend['time_to_failure'] + second['time_to_failure'])/2\nblend.to_csv('frontback.csv', index=False)\nblend.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title('blended')\nsns.distplot(blend['time_to_failure'].values)","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}