{"cells":[{"metadata":{"_uuid":"eb7a9f0e9cea01af54767003eaf78388885ab10d"},"cell_type":"markdown","source":"# LANL Earthquake Prediction\nA *Short-Time Fourier Transform* analysis introduction"},{"metadata":{"_uuid":"e2339054b9fbeb44bc1a89a88077fe89df940fdb"},"cell_type":"markdown","source":"-----------------------"},{"metadata":{"_uuid":"f4d1d1e5bd638c29cee580804589f1201bcd0c42"},"cell_type":"markdown","source":"## 1. Introduction\nIn this notebook I seek to apply a simplified EDA on the available training data. Going further, and since earthquake data are strongly related to cyclic events, I intend to explore this data on the Fourier domain. "},{"metadata":{"_kg_hide-input":true,"_uuid":"e314f7ce5bd1641e379e736789e896ec768c01d8","trusted":true},"cell_type":"code","source":"### Importing libraries\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.colors as colors\nimport matplotlib.pyplot as plt\n\n%matplotlib inline\n\nimport scipy.signal as signal\nimport os\nimport gc","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8b4a8c49433f5046d90f6867bc25020fd892b1d4"},"cell_type":"markdown","source":"## 2. Analyzing the data"},{"metadata":{"_uuid":"080380bb0483858db7df553034bc9565f9bd9b76","scrolled":true,"trusted":true},"cell_type":"code","source":"## Listing the files in the directory\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1b4dd13fcf152146079568366057bcea44c918c8"},"cell_type":"markdown","source":"In order to make it easy to deal with, we will start reading just a small portion of the data. To do so we will use the *skiprows* and *nrows* parameters:"},{"metadata":{"_uuid":"e6fcb927f8a42ba28c27f19127b6855a9efe3883","trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv', skiprows=range(1,4000000), nrows=2000000, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f4523a3ae6ecdc6899363937baa931c062090ed3"},"cell_type":"markdown","source":"Notice I had to skip the first 4M samples due to memory overuse here in Kaggle's cloud environment. Running on my local computer, it's better setting the *skiprows* parameter with the value 0, so the time-frequency analysis will be shown with the right time scale. In the next versions of this kernel a short function can be defined to treat this correctly. "},{"metadata":{"_uuid":"d3997f234dfdc94c6a57a1ed0aa0e94aba24ee18","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"## Checking the number of samples in the file\nprint(\"Dataset size:\", train.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"82a19d7e8b66ee3494d3625422441c6f27e4246a","trusted":true},"cell_type":"code","source":"## Showing the first data samples with high precision\npd.set_option(\"display.precision\", 10)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6c6c7c42dd216687ae3f9f53a2eb1fd18858b9e1"},"cell_type":"markdown","source":"An important aspect of the data, when dealing with time-frequency analysis, is the sampling rate. Accordingly to the [additional information](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77526), *the input is a chunk of 0.0375 seconds of seismic data (ordered in time), which is recorded at 4MHz, hence 150'000 data points, and the output is time remaining until the following lab earthquake, in seconds.*  "},{"metadata":{"_uuid":"4e46560350a752e660af8824a7bb447f8496965a"},"cell_type":"markdown","source":"Given the sampling frequency ($4{MHz}$) and the number of samples ($N=1e7$), we can estimate the time length of the signal:"},{"metadata":{"_uuid":"d9bcbc642f9d8e07827ef81907bc50d0101f0ab6","trusted":true},"cell_type":"code","source":"(1/4e6)*1e7","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c919a822c62735ca2de4fb6b4fc8839c1efe9869"},"cell_type":"markdown","source":"Notice that since we are skipping the first $4e6$ samples, we have a resulting sinal length of 0.5s that is equivalent to the interval from 1.0 to 1.5 seconds of the original signal."},{"metadata":{"_uuid":"a4ed5af674af8fc0c67086ab93fa0f8d1a6ed4e9"},"cell_type":"markdown","source":"\n### 2.1 Short-Time Fourier Transform"},{"metadata":{"_uuid":"0dbac90b9234e588596e8c53130d892b11246b27"},"cell_type":"markdown","source":"The following function will be used to calculate the *Short-Time Fourier Transform* for a given signal."},{"metadata":{"_uuid":"6d48adc13e9cef3e75074295a5d39507147b3d3c","trusted":true},"cell_type":"code","source":"def calcSTFT_norm(inputSignal, samplingFreq, window='hann', nperseg=256, figsize=(9,5), cmap='magma', ylim_max=None):\n    '''Calculates the STFT for a time series:\n        inputSignal: numpy array for the signal (it also works for Pandas.Series);\n        samplingFreq: the sampling frequency;\n        window : str or tuple or array_like, optional\n            Desired window to use. If `window` is a string or tuple, it is\n            passed to `get_window` to generate the window values, which are\n            DFT-even by default. See `get_window` for a list of windows and\n            required parameters. If `window` is array_like it will be used\n            directly as the window and its length must be nperseg. Defaults\n            to a Hann window.\n        nperseg : int, optional\n            Length of each segment. Defaults to 256.\n        figsize: the plot size, set as (6,3) by default;\n        cmap: the color map, set as the divergence Red-Yellow-Green by default;\n        ylim_max: the max frequency to be shown. By default it's the half sampling frequency.'''\n    ##Calculating STFT\n    f, t, Zxx = signal.stft(inputSignal, samplingFreq, window=window, nperseg=nperseg)\n    ##Plotting STFT\n    fig = plt.figure(figsize=figsize)\n    ### Different methods can be chosen for normalization: PowerNorm; LogNorm; SymLogNorm.\n    ### Reference: https://matplotlib.org/tutorials/colors/colormapnorms.html\n    spec = plt.pcolormesh(t, f, np.abs(Zxx), \n                          norm=colors.PowerNorm(gamma=1./8.),\n                          #norm=colors.LogNorm(vmin=np.abs(Zxx).min(), vmax=np.abs(Zxx).max()),\n                          #norm=colors.SymLogNorm(linthresh=0.13, linscale=1,\n                          #                       vmin=-1.0, vmax=1.0),\n                          cmap=plt.get_cmap(cmap))\n    cbar = plt.colorbar(spec)\n    ##Plot adjustments\n    plt.title('STFT Spectrogram')\n    ax = fig.axes[0]\n    ax.grid(True)\n    ax.set_title('STFT Magnitude')\n    if ylim_max:\n        ax.set_ylim(0,ylim_max)\n    ax.set_ylabel('Frequency [Hz]')\n    ax.set_xlabel('Time [sec]')\n    fig.show\n    return","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"becee1f23a17957a1c314081bdc35076d515650d"},"cell_type":"markdown","source":"There's no specific rule to determine the number of samples per segment. However, it's useful to consider 1/4 of the sampling rate. E.g. for a 1kHz sampling rate, the nperseg parameter is set by default to 256. If we consider 4MHz, we can estimate the ideal nperseg by $log_2(4e6 / 4)$."},{"metadata":{"_uuid":"8d2c2bcdf4f6c3146fa28f5ecbc437442777be9e","trusted":true},"cell_type":"code","source":"calcSTFT_norm(train['acoustic_data'], 4e6, nperseg=1048576, ylim_max=300000)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"daad7b038548e6dcc2e635c9031917bafb76abe4"},"cell_type":"markdown","source":"From the plot above, we can see there's a sharp signal (i.e. goes to a high frequency range) around 1.2 seconds of the signal."},{"metadata":{"_uuid":"a8bd1e3b7c87ae9269d9349a9af9e9b3a6d594e4"},"cell_type":"markdown","source":"## Time-frequency analysis in near-failure periods:"},{"metadata":{"_uuid":"155fe1c31dd1294ffa5fb07b005282b0340aaea0"},"cell_type":"markdown","source":"As stated in the previous version of this kernel, we cannot affirm there is any relation between the pulses in the frequency domain and the seismic activities: the event shown in time 0.2s (or 1.2s) could be due to a measurement error, for example.  \n\nIn order to obtain more solid evidences of this hypothetical relationship, I propose the following test: select the positions in the dataset where the *time_to_failure* attribute is close to null (i.e. where the failure is occurring is a local minimum) and check the domain frequency for a given interval before it."},{"metadata":{"_uuid":"50e852845cb4626b23d5fb8cc68545151d846949"},"cell_type":"markdown","source":"### Finding the failure positions and regions\nSince we are working with a huge dataset (at least when running the kernel on Kaggle's cloud), the most efficient way to process the data is reading it in chunks, from where we will look for the local minima, or the \"failure time\". Then, for each failure time in each chunk we will extract the **failure region**, defined as the 1.0s before the failure and the 0.25s after the failure. The region limits are calculated by the sampling frequency, which results that *1e6* samples are equivalent to 0.25s."},{"metadata":{"trusted":true,"_uuid":"1c1990654f0d3127f35f4837d25a32f4c04f533e"},"cell_type":"code","source":"## Creating a list of positions where the \"time_to_failure\" is null\nfailure_regions = {}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3d1dc3c124781bc4aae3485521a8a03e5162cb89"},"cell_type":"code","source":"## Defining a function to find the local minima positions\ndef is_failure(df):\n    failures = df.index[(df.time_to_failure.shift(1) > df.time_to_failure) & (df.time_to_failure.shift(-1) > df.time_to_failure)].tolist()\n    return failures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0b513679829ae1b542b563f1560e470159c83622"},"cell_type":"code","source":"## Returns a signal region for a given failure position in a pandas.chunk:\n## REMEMBER: a chunk is a pandas.Dataframe.\ndef region_readChunk(pos, chunk):\n    ## Setting the region limits:\n    beforePos = pos - 4e6\n    afterPos = pos + 1e6\n    #In some cases it can be shorter than the other regions\n    if beforePos < chunk.index.min():\n        beforePos = chunk.index.min()\n    if afterPos > chunk.index.max():\n        afterPos = chunk.index.max()\n    ## Extracting the region data:\n    data = chunk.loc[int(beforePos):int(afterPos)]\n    return data","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"e2809e6847f6d0b8702158e1767718dff7b72344"},"cell_type":"code","source":"## Processing the dataset in chunks:\n## For each chunk, find the local minima and extract the respective region.\nfor chunk in pd.read_csv('../input/train.csv', chunksize=5e7, \n                         dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64}):\n    failure_pos = is_failure(chunk)\n    for pos in failure_pos:\n        data = region_readChunk(pos, chunk)\n        failure_regions[pos] = data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2156ab17eaaf5c616368f28a40b3947b953d7d56"},"cell_type":"code","source":"## Garbage collection (to reduce used memory)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7c623632a0b6b162e9fcbf623457885da45da362"},"cell_type":"markdown","source":"Once the data wrangling process is finished, we then obtain a dictionary (it could be stored in any other data structure) where the keys are the failure positions and the items correspond to the failure regions. "},{"metadata":{"trusted":true,"_uuid":"1c8ec46e768ec765f3c0a44e6806c690f341d55e"},"cell_type":"code","source":"failure_regions.keys()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9c95a2d6e8bace2c25eaf9c557354e88d3a6f28e"},"cell_type":"markdown","source":"### Choosing the right chunk size\nI initialy set the chunksize attribute as 5e6. However, this lead to many failure positions starting almost in the beginning of the chunk. To choose it more appropriately, I will now check if is there a regular space among them:"},{"metadata":{"trusted":true,"_uuid":"87abf6bf28bcb914972b769fed9c6711adb66385"},"cell_type":"code","source":"posFailures = pd.DataFrame(list(failure_regions.keys()))\nposFailures.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bed44deb2add1d48adbe5d1b3150199f67f6e540"},"cell_type":"markdown","source":"We can notice that the failure positions are around 3e8 with a standard deviation about 2e8. From this we can choose 5e8 as an appropriate chunk size."},{"metadata":{"_uuid":"1a691ff934b6ac764fef6b88a1c809333af4137c"},"cell_type":"markdown","source":"### Analyzing the failure regions"},{"metadata":{"trusted":true,"_uuid":"44b91451c6c202c6b9aeafc58fe218a9cc5a4cdd"},"cell_type":"code","source":"def region_plot(df):\n    data = df.copy()\n    ##Aplying a 1e3 gain to the \"time_to_failure\" in order to make the failure visible\n    data['time_to_failure'] = data['time_to_failure']*100\n    ##Calculate elapsed time through the index and sampling frequency (4e6 MHz)\n    data['time'] = data.index\n    data['time'] = data['time']*(1/4e6)\n    ##Normalize time information (0.0s  to 1.25s)\n    data['Time [sec]'] = data['time'] - data['time'].min()\n    ##Plot lines\n    data[['acoustic_data','time_to_failure','Time [sec]']].plot(x='Time [sec]', figsize=(8,5))\n    return","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4534a6b78f5b664a7230fed74c468d09c1ba8fe3"},"cell_type":"markdown","source":"The code above could be use to analyze the signal region of each failure point. However, since it uses too much memory, we will choose a random one:"},{"metadata":{"scrolled":false,"trusted":true,"_uuid":"c486d51c4ec9cbd9ef133f5f4833184bd236b37a"},"cell_type":"code","source":"#for key, value in failure_regions.items():\n#    region_plot(value)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cec5f273d87cda124737d0283c9c4381ce66acdb"},"cell_type":"code","source":"region_plot(failure_regions[419368879])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b41c649350053106ce5310f0f7a509ae16dda0ad"},"cell_type":"code","source":"calcSTFT_norm(failure_regions[419368879]['acoustic_data'], 4e6, nperseg=1048576, ylim_max=300000)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6022588de78444d75097bd3cea3cbd0acf1d859a"},"cell_type":"markdown","source":"## Discussion and the next steps"},{"metadata":{"_uuid":"1c00cb25b2c6a0c664a96098f8242c20a0ec5d74"},"cell_type":"markdown","source":"In the previous version of this kernel we could not affirm that a time-frequency analysis would be useful for this problem. However, since now we extend it to all the failure points, there are some evidences that time-frequency analysis can be useful for interpreting or even for feature extraction."},{"metadata":{"trusted":false,"_uuid":"5d6f7b02832fa7901123cd89c4b1b38906f4b13f"},"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.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}