{"cells":[{"metadata":{},"cell_type":"markdown","source":"# LANL Earthquakes Competition: why do spectrograms fail?"},{"metadata":{},"cell_type":"markdown","source":"## Summary\n\nAs I went over the first lectures of FastAI Deep-Learning course, and as I've previously encountered the [use of Convolutional NN for analysis of spectrograms of audio signals](https://towardsdatascience.com/audio-classification-using-fastai-and-on-the-fly-frequency-transforms-4dbe1b540f89), I decided to apply this approach on the [LANL earthquakes competition](https://www.kaggle.com/c/LANL-Earthquake-Prediction/overview).\n\nSeveral people already referred to this approach in the context of the LANL competition (e.g. [[1](https://www.kaggle.com/michael422/spectrogram-convolution)],[[2](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/91776#latest-529254)]), but I did not notice any successful reported results (which of course may simply indicate my ignorance).\n\n**I didn't manage to build any impressive spectrogram-based CNN model either, and I'd like to try to shed some light on the reasons, focusing on mel-spectrograms which are based on librosa library.**\n\n**In short, it seems that librosa applies by default certain manipulations on the spectrograms (e.g. choice of frequency-scaling and power-scaling), which might have made sense for audio classification, but are apparently less relevant to the LANL competition.**\n\n**Tuning librosa's spectrogram methods to correspond to LANL dataset seems quite tricky, thus I believe that any attempt of spectrograms in the scope of the competition would benefit from avoiding such general audio tools - unless the user knows exactly how these tools work and why.**\n\n**Having that said, since the different spectrograms variants that I tried led to similar and not-so-good prediction results, I currently don't see spectrogram-based models as a promising direction.\nI'd be glad to hear any observations and thoughts regarding this issue.**\n\nFor whoever interested in the librosa-based spectrograms of the train data nonetheless, the code that generates them appears at the bottom of this notebook (took me ~40 minutes to run on my laptop, but for some reason I can't read through the end of train.csv within this kernel)."},{"metadata":{},"cell_type":"markdown","source":"#### Update (21.5.19)\nI tried to generate manual spectrograms. They look reasonable (BW images with the color directly representing the power of the signal). However, the results haven't changed too much, so currently the whole direction doesn't look promising to me."},{"metadata":{},"cell_type":"markdown","source":"## Initialization"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"%reload_ext autoreload\n%autoreload 2\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"import numpy as np\nimport scipy as sp\nimport scipy.stats as stats\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport seaborn as sns\nimport re\nfrom pathlib import Path\nimport pickle as pkl\nfrom time import time\nfrom tqdm import tqdm, tnrange, tqdm_notebook\nfrom pprint import pprint\nimport os, sys\nfrom warnings import warn\nimport itertools\nimport librosa\nimport librosa.display","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BASE = Path('../input')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Mel-spectrogram tuning"},{"metadata":{},"cell_type":"markdown","source":"Following several resources (e.g. [this one](https://medium.com/@etown/great-results-on-audio-classification-with-fastai-library-ccaf906c5f52)), I decided to **use librosa library** to generate the images.\n\nAfter several attempts, I decided (according to my judgement of the images' visuality, see below) to go with **log-scaled frequencies and log-scaled powers**.\nFollowing the same references mentioned above, I decided to do so through librosa's **mel-spectrograms**.\nThe main steps are simply applying *librosa.feature.melspectrogram()*, *librosa.amplitude_to_db()* and *librosa.display.specshow()* sequentially on the input signal."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Load some data segments\n\ndef get_segment(skiprows=None, samples_per_segment = 150000):\n    df = pd.read_csv(BASE/'train.csv', nrows=samples_per_segment, skiprows=skiprows,\n                 dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n    df.rename(columns={'acoustic_data': 'signal', 'time_to_failure': 'quaketime'}, inplace=True)\n    return df\n\n\ndf = get_segment()\ndf2 = get_segment(range(1,int(2.2e6)))\ndf_peak = get_segment(range(1,int(4.4e6)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_spec(df, n_fft=4096, hop_length=1024, sr=4096e3, n_mels=64, fmin=1480, fmax=640e3):\n    stft = librosa.stft(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length)\n    stft_magnitude, stft_phase = librosa.magphase(stft)\n    stft_magnitude_db = librosa.amplitude_to_db(stft_magnitude)\n\n    fig,axs = plt.subplots(2,2, figsize=(12,8))\n    librosa.display.specshow(stft_magnitude, ax=axs[0,0], x_axis='time', y_axis='linear', \n                             sr=sr, hop_length=hop_length, fmin=fmin, fmax=fmax)\n    axs[0,0].set_ylim((fmin,fmax))\n    axs[0,0].set_title(f'Spectrogram ({stft_magnitude_db.shape}, hop={hop_length:d}, Linear)')\n\n    librosa.display.specshow(stft_magnitude_db, ax=axs[0,1], x_axis='time', y_axis='linear', \n                             sr=sr, hop_length=hop_length, fmin=fmin, fmax=fmax)\n    axs[0,1].set_ylim((fmin,fmax))\n    axs[0,1].set_title(f'Spectrogram ({stft_magnitude_db.shape}, hop={hop_length:d}, DB)')\n\n    mel_spec = librosa.feature.melspectrogram(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length,\n                                              n_mels=n_mels, sr=sr, power=1.0,\n                                              fmin=fmin, fmax=fmax)\n    librosa.display.specshow(mel_spec, ax=axs[1,0], x_axis='time', y_axis='mel', \n                             sr=sr, hop_length=hop_length, fmin=fmin, fmax=fmax)\n    axs[1,0].set_title(f'Spectrogram ({mel_spec.shape}, hop={hop_length:d}, Linear)')\n\n    mel_spec = librosa.feature.melspectrogram(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length,\n                                              n_mels=n_mels, sr=sr, power=1.0,\n                                              fmin=fmin, fmax=fmax)\n    mel_spec = librosa.amplitude_to_db(mel_spec, ref=np.max)\n    librosa.display.specshow(mel_spec, ax=axs[1,1], x_axis='time', y_axis='mel', \n                             sr=sr, hop_length=hop_length, fmin=fmin, fmax=fmax)\n    axs[1,1].set_title(f'Spectrogram ({mel_spec.shape}, hop={hop_length:d}, DB)')\n    \n    plt.tight_layout()\n\nplot_spec(df)\nplot_spec(df_peak)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfirtunately, the pretty logarithmic spectrograms suffered from **annoying black rows in their images**.\n\nIt turns out that mel-spectrograms essentially sample the FFT of the signal on multiple bands of frequencies (***filter banks***), which become exponentially-wider in large frequencies. In addition, **the banks of mel-spectrograms use [universally-conventional values](https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html)**, which were probably chosen for applications of humanly-heardable audio sounds, apparently unlike the frequencies of a simulated earthquakethe.\nThe frequencies mis-match resulted in empty frequency-banks which were represented as black rows in the images.\n\nAs a result, I had to either change the whole banks scale (which I didn't easily find a trivial way to do through librosa library) or use fictive sample rate. I chose the latter, and instead of the true sample rate of ~4MHz, I fed librosa (using again my special visual judgment skills) with a fictive sample rate of 44KHz."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_spec(dfs, n_fft=4096, hop_length=None, sr=44e3, n_mels=256, fmin=0, fmax=128e3, power=2.0, lin=False):\n    if hop_length is None:\n        hop_length = n_fft // 4\n    fig,axs = plt.subplots(1+lin,len(dfs), figsize=(6*len(dfs),4+4*lin))\n    for i,df in enumerate(dfs):\n        ax = axs[0,i] if lin else axs[i]\n        mel_spec = librosa.feature.melspectrogram(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length,\n                                                  n_mels=n_mels, sr=sr, power=power)#, fmin=fmin, fmax=fmax)\n        mel_spec = librosa.amplitude_to_db(mel_spec, ref=np.max)\n        l = mel_spec.shape[0]\n        mel_spec = mel_spec[int(0.02*l):int(0.75*l),:]\n        librosa.display.specshow(mel_spec, ax=ax, x_axis='time', y_axis='mel', \n                                 sr=sr, hop_length=hop_length)#, fmin=fmin, fmax=fmax)\n        ax.set_title(f'Spectrogram ({mel_spec.shape}, hop={hop_length:d})')\n        \n        if lin:\n            ax = axs[1,i]\n            stft = librosa.stft(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length)\n            stft_magnitude, stft_phase = librosa.magphase(stft)\n            stft_magnitude_db = librosa.amplitude_to_db(stft_magnitude)\n            librosa.display.specshow(stft_magnitude_db, ax=ax, x_axis='time', y_axis='linear', \n                                     sr=sr, hop_length=hop_length, fmin=fmin, fmax=fmax)\n            ax.set_ylim((fmin,fmax))\n            ax.set_title(f'Spectrogram ({stft_magnitude_db.shape}, hop={hop_length:d})')\n\n\nplot_spec((df,df2,df_peak), sr=44e3, n_fft=4096, n_mels=256, power=2.0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This time, the problem was **similar scale of colors for all the spectrograms, despite of significant differences in the powers of the underlying signals**, which should carry important imformation regarding the desired time-to-failure (which is the number we need to predict in the competition).\n\nAfter going through the code, I found out that\n1. I used amplitude_to_db(..., ref=np.max), which chooses the DB scale according to the maximum of the current signal's FFT, hence I lost most of the differentiation between segments' scales.\n2. I used plt.imsave(filename, image) which uses an input array image instead of plt.savefig(filename) which uses the current open figure, hence I totally ignored the plotting function specshow().\n\nI tried to choose a constant reference to the DB scale, plot nicely-looking spectrograms with specshow(), and save them correctly with savefig().\n**I did receive a much more variat color scale** (see below), and I did recieve (over 20 randomly chosen segments) a **quite monotonous (though apparently not continuous) relation between the power of the signal and the power of the colors** in the spectrograms."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from scipy.misc import imread\nseg_len = int(150e3)\nn_fft=4096\nhop_length=1024\nsr=44e3\nn_mels=256\npower=2.0\nclass_limits=(0,1,2,3,4,5,6,7,8,10,12,99)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Get 20 random segments + a peak segment, generate their spectrogram-arrays and PNG files.\n\nnp.random.seed(0)\nsigs = []\nspecs = []\nmels = []\nims = []\nhops = list(np.random.randint(0,6000,20)) + [int(4.4e6/hop_length)]\nfor hop in hops:\n    # load signal\n    tmp = pd.read_csv(BASE/'train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32},\n                      nrows=seg_len, skiprows=range(1,1+hop*hop_length))\n    tmp.rename(columns={'acoustic_data': 'signal', 'time_to_failure': 'quaketime'}, inplace=True)\n    # extract signal\n    sigs.append(tmp.signal.values)\n    # generate spec\n    mel_spec = librosa.feature.melspectrogram(np.array(tmp.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length,\n                                          n_mels=n_mels, sr=sr, power=power)\n    specs.append(mel_spec)\n    mel_spec = librosa.amplitude_to_db(mel_spec, ref=10**8)\n    l = mel_spec.shape[0]\n    mel_spec = mel_spec[int(0.02*l):int(0.75*l),:]\n    mels.append(mel_spec)\n    # save as PNG and load image\n    fig = plt.figure(figsize=(2.6,2.6))\n    fig.subplots_adjust(left=0,right=1,bottom=0,top=1)\n    librosa.display.specshow(mel_spec, ax=plt.axes(), x_axis='time', y_axis='mel', vmin=-120, vmax=-20, # -50,50\n                             sr=sr, hop_length=hop_length)\n    plt.axes().set_axis_off()\n    plt.savefig('tmp_spec.png')\n    plt.close()\n    ims.append(imread('tmp_spec.png'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Look for relation between the signal \"power\" (measured by mean(|signal|))\n# and the PNG image color range (measured by sum(color^2) for each channel separately).\n\nfig, axs = plt.subplots(1,3, figsize=(15,5))\nfor channel in range(3):\n    axs[channel].plot([np.mean(np.abs(m.flatten())) for m in sigs],\n                      [np.mean(([t**2 for t in im[:,:,channel].flatten()])) for im in ims],'.')\n    axs[channel].set_xlabel('mean |signal|')\n    axs[channel].set_ylabel('sum(color^2)')\n    axs[channel].set_title(f'Channel {channel:d}')\n    axs[channel].set_xscale('log')\n    axs[channel].set_yscale('log')\n    axs[channel].grid()\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, axs = plt.subplots(7,3, figsize=(12,18))\n\nfor i,(s,im) in enumerate(zip(sigs,ims)):\n    ax = axs[i%7,i//7]\n    ax.imshow(im)\n    ax.set_axis_off()\n    ax.set_title(f'mean(|signal|) = {np.mean(np.abs(s.flatten())):.1f}')\n\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately, all the efforts **did not yield much success in predicting the next earthquake**.\n\nI tried to build a standard RESNET34 CNN using FastAI to predict the TTF of each segment using its spectrogram. The results on my Cross-Validation mechanism were quite bad, and I didn't bother applying the model on public test data to submit for LB (according to the empirically-known relation between my CV-score and the LB-score, I'd probably get around 1.6 on LB)."},{"metadata":{},"cell_type":"markdown","source":"**In conclusion, it looks like blackbox-based spectrograms (e.g. the mel-spectrograms of librosa library) are quite tricky to use for the LANL dataset, since they apply manipulations on the signal which are often unwanted, and are not trivial to understand for the common user.**\n\n**While it is not guaranteed that a spectrogram-based model can achieve good results on the LANL dataset, I'd focus any further attempts on manually-generated spectrograms, with careful tuning of scales of frequencies and powers.**"},{"metadata":{},"cell_type":"markdown","source":"## Generate spectrograms data"},{"metadata":{},"cell_type":"markdown","source":"If anyone still happens to be interested, the code below should generate a dataset of the latter variant of the spectrograms (as PNG files) - just make sure to adjust the paths below.\n\n(took about 40 minutes on my laptop)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"BASE = Path('../input')\nTRAIN = Path('../specs')\n\nn_all = 629145480 # df_full.shape[0]\nseg_len = int(150e3)\nhop_len = int(25e3)\nn_hops = int(n_all/hop_len)\nsegs_to_read = 1050\nhops_to_read = int(segs_to_read * seg_len/hop_len)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def save_spec(df, seg_index, meta, class_limits=(0,1,2,3,4,5,6,7,8,10,12,99), base_path=TRAIN,\n              n_fft=4096, hop_length=1024, sr=44e3, n_mels=256, power=2.0, fig=None):\n    \n    # get ttf and derive the filename\n    tf = df.quaketime.values[-1]\n    cls = np.where([a<=tf<b for a,b in zip(class_limits[:-1],class_limits[1:])])[0][0]\n    cls_nm = '-'.join((f'{class_limits[cls]:02d}',f'{class_limits[cls+1]:02d}'))\n    fname = str(seg_index) + '_' + cls_nm + '.png'\n    \n    # compute spectrogram\n    mel_spec = librosa.feature.melspectrogram(np.array(df.signal,dtype='float'), n_fft=n_fft, hop_length=hop_length,\n                                              n_mels=n_mels, sr=sr, power=power)\n    mel_spec = librosa.amplitude_to_db(mel_spec, ref=10**8)\n    l = mel_spec.shape[0]\n    mel_spec = mel_spec[int(0.02*l):int(0.75*l),:]\n    meta.loc[seg_index] = [fname, tf, cls_nm]\n    \n    # save as PNG\n    if fig is None:\n        fig = plt.figure(figsize=(2.6,2.6))\n    fig.subplots_adjust(left=0,right=1,bottom=0,top=1)\n    librosa.display.specshow(mel_spec, ax=plt.axes(), x_axis='time', y_axis='mel', vmin=-50, vmax=50,\n                             sr=sr, hop_length=hop_length)\n    plt.axes().set_axis_off()\n    plt.savefig(base_path/fname)\n    plt.clf()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Read train data in 4 chunks of 1K segments;\nsave each segment's spectrogram (as png file);\nand keep a table of segments and their quaketimes:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"GENERATE_SPECTS = False\n\nif GENERATE_SPECTS:\n    meta = pd.DataFrame(columns=('filename','time','class'))\n    fig = plt.figure(figsize=(2.6,2.6))\n    for i in tqdm_notebook(range(n_hops)):\n        if i % hops_to_read == 0:\n            tmp = pd.read_csv(BASE/'train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32},\n                              nrows=seg_len*(segs_to_read+1), skiprows=range(1,1+i*hop_len))\n            tmp.rename(columns={'acoustic_data': 'signal', 'time_to_failure': 'quaketime'}, inplace=True)\n        save_spec(tmp[(i%hops_to_read)*hop_len:(i%hops_to_read)*hop_len+seg_len], i, meta, fig=fig)\n\n    meta.to_csv(TRAIN/'train_spec_meta.csv', index=False)\n    \n    plt.figure(figsize=(10,6))\n    plt.plot(meta.time)\n    plt.xlabel('Training segment')\n    plt.ylabel('TTF')\n    meta.head()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"fastai v1","language":"python","name":"Python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.8"}},"nbformat":4,"nbformat_minor":1}