{"cells":[{"metadata":{"_uuid":"d615566253915a33500c554b96925869e56458ca"},"cell_type":"markdown","source":"## Predicting Earthquakes?\nPredicting earthquakes has long been thought to be near-impossible. But same has been said about many other things, I still remember how 15 years ago my IT teacher in high school was saying that speech recognition can never be done reliably by a computer, and today it's more reliable than a human. So who knows, maybe this competition will be a major step in solving another \"impossible\" problem. Importance of predicting earthquakes is hard to underestimate - selected years of our century claimed hundreds of thousands of causalities. Being able to predict earthquakes could allow us to better protect human life and property."},{"metadata":{"trusted":true,"_uuid":"c65dddbe77f899b5dcfdc9c627fb63edfaef59c7","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# please ignore this (some peculiarity of Kaggle virtualisation)\n!ln -s ../input/LANL-Earthquake-Prediction/sample_submission.csv ../input/sample_submission.csv\n!ln -s ../input/LANL-Earthquake-Prediction/test ../input/test\n!ln -s ../input/LANL-Earthquake-Prediction/train.csv ../input/train.csv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b820744598296f6b9998a8abfc588ba721fa7722","_kg_hide-input":true},"cell_type":"code","source":"# standard imports\n%matplotlib inline\nfrom matplotlib import pyplot as plt\nimport matplotlib as mpl\nimport pandas as pd\n# found this library when looking at this competiton - it proved worthy! :)\nimport librosa\nimport librosa.display\n\n# let's read 10M rown only for now (the dataset is HUGE! over 600M records)\ndf = pd.read_csv(\"../input/train.csv\", nrows=10_000_000)\n\n# let's rename columns to something easier to type: Signal, Time\ndf.columns=[\"s\",\"t\"]\n\n# let's identify indices of the \"earthquakes\"\n# (please look at the data description, if this doesn't make sense)\nquake_indices = df.index[df.t.diff()>0]\n\n# \"Fourier Transformations\" - just google that - it's a big, but fruitful topic\n# basically this allows us to \"see\" patterns of \"waves\" inside time-series data\nX = librosa.stft(df.s.values.astype(\"float\").clip(-50,50))\nXdb = librosa.amplitude_to_db(abs(X))\nplt.figure(figsize=(18, 6))\nax = librosa.display.specshow(Xdb, sr=1, x_axis=\"time\", y_axis=None)\nax.set_title('Spectrogram of the first 10M points of data provided - white line signifies the earthquake')\n\n# drawing white lines (just one), where earthquake occured\nfor quake in quake_indices:\n    plt.axvline(x=quake, color=\"white\")\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f04bb31cb262dbbabee8d9160f33683b07c63da0"},"cell_type":"markdown","source":"## Data Provided\nTraining data-set for this competition consists of a long time-series of measurements from some sort of an earthquake simulation in a laboratory. As a \"test\" dataset we are given bits of similar time-series, and the question we have to answer is \"when\". When the next earthquake (in a lab) is going to happen. The idea is that by analysing those time-series measurements and having some examples of when those lab earthquakes were registered in the past, we can make predictions about the future."},{"metadata":{"trusted":true,"_uuid":"2f8d933edc62469b7d771d49fbd7023db1eeba42","_kg_hide-input":true},"cell_type":"code","source":"# let's take every 200th line (including first), since \"train.csv\" is huge.\n# we store those lines in \"./train_every_Nth.csv\" for later reading with pandas\n!awk 'NR == 1 || NR % 200 == 0'  ../input/train.csv > ./train_every_Nth.csv\n\n# let's create DataFrame from those lines\ndf_n = pd.read_csv(\"./train_every_Nth.csv\", skiprows=1, names=[\"s\",\"t\"])\n\n# let's index every earthquake\nquake_indices = df_n.index[df_n.t.diff()>0]\nquake_indices\n\n# creating figure and clipping dataset\nplt.figure(figsize=(18, 6))\ndata = df_n.s.clip(-200,200)\n\n# a simple \"data.plot()\" would've suffice.. just to make it look a little bit better\nfor i in range(20):\n    data[i::20].plot(alpha=0.05,color=\"midnightblue\")\n\nplt.gca().set_title('Whole dataset - pink dashed lines are earthquakes')\n\n# drawing eartquake lines\nfor quake in quake_indices:\n    plt.axvline(x=quake, linestyle = \":\", color=\"salmon\")\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"08395081187dc53df3d5633dd8941531c9cf0f03"},"cell_type":"markdown","source":"## Hear it coming?\nWhen looking at the data in this form it looks suspiciously like sound waves. So my initial guess was to listen to it - to \"audiolise\" it, just like we're visualising graphical data. It was unknown if the earthquake can be predicted: \"by ear\". Just like in a dubstep song - the listener can hear the build-up before the \"drop\", so, perhaps, an earthquake can be felt in a similar fashion."},{"metadata":{"trusted":true,"_uuid":"fc132e913592df62cf9dae34ecd1811244e01f32","_kg_hide-input":true},"cell_type":"code","source":"# found these values by investigating tryint to include data around the first \"event\" \nSTART = 5_200_000\nEND = 5_900_000\n\n# taking an \"interesting\" sample of that data (one which includes an \"event\")\nsamp = df.iloc[START:END]\nquake_indices = df.index[df.t.diff()>0]\n\nplt.figure(figsize=(18, 6))\n\nax = plt.gca()\nax.set_title('Part of the raw data just around the earthquake - red line is when the earthquake occured')\n\n# lets plot sample, with event indicated\nfor i in range(1,10,1):\n    (samp.s*i/10).clip(-70,70).plot(alpha=(1-i/10), color=\"blue\")\n    plt.axvline(x=quake_indices[0]+1000*i, alpha=(1-i/10), color=\"red\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"40debf12d069cb0774817de79efe95af4a1c4da8"},"cell_type":"markdown","source":"## Is 44100 Hz enough?\nWhen interpreting signals as sounds, it's not quite obvious how one should go about it. Perhaps signal coming from experiment itself wasn't audible - for instance, it was infrasound. Either way, one approach was to just treat the data as a digital audio signal at 44100 Hz frequency. 44100 was chosen because it's most common audio resolution (a bit like 720p for video)."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-output":false,"_kg_hide-input":true},"cell_type":"code","source":"# we'll be using this library to convert dataframe into sound\nfrom scipy.io import wavfile\n# HTML tag to show audio player\nfrom IPython.core.display import HTML\n\n\n# constant for sampling frequency\nAUDIO_RATE = 44100\n# amplification constant - so the sound is a bit louder\nAMP_CONST = 400\n# adding \"volume channel\" to dataframe\ndf[\"v\"] = df.s * AMP_CONST\n\n# take a sub-sample from initial dataframe\nSTART = 5_200_000\nEND = 5_900_000\nsamp = df.iloc[START:END]\n\n# creating audiowave and writing into a file\nwave = (samp.v.values).astype(\"int16\")\nwavfile.write(\"./sound.wav\", AUDIO_RATE, wave)\n\n# displaying a player\nsrc = \"\"\"<audio controls=\"controls\" style=\"width:600px\" >\n    <source src=\"./sound.wav\" type=\"audio/wav\" />\n      Audio NOT working!</audio>\"\"\"\ndisplay(HTML(src))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4328af86c0d1c769f0b41a59189d703772e8dc21"},"cell_type":"markdown","source":"## De-noising the signal\nThe signal seemed to be full of noise, so one thing to try was to try to remove that noise by using 'sox' tool, which is conveniently provided by Kaggle image. The sound became more \"audible\", though it was still hard to make sense of what's happening!\n"},{"metadata":{"trusted":true,"_uuid":"3e95b7a9cda1a84aa43123993a9e9b53d7db4380","_kg_hide-input":true},"cell_type":"code","source":"# filename to use\nAUDIO = \"./sound\"\n\n# this is part where almost no \"clicking\" happens - we're using it as a guide to define what \"noise\" is\n!sox {AUDIO}.wav {AUDIO}_noise.wav trim 7.383 0.451\n# creating noise profile before removing\n!sox {AUDIO}_noise.wav -n noiseprof {AUDIO}_noise.prof\n# here we're remoiving the noise\n!sox {AUDIO}.wav {AUDIO}_clean.wav noisered {AUDIO}_noise.prof 0.15\n\nsrc = \"\"\"<audio controls=\"controls\" style=\"width:600px\" >\n    <source src=\"./sound_clean.wav\" type=\"audio/wav\" />\n      Audio NOT working!</audio>\"\"\"\ndisplay(HTML(src))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5ec2187d2bd55af5b517fff76a09afc74b4a5fac"},"cell_type":"markdown","source":"## Animating the wave\nWhat's lacking in the previous sample, is the ability to know when the earthquake was registered. In order to have that - we need to have some visualisation."},{"metadata":{"trusted":true,"_uuid":"c56edf73224af736acba96e69c7af40b7756703a","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"from matplotlib import lines\n# Fourier Transformations again\nfrom numpy.fft import rfft\nimport numpy as np\n\n# setting video rate to 24, for 'smooth' animation\nVIDEO_RATE = 24\n\n# setting up sub-plots, first one is for Fourier, second one for Amplitude\nfig, ax = plt.subplots(1, 2,figsize=(18,6))\n\n# earthquake locations\nquake_indices = samp.index[samp.t.diff()>0]\nmin_index = samp.index.min()\n\n# sliding line to show which part is playing\ndef init_blue_line(axi=ax):\n    global blue_line\n    axi.plot(samp.v, color=\"skyblue\")\n    for quake in quake_indices:\n        axi.axvline(x=quake, color=\"red\")\n    blue_line, = axi.plot([min_index, min_index], [0, 0], color=\"blue\")\n\n    return axi,\n\ndef animate_blue_line(num,axi=ax):\n    blue_line.set_data([num, num], [+100000, -100000])\n    return axi,\n\n# Fourier Transformation to show what's \"happening\" at that time\ndef init_rfft(axi=ax):\n    global rfft_line\n    fourier = rfft(samp.v.loc[min_index:min_index+10000].values)[1:1001]\n    fourier = abs(fourier)\n    axi.set_xlim(0, 1000)\n    axi.set_ylim(-500_000, 3_000_000)\n    rfft_line, = axi.plot([],[], color=\"blue\")\n\ndef animate_rfft(num, axi=ax):\n    fourier = rfft(samp.v.loc[num:num+10000].values)[1:1001]\n    fourier = abs(fourier)\n    rfft_line.set_data(np.arange(0,len(fourier)), fourier)\n\ndef init_all():\n    init_rfft(ax[0])\n    init_blue_line(ax[1])\n\n# this is just how matplotlib animation works - we need one function to update the whole frame\ndef animate_all(num):\n    animate_rfft(num, ax[0])\n    animate_blue_line(num, ax[1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c876596333ee15ccbe49daa0dc44b739d7ca60d5","_kg_hide-input":true},"cell_type":"code","source":"# imports and setting for animation\nfrom matplotlib import animation, rc\nrc('animation', html='jshtml')\n\n# defining frames\nframes = np.arange(samp.index.min(), samp.index.max(), AUDIO_RATE/VIDEO_RATE)#[:24]\n\n# creating animation\nanim = animation.FuncAnimation(fig, animate_all, init_func=init_all, frames=frames, interval=1000/VIDEO_RATE, blit=False)\nanim","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70f089e6b303f425d36d3bbe00ffadb72f34b7e6"},"cell_type":"markdown","source":"## Putting two together\nNow all that's left is to combine sound with video. This is not a particularly fast or easy process in a Kaggle notebook; be sure you are ready to wait if you want to re-run it. Also, some libraries are missing from Kaggle image - so the next cell does a little \"fix\"."},{"metadata":{"trusted":true,"_uuid":"02c19d3d2b4432959d555fa1868b935d6380601f","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# installing video libraries\n# this might take several minutes\n!mkdir -p $HOME/.imageio/ffmpeg/\n!cp ../input/skrillex-codecs/ffmpeg-linux64-v3.3.1   $HOME/.imageio/ffmpeg/ffmpeg-linux64-v3.3.1\n!cp ../input/skrillex-codecs/moviepy-0.2.3.5.tar.gz_ ./moviepy-0.2.3.5.tar.gz\n!pip install ./moviepy-0.2.3.5.tar.gz\n!rm ./moviepy-0.2.3.5.tar.gz","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0c7e6ae06e87c270328d1c8ae6a566743d4a5d37","_kg_hide-input":true},"cell_type":"code","source":"# importing movie library\nimport moviepy.editor as mpe\n\n# creating animation with sound - this takes literally forever! (expect 5 mins or more)\nfig, ax = plt.subplots(1, 2,figsize=(18,6))\nanim = animation.FuncAnimation(fig, animate_all, init_func=init_all, frames=frames, interval=1000/VIDEO_RATE, blit=False)\nanim.save('animation.gif', fps=VIDEO_RATE)\nmy_clip = mpe.VideoFileClip('./animation.gif')\naudio_background = mpe.AudioFileClip('./sound.wav')\nfinal_audio = mpe.CompositeAudioClip([audio_background])\n# this looks like a hack and it is.. it should be that \"my_clip.fps==VIDEO_RATE\"\n# but for some reason it's not. Matplotlib isn't the best video saver..\n# or moviepy isn't the best movie reader.. it's one or the other.\n# basically when we're writing, we're setting fps to 24..\n# but when we're reading, it's 25 fps all of a sudden\n# Either this \"anim.save('animation.gif', fps=VIDEO)\" doesn't respect the FPS\n# or this \"mpe.VideoFileClip('./animation.gif')\" doesn't read the file properly\n# doesn't relly matter who... doing \"speedx\" does a \"quick fix!\"\nfinal_clip = my_clip.speedx(VIDEO_RATE/my_clip.fps).set_audio(final_audio)\nfinal_clip.write_videofile(\"./animation_with_sound.mp4\")\n\n# finally display produced animation\nHTML(\"\"\"\n<div>\n<video width=\"760\" height=\"440\" controls>\n  <source src=\"./animation_with_sound.mp4\" type=\"video/mp4\">\n</video>\n</div>\"\"\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a4435ab6d430a1b10283ed0a8dd6dd32925c43ac"},"cell_type":"markdown","source":"## Compare Laboratory Earthquake to Dubstep Drop\nTo make sure we did everything correctly, let's do the exact same procedure for a piece of music - \"Nero & Skrillex - Promises\". The used sample is under 15 seconds, so I hope \"fair use\" applies!"},{"metadata":{"trusted":true,"_uuid":"cb05224da6bb7b416b0f51c159dc6a52ae1a4bcd","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"x = mpe.AudioFileClip('../input/skrillex-codecs/promises_short.wav').to_soundarray()[:,0]\n\nsamp = pd.DataFrame()\nsamp[\"s\"] = x\nsamp[\"v\"] = x*7_000\n\nDROP = 8.314 * AUDIO_RATE\nquake_indices = [DROP]\nmin_index = samp.index.min()\nfig, ax = plt.subplots(1, 2,figsize=(18,6))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cb05224da6bb7b416b0f51c159dc6a52ae1a4bcd","_kg_hide-input":true},"cell_type":"code","source":"frames = np.arange(samp.index.min(), samp.index.max(), AUDIO_RATE/VIDEO_RATE)\n\nanim = animation.FuncAnimation(fig, animate_all, init_func=init_all, frames=frames, interval=1000/VIDEO_RATE, blit=False)\nanim.save('animation_drop.gif')\n\nmy_clip = mpe.VideoFileClip('./animation_drop.gif')\naudio_background = mpe.AudioFileClip('../input/skrillex-codecs/promises_short.wav')\nfinal_audio = mpe.CompositeAudioClip([audio_background])\nfinal_clip = my_clip.speedx(VIDEO_RATE/my_clip.fps).set_audio(final_audio)\nfinal_clip.write_videofile(\"./animation_with_dubstep.mp4\")\n\nHTML(\"\"\"\n<div>\n<video width=\"760\" height=\"440\" controls>\n  <source src=\"./animation_with_dubstep.mp4\" type=\"video/mp4\">\n</video>\n</div>\"\"\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4dd4fe0116db4e077bf6821a1d6ad83472387268"},"cell_type":"markdown","source":"*Thank you for your attention - I know that this Notebook is a bit of a disorganised dump of ideas, but I really just wanted to put something out already! Any feedback greatly appreciated, and perhaps I could improve on it in the future!*"}],"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}