{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"## **Can convolution extract useful features from a spectrogram of seismic measurements?**\n\n*[See [here](https://www.kaggle.com/michael422/tf-lstm-on-acoustic-signals) for a related article on using [Recurrent Neural Networks](https://www.kaggle.com/michael422/tf-lstm-on-acoustic-signals) on seismic data.]*\n\n[Spectrograms](https://en.wikipedia.org/wiki/Spectrogram) evaluate a 1-dimensional signal input, and produce a 2-dimensional output which is a [Fourier transform](https://en.wikipedia.org/wiki/Fourier_transform) moving through time.  The approach is to divide the signal into frames, then apply a discrete Fourier transform ([FFT](https://en.wikipedia.org/wiki/Fast_Fourier_transform)) to each of the frames using a smoothing kernel or window function.  The output indicates the power of the signal in each of a fixed number of frequency bins (example below).\n\n1-dimensional convolutional neural networks are often used in signal processing.  2-dimensional convolution can also be used with spectrograms ([example](https://arxiv.org/abs/1701.02720)).  In the context of speech recognition, a variety of additional preprocessing techniques suited to the audio spectrum of human speech are applied.  Here, we will attempt convolution with minimal preprocessing on log spectrograms generated from the LANL labquake data.\n\nWe use the first 400MM records in order to process efficiently in-memory:"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true,"scrolled":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\nfrom scipy import signal\nimport matplotlib.pyplot as plt\nnp.warnings.filterwarnings('ignore')\n\nrng = np.random.RandomState(datetime.now().microsecond)\n\ndata_path = '../input/LANL-Earthquake-Prediction'\ntf_path = os.path.join(data_path, 'tf')\n\nnrows = 400000000   # 629145480\ndata = (pd.read_csv(os.path.join(data_path, 'train.csv'), nrows=nrows,\n                   dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\n            .values)\nprint(f'input data shape: {data.shape}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"abc0ee665f6327f9571b1934ac2c90d206996aa7"},"cell_type":"markdown","source":"First, downsample the training data.  This will aid performance and reduce some of the noise in the raw signal.  The [scipy.signal.decimate](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.decimate.html) function applies a low-pass (noise-reducing) filter, and samples the data every n-th value based on the downsample rate"},{"metadata":{"trusted":true,"_uuid":"4b08d83b36e3e3748f386dbc00296ad2e1c1c70e","_kg_hide-input":true},"cell_type":"code","source":"def downsample(data, downsample_rate):\n    trunc = int((data.shape[0] // downsample_rate) * downsample_rate)\n    data = data[:trunc, :]\n    x = signal.decimate(data[:, 0].ravel(), downsample_rate)\n    y_ = data[:, 1].reshape((-1, downsample_rate)).T\n    y = y_.mean(axis=0)\n    return x, y\n\nRAW_SEQ_LEN = 150000\nDOWNSAMPLE_RATE = 4\n\nassert (RAW_SEQ_LEN / DOWNSAMPLE_RATE) % 1 == 0\ndownsampled_seq_len = int(RAW_SEQ_LEN/DOWNSAMPLE_RATE)\n\nX_train_, y_train_ = downsample(data, DOWNSAMPLE_RATE)\n\nprint(f'downsample rate: {DOWNSAMPLE_RATE}')\nprint(f'downsampled seq len: {downsampled_seq_len}')\nprint(f'X shape: {X_train_.shape} y shape: {y_train_.shape}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"47c33cc2347b39419782f003a6d75ba7f27ec2c8"},"cell_type":"markdown","source":"Next, define a spectrogram function using [scipy.signal.spectrogram](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html).  The choice of window function, segment length, and overlap all affect the output.  In this case, we use a triangular smoothing window, a frame length of 256, and an overlap of 64.  The function will return a log spectrogram, standardized by arguments we will determine next.\n\nFor more background on applications of spectrograms to machine learning, see Haytham Fayek's [blog](https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html) post."},{"metadata":{"trusted":true,"_uuid":"3ee93522afc36d250d7b0f04a244081e976f59cc"},"cell_type":"code","source":"def spectrogram(sig_in, dsamp):\n    nperseg = 256 # default 256\n    noverlap = nperseg // 4 # default: nperseg // 8\n    fs = 4000000 // dsamp # raw signal sample rate is 4MHz\n    window = 'triang'\n    scaling = 'density' # {'density', 'spectrum'}\n    detrend = 'linear' # {'linear', 'constant', False}\n    eps = 1e-11\n    f, t, Sxx = signal.spectrogram(sig_in, nperseg=nperseg, noverlap=noverlap,\n                                   fs=fs, window=window,\n                                   scaling=scaling, detrend=detrend)\n    return f, t, np.log(Sxx + eps)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"69633f569083fc528016330134c9f4d437bacd73"},"cell_type":"markdown","source":"Visualize a sample signal and the log spectrogram generated from it.  The vertical axis is a range of frequency bands ranging 0-500kHz.  The axis is inverted so that high frequencies are at the top."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"da769f9e79985b4e9d45b0db11e8e45c31490e97"},"cell_type":"code","source":"start = rng.randint(0, X_train_.shape[0] - downsampled_seq_len)\nsig_in = X_train_[start:start+downsampled_seq_len]\n\nf, t, Sxx_out = spectrogram(sig_in, DOWNSAMPLE_RATE)\n\nfig = plt.figure(figsize=(18, 8))\nax = fig.add_subplot(2, 1, 1)\nax.margins(x=0.003)\nplt.plot(X_train_[start:start+downsampled_seq_len])\nplt.title('downsampled signal:', fontsize=18, loc='left')\nplt.axis('off')\n\nax = fig.add_subplot(2, 1, 2)\ncmap = plt.get_cmap('magma')\nspec = plt.pcolormesh(t, f, Sxx_out, cmap=cmap)\nplt.title('normalized log spectrogram (aspect stretched):',\n          fontsize=18, loc='left')\nplt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3a47963413537cd307315dff0914add7798c8a6f"},"cell_type":"markdown","source":"Using convolution on spectrographs is appealing on the surface, since spectrograms are 2-dimensional and visually intuitive.  However, the height and width of a photograph measure the same thing: spatial distance (3-dimensional distance projected to 2 dimensions).  2D convolution leverages the notion that objects in a photograph can shift, rotate, and warp (due to perspective), but the object is still the same.  Spectrograms are different.  The horizontal axis measures time, while the vertical axis corresponds to frequency.  See [this article](https://towardsdatascience.com/whats-wrong-with-spectrograms-and-cnns-for-audio-processing-311377d7ccd) for a more detailed discussion.\n\nIn fact, a spectrogram is a just time series of frequency measurements.  Local features (periodic, repeating signals) are present in most time series on multiple scales.  Convolution is one way to extract these recurring patterns if they exist.\n\nFollowing (https://arxiv.org/abs/1701.02720) we will augment the features with a delta spectrogram (representing the change in FFT response between time frames) and the delta-delta spectrogram, which is the second-order change. Since these 2D matrices are all indexed by the same frequency bands and time periods, we will stack the matrices along a third dimension, creating features that are h x w x 3."},{"metadata":{"trusted":true,"_uuid":"518661d4c7b4a9cd346e6ea8bdfbd60d9b03ee3a"},"cell_type":"code","source":"N_CHANNELS = 3\n\ndef get_3d_spec(Sxx_in, moments=None):\n    if moments is not None:\n        (base_mean, base_std, delta_mean, delta_std,\n             delta2_mean, delta2_std) = moments\n    else:\n        base_mean, delta_mean, delta2_mean = (0, 0, 0)\n        base_std, delta_std, delta2_std = (1, 1, 1)\n    h, w = Sxx_in.shape\n    right1 = np.concatenate([Sxx_in[:, 0].reshape((h, -1)), Sxx_in], axis=1)[:, :-1]\n    delta = (Sxx_in - right1)[:, 1:]\n    delta_pad = delta[:, 0].reshape((h, -1))\n    delta = np.concatenate([delta_pad, delta], axis=1)\n    right2 = np.concatenate([delta[:, 0].reshape((h, -1)), delta], axis=1)[:, :-1]\n    delta2 = (delta - right2)[:, 1:]\n    delta2_pad = delta2[:, 0].reshape((h, -1))\n    delta2 = np.concatenate([delta2_pad, delta2], axis=1)\n    base = (Sxx_in - base_mean) / base_std\n    delta = (delta - delta_mean) / delta_std\n    delta2 = (delta2 - delta2_mean) / delta2_std\n    stacked = [arr.reshape((h, w, 1)) for arr in (base, delta, delta2)]\n    return np.concatenate(stacked, axis=2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"60a75294fd91486d59123784d034143592ce5dd9"},"cell_type":"markdown","source":"Sample the training data and process spectrograms for the samples in order to get factors for standardizing inputs to the model.  There are three pairs of moments required to standardized the inputs: a mean and sigma for the spectrograms as well as the deltas and delta-deltas:"},{"metadata":{"trusted":true,"_uuid":"51027a0001afd24f17c55460bf601c7a61e5357d","_kg_hide-input":true},"cell_type":"code","source":"def get_moments(X, downsampled_seq_len, samp=1000):\n    nrows = X.shape[0]\n    start_list = [rng.randint(0, nrows-downsampled_seq_len)\n                  for i in range(samp)]\n    Sxx_samples = [spectrogram(X[start:start+downsampled_seq_len],\n                               dsamp=DOWNSAMPLE_RATE)[2]\n                   for start in start_list]\n    sxx_h, sxx_w = Sxx_samples[0].shape\n    Sxx_3d_samples = np.array([get_3d_spec(Sxx) for Sxx in Sxx_samples])\n    base_mean = Sxx_3d_samples[:, :, :, 0].mean()\n    base_std = Sxx_3d_samples[:, :, :, 0].std()\n    delta_mean = Sxx_3d_samples[:, :, :, 1].mean()\n    delta_std = Sxx_3d_samples[:, :, :, 1].std()\n    delta2_mean = Sxx_3d_samples[:, :, :, 2].mean()\n    delta2_std = Sxx_3d_samples[:, :, :, 2].std()\n    return (sxx_h, sxx_w), (base_mean, base_std, delta_mean,\n                            delta_std, delta2_mean, delta2_std)\n\n(sxx_h, sxx_w), moments = get_moments(X_train_, downsampled_seq_len)\nprint(f'spectrogram dims: {sxx_h}x{sxx_w}')\nprint(f'base spectrogram mean, sigma: {moments[0]:.4f}, {moments[1]:.4f}')\nprint(f'delta spectrogram mean, sigma: {moments[2]:.4f}, {moments[3]:.4f}')\nprint(f'delta-delta spectrogram mean, sigma: {moments[4]:.4f}, {moments[5]:.4f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b98c9582d294a999c34e5ba0b9f62ad1e4f4fa3d"},"cell_type":"markdown","source":"Visualize a sample  spectrogram as well as the delta and delta-delta spectrograms:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"d3fcf94775b827ac4d3ec6bf2eca0dd8c6471c24"},"cell_type":"code","source":"start = 21307218\nsig_in = X_train_[start:start+downsampled_seq_len]\n\nf, t, Sxx = spectrogram(sig_in, DOWNSAMPLE_RATE)\ns3d = get_3d_spec(Sxx, moments)\n\nfig = plt.figure(figsize=(18, 16))\nax = fig.add_subplot(4, 1, 1)\nax.margins(x=0.003)\nplt.plot(X_train_[start:start+downsampled_seq_len])\nplt.title('downsampled signal:', fontsize=14, loc='left')\nplt.axis('off')\n\nax = fig.add_subplot(4, 1, 2)\ncmap = plt.get_cmap('magma')\nspec = plt.pcolormesh(t, f, s3d[:, :, 0], cmap=cmap)\nplt.title('log spectrogram (aspect stretched):',\n          fontsize=14, loc='left')\nplt.axis('off')\n\nax = fig.add_subplot(4, 1, 3)\nspec = plt.pcolormesh(t, f, s3d[:, :, 1], cmap=cmap)\nplt.title('delta log spectrogram:', fontsize=14, loc='left')\nplt.axis('off')\n\nax = fig.add_subplot(4, 1, 4)\nspec = plt.pcolormesh(t, f, s3d[:, :, 2], cmap=cmap)\nplt.title('delta-delta log spectrogram:', fontsize=14, loc='left')\nplt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a39f8259daf4bc693b199764018db951e3e83e85"},"cell_type":"markdown","source":"The delta spectrograms add useful features that emphasize abrupt changes in the signal at different frequencies.\n\nIn this application, we do not use filter banks, which transform the frequency band of the spectrogram to a log scale, leaving constant spacing between harmonic frequencies.  Filter banks are frequently used in analysis of speech and other audible signals.  However, the signals in this dataset are mostly well above 50kHz (see below), while human hearing is limited to about 20kHz, so Mel-frequency banks are not suitable.  A custom log-scale filter bank might bring out the harmonics in this particular spectrum, but we will keep things simple for now."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"d34828e4ef673381534d7874bf78325cf289c0d8"},"cell_type":"code","source":"plt.figure(figsize=(18, 4))\ncmap = plt.get_cmap('magma')\nax = plt.pcolormesh(t, f/1000, s3d[:, :, 0], cmap=cmap)\nplt.ylabel(r'frequency (kHz)', fontsize=14)\nax.axes.tick_params(axis='y', labelsize=14)\nax.axes.get_xaxis().set_visible(False);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0edbb256fe6b278c4b940679b6be37768c0b39f0"},"cell_type":"markdown","source":"Split the data into train and validation sets.  The train and validation data have different distributions, which will lead to unexpected validation results below."},{"metadata":{"trusted":true,"_uuid":"060d759dc3ddd8ac59894af6ddee9684b6ced3bd"},"cell_type":"code","source":"TRAIN_VAL_RATIO = 0.7\n\ntrain_size = int(X_train_.shape[0] * TRAIN_VAL_RATIO)\nX_train = X_train_[:train_size]\nX_val = X_train_[train_size:]\ny_train = y_train_[:train_size]\ny_val = y_train_[train_size:]\n\nprint(f'train mean ttf: {y_train.mean():.4f}, val mean ttf: {y_val.mean():.4f}')\nprint(f'train 1st quartile ttf: {np.quantile(y_train, 0.25):.4f}, ' +\n      f'val 1st quartile ttf: {np.quantile(y_val, 0.25):.4f}')\nprint(f'train 3rd quartile ttf: {np.quantile(y_train, 0.75):.4f}, ' +\n      f'val 3rd quartile ttf: {np.quantile(y_val, 0.75):.4f}')\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9844fe51c690e4ffb805bbd0615377e52c445fc3"},"cell_type":"markdown","source":"Define a batching function to sample the input signal.  2d convolution takes a rank-4 tensor as input:\n\n    (batch x image height x image width x channels)\nImage files have three color channels.  The spectrograms stacked with their deltas and delta-deltas also have a final dim of 3."},{"metadata":{"trusted":true,"_uuid":"ff808072ed1c369c9e1b66ed7cb178944999df8f"},"cell_type":"code","source":"def get_batch(X, y, n_samples, seq_len, moments, sxx_h, sxx_w, dsamp_rate):\n    start_list = [rng.randint(0, X.shape[0]-seq_len) for i in range(n_samples)]\n    X_batch = []\n    for start in start_list:\n        sig_in = X[start:start+seq_len]\n        _, _, Sxx = spectrogram(sig_in, dsamp_rate)\n        \n        Sxx3d = get_3d_spec(Sxx, moments)\n        X_batch.append(Sxx3d)\n    n_channels = X_batch[0].shape[-1]\n    X_batch = np.array(X_batch).reshape(n_samples, sxx_h, sxx_w, n_channels)\n    y_batch = [y[start+seq_len-1] for start in start_list]\n    return X_batch, y_batch","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6c325257011be25421fb8113f04317dfe448f337"},"cell_type":"markdown","source":"Nex, build an inner layer for the convnet.  We will use three layers of convolution, followed by a pooling layer.  This inner layer includes options for several regularization and optimization features:\n\n1. [Batch normalization](https://arxiv.org/abs/1502.03167).  This aids learning in deeper neural networks.  We will not use batch norm in this kernel.\n\n1. [Local response regularization](http://papers.nips.cc/paper/4824-imagenet-classification-with-deep-convolutional-neural-networks) (LRN):  This layer encourages adjacent convolutional kernels to train on different features, but it introduces another set of hyperparameters to be tuned.\n\n1. Residual connections: This allows information from the inputs to reach the pooling layer with only one layer of convolution, using a \"skip\" layer.  The skip layer is added element-wise to the outputs of the final convolutional sub-layer.  This is particularly helpful when there are many iterations of the interior layer.  Residual connections are used in [ResNet](https://arxiv.org/abs/1512.03385) and in many other deep learning models.\n\nFor the purpose of demonstration, we will use LRN on the pooling layers only.  We will also enable residuals.  In practice, these elements  were designed for use with large, deep image classification networks, and may not be suitable for this particular learning task.  In fact, even though the data from the spectrogram is visually informative, its structure is very different from the photographs that convnets are typically trained on."},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"5d7ecc7cd5e3fdd3fcf0ee30f67a1cc248be8d0b"},"cell_type":"code","source":"\ndef layer(inputs, filt, kn, stride, act, ptype, psize, pstride, ppad,\n          use_lrn_conv, use_lrn_pool, use_bn=False, use_res=False, name=None):\n    '''Internal layer composed of 2x conv1d, then pool.\n    Args:\n        inputs: tf.Tensor rank 4\n        filt: int, number of maps to use\n        kn: list of tuples (h, w), kernel size\n        stride: list of tuples (h, w), stride for each conv layer\n        act: tf.nn object or None, type of activation after each conv layer\n        ptype: str in {'max', 'avg'}: type of pooling layer\n        psize: tuple (h, w), size of pooling layer\n        pstride: tuple (h, w), stride of pooling layer\n        ppad: str in {'same', 'valid'}, padding type for pooling layer\n        use_lrn: bool, optional local response normalization\n        use_bn: bool, optional batch norm\n        use_res: bool, optional residual connection around conv layers\n        name: str, name for the layer\n\n    Returns: tf.Tensor rank 4, layer output\n    '''\n    with graph.as_default():\n        batch_norm = partial(tf.layers.batch_normalization, training=train_flag,\n                     momentum=bn_momentum)\n        lrn = partial(tf.nn.local_response_normalization, depth_radius=5, bias=1,\n                      alpha=1, beta=0.5)\n\n        conv1 = conv2d(inputs, filters=filt, kernel_size=kn[0],\n                       strides=stride[0], name=name+'_conv1')\n        if use_lrn_conv[0]:\n            conv1 = lrn(conv1, name=name+'_conv1_lrn')\n        if use_bn: conv1 = batch_norm(conv1, name=name+'_conv1_bn')\n        if act is not None:\n            conv1 = act(conv1, name=name+'_conv1_act')\n\n        conv2 = conv2d(conv1, filters=filt, kernel_size=kn[1],\n                       strides=stride[1], name=name+'_conv2')\n        if use_lrn_conv[1]:\n            conv2 = lrn(conv2, name=name+'_conv2_lrn')\n        if use_bn: conv2 = batch_norm(conv2, name=name+'_conv2_bn')\n\n        if use_res:\n            stride_h = np.prod([i[0] for i in stride])\n            stride_w = np.prod([i[1] for i in stride])\n            kn_h = np.prod([i[0] for i in kn])\n            kn_w = np.prod([i[1] for i in kn])\n            resid = conv2d(inputs, filters=filt, kernel_size=(kn_h, kn_w),\n                           strides=(stride_h, stride_w), name=name+'_skip')\n\n            conv_final = tf.add(conv2, resid, name=name+'_resid')\n        else:\n            conv_final = tf.identity(conv2, name=name+'_conv_final')\n\n        if act is not None:\n            conv_final = act(conv_final, name+'_conv2_act')\n\n        if ptype == 'max':\n            pool = tf.layers.max_pooling2d(conv_final, pool_size=psize,\n                        strides=pstride, padding=ppad, name=name+'_max_pool')\n        elif ptype == 'avg':\n            pool = tf.layers.average_pooling2d(conv_final, pool_size=psize,\n                        strides=pstride, padding=ppad, name=name+'_avg_pool')\n        else:\n            pool = tf.identity(conv_final, name=name+'_pool_passthru')\n        \n        if use_lrn_pool:\n            pool = lrn(pool, name=name+'_pool_lrn')\n\n        return pool","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f708b39d32157c1473ec9c882348da2d2eaa1e9a"},"cell_type":"markdown","source":"A graph of the internal layer (from Tensorboard) makes this a bit clearer.  The two nodes between each convolutional layer are the LRN layer and the activation:"},{"metadata":{"trusted":true,"_uuid":"24bdac085dc33ed759cb1a6b9dca905f3fee38f8","_kg_hide-input":true},"cell_type":"code","source":"from IPython.display import Image\nImage(filename=r'../input/images/graph.JPG', height=300, width=300)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"532681c5be4d271481417be08430677fa84f5e61"},"cell_type":"markdown","source":"Next, set hyperparameters for the model.  For information about convnet design, see this [Stanford class](https://cs231n.github.io/convolutional-networks/):  For more about implementing convnets in TensorFlow, see [here](https://www.tensorflow.org/tutorials/estimators/cnn).\n\nConvolutional neural nets have a lot of hyperparameters.  Kernel, stride, and pool dimensions all have a big impact on the quality of the model.\n\nWe use mean squared error (MSE) as the objective function.  Mean absolute error (MAE) is also an option, but MSE more heavily penalizes predictions that are far from the actual value.  (Symmetric) exponential error is another option, but can lead to overfitting because outliers have an outsized effect on training loss."},{"metadata":{"trusted":true,"_uuid":"2a02b7f4c9ff849ef2560518541718c2fcfdf0d1"},"cell_type":"code","source":"init_mode = 'FAN_AVG'\ninit_uniform = True\n\n# convolutional layer(s):\nn_layers = 2\nfilters=[64, 128]\nkernels = [[(3, 6), (4, 8)],\n           [(3, 6), (4, 8)]\n          ]\nstrides = [[(1, 2), (2, 2)],\n           [(1, 2), (2, 2)]\n          ]\nconv_activation = tf.nn.elu\nuse_conv_bn = [0, 0]\nbn_momentum = 0.9999 # only if use_conv_bn\nuse_lrn_conv = [[0, 0, 0],\n                [0, 0, 0]\n               ]\nuse_lrn_pool = [1, 1]\nuse_residuals = [1, 1]\nptypes = ['max', 'avg']\nppads = ['valid', 'valid']\npools = [(5, 5,), (5, 5,)]\npstrides = [(1, 1), (1, 1)]\n\n# feedforward layer:\ndim_ff = [1024]\nff_activation = tf.nn.relu\ndrop_rate = [0.4]\n\n# train config:\nn_epochs = 2500\neta = 2e-8\ntrain_batch_size = 50\neval_batch_size = 50\nn_eval_batches = 10\noptimizer_name = 'adam'  # ['sgd', 'adam', 'adagrad', 'adadelta']\nobjective = 'MSE'  # ['MSE', 'MAE', 'exp_error']\ntensorboard_eval_interval = 20\nstdout_eval_interval = 100\n\noptimizers = {\n    'sgd': tf.train.GradientDescentOptimizer,\n    'adam': tf.train.AdamOptimizer,\n    'adagrad': tf.train.AdagradOptimizer,\n    'adadelta': tf.train.AdadeltaOptimizer\n    }\n\ninit = tf.contrib.layers.variance_scaling_initializer(mode=init_mode,\n                                                      uniform=init_uniform)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f3c53b5298bdfbf327de8911abbe1e132c8a961"},"cell_type":"markdown","source":"Now construct the model graph.  This particular model has the following structure:\n* one or more of the inner layers defined above, each of which has 2x convolutional sub-layers followed by a pooling layer.\n* a Flatten() layer to collapse the 2-dimensional outputs from convolution to a 1-d vector\n* a [dropout](http://jmlr.org/papers/v15/srivastava14a.html) layer for regularization\n* 1x feedforward layer\n* a linear layer to produce a scalar prediction for the target value (time to failure).\n\nThis structure is verly loosely based on features of AlexNet and ResNet, and is not tailored to this particular task."},{"metadata":{"trusted":true,"_uuid":"c2b1eba120161cb4234bcf0c5587e9badd24ecc6"},"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, sxx_h, sxx_w, N_CHANNELS),\n                       name='X')\n    y = tf.placeholder(tf.float32, shape=(None), name='y')\n    train_flag = tf.placeholder_with_default(False, shape=(), name='train_flag')\n\n    conv2d = partial(tf.layers.conv2d, padding='same',\n                     activation=conv_activation,\n                     kernel_initializer=init)\n\n    lrn = partial(tf.nn.local_response_normalization,\n                  depth_radius=5, bias=1, alpha=1, beta=0.5)\n    \n    # unroll internal (3x conv1d+pool) layers\n    layers = [X]\n    for i in range(n_layers):\n        layers.append(\n            layer(layers[i], filt=filters[i], kn=kernels[i], stride=strides[i],\n                  act=conv_activation, ptype=ptypes[i], psize=pools[i],\n                  pstride=pstrides[i], ppad=ppads[i],\n                  use_lrn_conv=use_lrn_conv[i], use_lrn_pool=use_lrn_pool[i],\n                  use_bn=use_conv_bn[i], use_res=use_residuals[i], name='L'+str(i+1)))\n\n    # feedforward layer(s)\n    dense = partial(tf.layers.dense,\n                    activation=ff_activation,\n                    kernel_initializer=init)\n\n    dropout = partial(tf.layers.dropout,\n                      training=train_flag)\n    \n    flat = tf.layers.Flatten()(layers[-1])\n    ff_layers = [flat]\n    for i, width in enumerate(dim_ff):\n        drop_layer = dropout(ff_layers[-1], rate=drop_rate[i], name='drop_' + str(i+1))\n        ff_layer = dense(drop_layer, width, name='ff_' + str(i+1))\n        ff_layers.append(ff_layer)\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    init = tf.global_variables_initializer()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"744974d9024e06fff50f1d73d1e5bd76bf0efe86"},"cell_type":"markdown","source":"Train and evaluate the model."},{"metadata":{"trusted":true,"_uuid":"d26b1037ed98078c9ea4ec43cd277f20e108e06c"},"cell_type":"code","source":"sess = tf.InteractiveSession(graph=graph)\ninit.run(session=sess)\nfor epoch in tnrange(n_epochs):\n    X_batch, y_batch = get_batch(X_train, y_train, train_batch_size, downsampled_seq_len,\n                                 moments, sxx_h, sxx_w, DOWNSAMPLE_RATE)\n    sess.run(training_op, feed_dict={X: X_batch, y: y_batch, 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(X_train, y_train, train_batch_size,\n                    downsampled_seq_len, moments, sxx_h, sxx_w, DOWNSAMPLE_RATE)\n            evals_out['train']['mse'].append(\n                MSE.eval(feed_dict={X: X_train_batch, y: y_train_batch, train_flag:False}))\n            evals_out['train']['mae'].append(\n                MAE.eval(feed_dict={X: X_train_batch, y: y_train_batch, train_flag:False}))\n            X_val_batch, y_val_batch = get_batch(X_val, y_val, eval_batch_size,\n                    downsampled_seq_len, moments, sxx_h, sxx_w, DOWNSAMPLE_RATE)\n            evals_out['eval']['mse'].append(\n                MSE.eval(feed_dict={X: X_val_batch, y: y_val_batch, train_flag:False}))\n            evals_out['eval']['mae'].append(\n                MAE.eval(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":"c644c117e3db74ce56ba3c5b1b26814911d31281"},"cell_type":"markdown","source":"Uh oh, the model trains for ~1700 its, but then learning plateaus.  A cursory examination of the outputs shows that the model is severaly underfitting.  It's probably worth exploring different convolutional structures that are better suited to this task, plus quite a bit of hyperparameter tuning."},{"metadata":{"trusted":true,"_uuid":"ec30a95ac28863d23d708729c8931d851d303c67"},"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):\n    path = os.path.join(test_path, fname)\n    obs = np.array(pd.read_csv(path).values)\n    data = np.c_[obs, np.zeros(len(obs))]\n    x, y = downsample(data, DOWNSAMPLE_RATE)\n    _, _, Sxx = spectrogram(x, DOWNSAMPLE_RATE)\n    Sxx3d = get_3d_spec(Sxx, moments)\n    x_feed = Sxx3d.reshape((1, sxx_h, sxx_w, N_CHANNELS))\n    pred = float(outputs.eval(feed_dict={X:x_feed}))\n    preds['time_to_failure'].append(pred)\n    preds['seg_id'].append(fname.split('.')[0])\npreds_df = pd.DataFrame(preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"be2beaab373d303580b21f894a170ace06dd430e"},"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}