{"cells":[{"metadata":{"_uuid":"24be44442fc23e7235dc34444b33d3c0a018d46f"},"cell_type":"markdown","source":"## **Can RNNs extract information from seismic signals?**\n\n*[See [here](https://www.kaggle.com/michael422/spectrogram-convolution) for a related article on using [Convolutional Neural Networks](https://www.kaggle.com/michael422/spectrogram-convolution) on seismic data.]*\n\nRecurrent neural networks (RNNs) and variants such as long short-term memory ([LSTM](https://arxiv.org/abs/1402.1128)) are used in a range of applications to interpret sequential data.  Here we will use RNNs to try to predict the time of the next simulated earthquake (\"labquake\") based on the acoustic signals emitted by the laboratory setup.\n\nWe use the first 400MM records in order to process efficiently in-memory:"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom datetime import datetime\nfrom functools import partial\nimport tensorflow as tf\nfrom tqdm import tqdm_notebook, tnrange, tqdm, trange\nfrom scipy import signal\nimport matplotlib.pyplot as plt\nnp.warnings.filterwarnings('ignore')\n\nrng = np.random.RandomState(datetime.now().microsecond)\n\ndata_path = r'../input'\ntf_path = os.path.join(data_path, 'tf')\n\nnrows = 400000000   # 629145480\ndata = (pd.read_csv(os.path.join(data_path, 'train.csv'),\n                    nrows=nrows,\n                    dtype={'acoustic_data': np.int16,\n                           'time_to_failure': np.float32})\n        .values)\nprint(f'input data shape: {data.shape}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f2a7c9892bedde22e40f159fa93feccdd24903b3"},"cell_type":"markdown","source":"Center and log-scale the raw signal data."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"obs = data[:, 0]\nobs_mean = obs.mean()\n\ndef preprocess_obs(obs, obs_mean):\n    obs = obs - obs_mean\n    obs = (np.log(obs * np.sign(obs) + 1)) * np.sign(obs)\n    return obs\n\nprint(f'data recentered about mean={obs_mean:.4f}')\nobs = preprocess_obs(obs, obs_mean)\ndata = np.c_[obs, data[:, 1]]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9f18ef6364c1a345f5032d1e05afaec3bef84ec7"},"cell_type":"markdown","source":"Split train and validation data.  Notice the difference in distribution between the two sets."},{"metadata":{"trusted":true,"_uuid":"afce79799edb2560021bd292915f5f7a874de94f"},"cell_type":"code","source":"TRAIN_VAL_SPLIT = 0.75\n\ntrain_size = int(data.shape[0] * TRAIN_VAL_SPLIT)\ntrain_data = data[:train_size, :]\nval_data = data[train_size:, :]\n\nprint(f'train mean ttf: {train_data[:, 1].mean():.4f}, ' +\n      f'val mean ttf: {val_data[:, -1].mean():.4f}')\nprint(f'train 1st quartile ttf: {np.quantile(train_data[:, -1], 0.25):.4f}, ' +\n      f'val 1st quartile ttf: {np.quantile(val_data[:, -1], 0.25):.4f}')\nprint(f'train 3rd quartile ttf: {np.quantile(train_data[:, -1], 0.75):.4f}, ' +\n      f'val 3rd quartile ttf: {np.quantile(val_data[:, -1], 0.75):.4f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ebb65103b5dcde95a79bb725783b88d87c40286a"},"cell_type":"markdown","source":"The raw signal samples have dimensionality 150,000.  The signals must be reduced to a length that is manageable by an RNN.  We will use a reduced sequence length of 40, which is an integer factor of the test sample length.  Then define a few simple feature extractors that downsample the data by dividing it into frames, returning a sequence of frame-wise aggregates:\n1. the range (max - min) of the log amplitudes\n1. the 10-90 percentile ranges of log amplitudes (more robust than range)\n1. an estimate/proxy for the energy level of the signal, using the sum of the log amplitudes"},{"metadata":{"trusted":true,"_uuid":"dd7b7c73d5592e54c90f3c3c2bcfe63a19ee0b2c"},"cell_type":"code","source":"RAW_SEQ_LEN = 150000\n# integer factors: {20, 25, 30, 40, 50, 60, 75, 80, 100, 120, 150}\nDSAMP_SEQ_LEN = 40\n\nassert (RAW_SEQ_LEN / DSAMP_SEQ_LEN) % 1 == 0\ndownsample_rate = int(RAW_SEQ_LEN / DSAMP_SEQ_LEN)\n\ndef range_downsample(x, dsamp_rate):\n    mat = x.reshape(-1, dsamp_rate).T\n    return mat.max(axis=0) - mat.min(axis=0)\n\ndef range80_downsample(x, dsamp_rate):\n    mat = x.reshape(-1, dsamp_rate).T\n    q90 = np.quantile(mat, 0.9, axis=0)\n    q10 = np.quantile(mat, 0.1, axis=0)\n    return q90 - q10\n\ndef energy_downsample(x, dsamp_rate):\n    mat = x.reshape(-1, dsamp_rate).T\n    return np.fabs(mat).sum(axis=0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"110fb4540584b1a912adde9734a35b4a3c822d10"},"cell_type":"markdown","source":"Calculate the mean and standard deviation of the features in order to standardize inputs to model.  We will reframe the full dataset to generate these scaling factors."},{"metadata":{"trusted":true,"_uuid":"1e1ee75fedbb3db6f037d82711ed2a259431b43b"},"cell_type":"code","source":"def get_scaling_factors(data, dsamp_rate):\n    trunc = (len(data) // dsamp_rate) * dsamp_rate\n    mat = data[:trunc].reshape((-1, dsamp_rate)).T\n\n    min_ = mat.min(axis=0)\n    max_ = mat.max(axis=0)\n    range_mean = (max_ - min_).mean()\n    range_std = (max_ - min_).std()\n\n    q10 = np.quantile(mat, 0.1, axis=0)\n    q90 = np.quantile(mat, 0.9, axis=0)\n    range80_mean = (q90 - q10).mean()\n    range80_std = (q90 - q10).std()\n\n    energy_mean = np.fabs(mat).sum(axis=0).mean()\n    energy_std = np.fabs(mat).sum(axis=0).std()\n\n    return (range_mean, range80_mean, energy_mean,\n            range_std, range80_std, energy_std)\n\nrange_mean, range80_mean, energy_mean,\\\n    range_std, range80_std, energy_std =\\\n    get_scaling_factors(obs, downsample_rate)\n\nprint(f'downsampled sequence length: {DSAMP_SEQ_LEN}')\nprint(f'downsample rate: {downsample_rate}')\nprint(f'range mean: {range_mean: .4f}, ' +\n      f'range sigma: {range_std: .4f}')\nprint(f'range80 mean: {range80_mean: .4f}, ' +\n      f'range80 sigma: {range80_std: .4f}')      \nprint(f'energy mean: {energy_mean: .4f}, ' +\n      f'energy sigma: {energy_std: .4f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"414a5863a117cc690f8ab57c0e57f4392076716e"},"cell_type":"markdown","source":"Define a function to retrieve sample features from a sample of raw signal, and a function to generate batches of these samples."},{"metadata":{"trusted":true,"_uuid":"0ec0c22ec4d3a2014be063f09d54c7ccd0e4742d"},"cell_type":"code","source":"def get_feats(args_in):\n    sample_in, dsamp_rate = args_in\n    range_ = ((range_downsample(sample_in, dsamp_rate)\n                - range_mean)\n                / range_std)\n    range80 = ((range80_downsample(sample_in, dsamp_rate)\n                - range80_mean)\n               / range80_std)\n    energy = ((energy_downsample(sample_in, dsamp_rate)\n               - energy_mean)\n              / energy_std)\n    return np.c_[range_, range80, energy]\n\n\ndef get_batch(data_in, n_samples, dsamp_rate, seq_len=150000):\n    n_rows = data_in.shape[0]\n    start_list = [rng.randint(0, n_rows-seq_len)\n                  for i in range(n_samples)]\n    X = map(get_feats,\n            ((data_in[start:start+seq_len, 0], dsamp_rate)\n                 for start in start_list))\n    y = [data_in[start+seq_len, 1] for start in start_list]\n    return list(X), y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f4f5647feab37e46f0891a7b80c0b98019770ef5"},"cell_type":"markdown","source":"Compare the downsampled features with the original log signal.  A lot of information is lost, but the coarse levels of the features correspond to amplitude patterns in the original signal.  Also, the features are very highly correlated.  Either the 10-90 range or the energy would probably be sufficient as a single input to the RNN, but we will use all three features in this demonstration."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"b44fc70f59f090d5572288e328265d6566e7c944"},"cell_type":"code","source":"start = 55000000\nsig_in = data[start:start+RAW_SEQ_LEN, 0]\nfeats = get_feats((sig_in, downsample_rate))\n\ndata_list = [sig_in] + [feats[:, i] for i in range(3)]\ntitles = [_ + ':' for _ in ['original log signal', 'range',\n                            '10-90 percentile range', 'energy']\n         ]\n\nfig = plt.figure(figsize=(14, 8))\nfor i in range(4):\n    ax = fig.add_subplot(4, 1, i+1)\n    plt.plot(data_list[i])\n    plt.title(titles[i], fontsize=12, loc='left')\n    ax.margins(x=0)\n    plt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"09287b421456b7bd47fe642850cdfdadc9e54393"},"cell_type":"markdown","source":"Next, set model parameters.  The primary options are:\n1. Number of layers in the RNN\n1. Number of hidden dimensions in each RNN layer\n1. Forward-only RNN, or bi-directional\n1. Type of RNN cell ([LSTM](https://arxiv.org/abs/1402.1128) or [GRU](https://arxiv.org/abs/1412.3555))\n\nAfter the RNN, we concatenate the final hidden states from each layer into a single 1-dimensional vector, and feed that vector into a densely-connected network.  Using all of the final hidden states preserves multiple levels of extracted features (fine-grained signal information and larger-scale features) to pass into the dense layers.\n\nIn the case of a bi-directional RNN, concat the final states from the forward and reverse RNNs into a single vector."},{"metadata":{"trusted":true,"_uuid":"d24cd12a9f47aae68224a94d69239496d37c81cd"},"cell_type":"code","source":"n_features = 3\n\ninit_mode = 'FAN_AVG'\ninit_uniform = True # False for normal init\n\n# RNN config:\ndim_h = [256, 256, 256] # dimensions of RNN internal layers\nmodel_type = 'bi'  # {'bi', 'fwd'}: bidirectional or forward-only\ncell_type = 'LSTM' # {'LSTM', 'GRU'}\nrnn_activation = tf.nn.elu\nlstm_use_peepholes = True\n\n# feedforward layers:\ndim_ff = [512, 512, 512] # dimensions of feedforward layers\nffn_activation = tf.nn.relu\nffn_drop_rate = 0.2\n\nn_epochs = 2000\neta = 0.000005\ntrain_batch_size = 100\nobjective = 'MSE'  # {'MSE', 'MAE', 'exp_error'}\noptimizer_name = 'adam'# ['sgd', 'adam', 'adagrad', 'adadelta']\neval_batch_size = 100\nn_eval_batches = 10\nstdout_eval_interval = 100\ntensorboard_eval_interval = 20\n\noptimizers = {\n    'sgd': tf.train.GradientDescentOptimizer,\n    'adam': tf.train.AdamOptimizer,\n    'adagrad': tf.train.AdagradOptimizer,\n    'adadelta': tf.train.AdadeltaOptimizer\n    }","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"00d697512bc5aff76665cd3f16961cee75e1b9ff"},"cell_type":"markdown","source":"Build the computational graph."},{"metadata":{"trusted":true,"_uuid":"f937d3edde6fd944db8e6fb882cb17673dd417df"},"cell_type":"code","source":"try:\n    sess.close()\nexcept: pass\n\nglobal graph\ngraph = tf.Graph()\ntf.reset_default_graph()\nwith graph.as_default():\n\n    X = tf.placeholder(tf.float32,\n        shape=(None, DSAMP_SEQ_LEN, n_features), name='X')\n    y = tf.placeholder(tf.float32, shape=(None), name='y')\n    X_r = tf.reverse(X, [1], name='X_reversed')\n    train_flag = tf.placeholder_with_default(False, shape=(), name='train_flag')\n\n    init = tf.contrib.layers.variance_scaling_initializer(\n        mode=init_mode, uniform=init_uniform)\n\n    if cell_type == 'LSTM':\n        cell = partial(tf.nn.rnn_cell.LSTMCell,\n                       use_peepholes=lstm_use_peepholes,\n                       activation=rnn_activation,\n                       initializer=init)\n\n    elif cell_type == 'GRU':\n        cell = partial(tf.nn.rnn_cell.GRUCell,\n                       activation=rnn_activation,\n                       kernel_initializer=init)\n\n    # unroll RNN\n    with tf.name_scope('forward_RNN'):\n        fwd_cells = [cell(num_units=layer, name='l' + str(i+1) + '_fwd')\n                    for i, layer in enumerate(dim_h)]\n\n    with tf.name_scope('reverse_RNN'):\n        rev_cells = [cell(num_units=layer, name='l' + str(i+1) + '_rev')\n                    for i, layer in enumerate(dim_h)]\n\n    multi_rnn = tf.contrib.rnn.MultiRNNCell(fwd_cells)\n    multi_rnn_r = tf.contrib.rnn.MultiRNNCell(rev_cells)\n\n    rnn_out, states = tf.nn.dynamic_rnn(multi_rnn, X, dtype=tf.float32)\n    rnn_out_r, states_r = tf.nn.dynamic_rnn(multi_rnn_r, X_r, dtype=tf.float32)\n\n    # concat the final states of each RNN layer\n    # rnn_out = batch x seq_len x dim_h\n    # states = (layers x (c, h)) x batch x dim_h\n    concat_fwd = tf.concat([layer[1] for layer in states], 1)\n    concat_rev = tf.concat([layer[1] for layer in states_r], 1)\n\n    # concat the forward and reverse outputs (if bidirectional)\n    if model_type == 'bi':\n        concat = tf.concat([concat_fwd, concat_rev], 1)\n    elif model_type == 'fwd':\n        concat = tf.identity(concat_fwd)\n\n    # unroll feedforward layers\n    ff_layers = [concat]\n    dropout = partial(tf.layers.dropout, rate=ffn_drop_rate, training=train_flag)\n    for i, layer in enumerate(dim_ff):\n        ff_layers.append(dropout(tf.layers.dense(ff_layers[i], dim_ff[i],\n                                                 activation=ffn_activation,\n                                                 name='feedforward_' + str(i+1))))\n\n    outputs = tf.layers.dense(ff_layers[-1], 1, name='outputs')\n\n    MAE = tf.reduce_mean(tf.abs(outputs - y))\n    MSE = tf.reduce_mean(tf.square(outputs - y))\n    exp_error = tf.reduce_mean(tf.subtract(tf.exp(tf.abs(outputs - y)), 1))\n\n    optimizer_fn = optimizers[optimizer_name]\n\n    if objective == 'MAE':\n        training_op = optimizer_fn(learning_rate=eta).minimize(MAE)\n    elif objective == 'MSE':\n        training_op = optimizer_fn(learning_rate=eta).minimize(MSE)\n    elif objective == 'exp_error':\n        training_op = optimizer_fn(learning_rate=eta).minimize(exp_error)\n\n    with tf.name_scope('MAE'):\n        train_mae = tf.summary.scalar('train', MAE)\n        val_mae = tf.summary.scalar('val', MAE)\n    with tf.name_scope('MSE'):\n        train_mse = tf.summary.scalar('train', MSE)\n        val_mse = tf.summary.scalar('val', MSE)\n\n    init = tf.global_variables_initializer()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8bde3a33cfcddcc42312209c624e21adbdd5df17"},"cell_type":"markdown","source":"Train and evaluate the model."},{"metadata":{"trusted":true,"_uuid":"c4e47c6102c1bf38717414ad63097c4f057c64a0"},"cell_type":"code","source":"sess = tf.InteractiveSession(graph=graph)\ninit.run(session=sess)\n\nfor epoch in tnrange(n_epochs):\n    X_train, y_train = get_batch(\n        train_data, train_batch_size, downsample_rate, RAW_SEQ_LEN)\n    sess.run(training_op, feed_dict={X:X_train, y:y_train, train_flag:True})\n    if (epoch+1) % stdout_eval_interval == 0:\n        evals_out = {'train': {'mse': [], 'mae': []}, 'eval': {'mse': [], 'mae': []}}\n        for i in range(n_eval_batches):\n            X_train_batch, y_train_batch = get_batch(\n                train_data, train_batch_size, downsample_rate, RAW_SEQ_LEN)\n            evals_out['train']['mse'].append(MSE.eval(\n                feed_dict={X: X_train_batch, y: y_train_batch, train_flag:False}))\n            evals_out['train']['mae'].append(MAE.eval(\n                feed_dict={X: X_train_batch, y: y_train_batch, train_flag:False}))\n            X_val_batch, y_val_batch = get_batch(\n                val_data, eval_batch_size, downsample_rate, RAW_SEQ_LEN)\n            evals_out['eval']['mse'].append(MSE.eval(\n                feed_dict={X: X_val_batch, y: y_val_batch, train_flag:False}))\n            evals_out['eval']['mae'].append(MAE.eval(\n                feed_dict={X: X_val_batch, y: y_val_batch, train_flag:False}))\n        print(f\"round {epoch+1} \" +\n              f\"train MSE: {sum(evals_out['train']['mse'])/n_eval_batches:.4f} \" +\n              f\"train MAE: {sum(evals_out['train']['mae'])/n_eval_batches:.4f} \" +\n              f\"val MSE: {sum(evals_out['eval']['mse'])/n_eval_batches:.4f} \" +\n              f\"val MAE: {sum(evals_out['eval']['mae'])/n_eval_batches:.4f} \")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ac363e77bf013a19b445c85c2239e3719ad353a4"},"cell_type":"code","source":"test_path = os.path.join(data_path, 'test')\nfiles = os.listdir(test_path)\nprint('total files', len(files)) # 2624\npreds = {'seg_id': [], 'time_to_failure': []}\nfor fname in tqdm_notebook(files[:10]):\n    path = os.path.join(test_path, fname)\n    obs = pd.read_csv(path).values.ravel()\n    obs = preprocess_obs(obs, obs_mean)\n    feats = get_feats((obs, downsample_rate)).reshape((1, DSAMP_SEQ_LEN, n_features))\n    pred = float(outputs.eval(feed_dict={X:feats}))\n    preds['time_to_failure'].append(pred)\n    preds['seg_id'].append(fname.split('.')[0])\npreds_df = pd.DataFrame(preds)\npreds_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"53055c572ad795b5877399a6b2304b8db55184bc"},"cell_type":"code","source":"preds_df.to_csv('submissions.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}