{"cells":[{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport os\nimport seaborn as sns\nprint(os.listdir(\"../input\"))\nPATH = '../input/'\n# Any results you write to the current directory are saved as output.\nnp.random.seed(512)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true},"cell_type":"markdown","source":"# Earthquake prediction \n\nWell my first kernel, and just some work in progress (many a fast evening hacks ...)\n\nIdea collection and some stupid(?) modeling ideas.\n\nI had a look before at  https://www.kaggle.com/jsaguiar/seismic-data-exploration by aguiar, who points out many things, that could be an issue for this kernel... And I will sometimes refer to the quite nice analysis in the kernel. \n\n\nFirst things first:\n1. We got one big training file\n    * It's possible a concatenation of different experiments (see other kernels, different time steps?\n2. We got many small test sets.\n    * Of 150.000 data points of acoustic data\n\nLet's have a look at it.\nTo lower the memory burden a little bit, I am only reading every second row for now. Not sure about the effect on the data yet. (Doesn't work on the kernel, higher memory load, than without...)\n\n## TODO:\nGet that timing clear!"},{"metadata":{"_uuid":"69dde763bfaeebfe48fe889a9a38b47228b9cae0","trusted":true},"cell_type":"code","source":"# https://nikgrozev.com/2015/06/16/fast-and-simple-sampling-in-pandas-when-loading-data-from-files/ \n# Following the approach to index the rows in pandas skiprows. Trust me with the lenght or run:\n# train_lines = sum(1 for l in open(f'{PATH}train.csv'))\ntrain_lines = 629145481 - 1 # there needs to be some substraction (0 indexing)\n# if not using pandas with skiprows:\n#skridx = np.arange(0, train_lines, 2)\nsamp_length = 75000","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4f35ceffc4322fe8f065d98e025f24a8a9f6e20c","trusted":true},"cell_type":"code","source":"# I used the dtype options from aguiars kernel, as this apparently solved the memory overflow in kaggle kernel\ntrain = pd.read_csv(f'{PATH}train.csv', dtype={'acoustic_data': np.int16, \n                                               'time_to_failure': np.float64}) # something less than 5gb, actually - things were getting worse with skiprows, so now ... \ntest_ids = pd.read_csv(f'{PATH}sample_submission.csv') \n\n# Well decimating now here:\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"16276e51d5b7ba5a951f1204175dd6ba1eac772b","trusted":true},"cell_type":"code","source":"# skridx = np.arange(0, train_lines, 2) # Decimate and hoping not to kill the memory again. \ntrain = train.iloc[::2, :]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2b69d9c1a2b02d365eed5bfa54baa6f889c9e2c5","trusted":true},"cell_type":"code","source":"# Short overview over the file format, note everything is decimated a bit:\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6750bbcb5dfb06f3275f3651934d314b8673583d","trusted":true},"cell_type":"code","source":"# Many many samples\ntrain.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8265d685db190de7d027f1002467a618f4fbd2c2","trusted":true},"cell_type":"code","source":"# We need to predict the time_to failure for simple snippets.\ntest_ids.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68aceaee5d42e9ff07eb352c1ae3198661349b33","trusted":true},"cell_type":"code","source":"# 2624 test files.\ntest_ids.shape ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5b8d1dadd3563108cca6b50a7bb522332eab880a"},"cell_type":"markdown","source":"Immediate question before going on with the training data: What is the length of each test set?\n``test_length = [pd.read_csv(f'{PATH}test/{i}.csv').shape[0] for i in test_ids['seg_id']]``\n\nWell ... it's 150000, which can be found in many of the other kernels (i.e. aguiars). \nRunning the line above would be a waste of time :P "},{"metadata":{"_uuid":"43a703c1b5130fcfe9c11f35db0b5eaf291f5dd0"},"cell_type":"markdown","source":"# Data overview"},{"metadata":{"_uuid":"b5938e124e82d119b8a1601f09e33e06dce09143","trusted":true},"cell_type":"code","source":"train['time_to_failure'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c82d2660ba72f7f548711137b01123982b8754a6","trusted":true},"cell_type":"code","source":"train['acoustic_data'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"901d635603762504571c27abce61feb1c8e8ad25"},"cell_type":"markdown","source":" Next we have a look at another, question which is interesting for me: What is the sampling frequency. \n How many data points are in the acoustic data before a change in time_to_failure occurs?"},{"metadata":{"_uuid":"f374d3d6cacab9df5a977127d14a5f6ed0d99dc9","trusted":true},"cell_type":"code","source":"unique_sampling_steps = np.unique(np.diff(np.round(train.time_to_failure, 9)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f374d3d6cacab9df5a977127d14a5f6ed0d99dc9","trusted":true},"cell_type":"code","source":"# We know there are several events. So a slightly different plot:\nplt.figure(figsize=(10,5))\nplt.subplot(121)\nplt.hist(unique_sampling_steps)\nplt.title('Steps with events')\nplt.subplot(122)\nplt.hist(unique_sampling_steps[unique_sampling_steps < 2])\nplt.title('Removing the long events.')\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a940891f9b3757cf28f858e7531186d738c8812d"},"cell_type":"markdown","source":"But apparently, there are different  sampling periods (see the referred kernel...).  \nWe can also see this above.\n\nSo maybe later: Discretize the data a bit more (well here already half the data was skipped anyways...)\n\nStupid other idea: Strict discretization and clustering? "},{"metadata":{"_uuid":"e05bcb473989bfcf2354a9138ea86ef16f5d0346","trusted":true},"cell_type":"code","source":"print(np.where(np.diff(np.round(train.time_to_failure[:100000], 5)))) # I.e. decreasing precision ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3a1352a9c109c0fc9c200c611abbb3360f9c8d50","trusted":true},"cell_type":"code","source":"# Now lets have a look at the audio signal around failures, with the index we are not really caring about being one off..\n# looking for increasing times as there's a restart in the signal\nfailures = np.where(np.diff(train.time_to_failure) > 0)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5d8d211f29fe64e564bf07be662d4c858d9b5497","trusted":true},"cell_type":"code","source":"# Seems like there are 16 failures:\nprint(failures[0].shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ceace4ab02ca5a23b403d2717a25d4668fd94a8f"},"cell_type":"markdown","source":"## Looking at the 150000 samples before each event:"},{"metadata":{"_uuid":"ff8fa9c2aa3b7770ce8a82fb4a4539a184c0e938","trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\nfor ii, f_idx in zip(range(failures[0].shape[0]), failures[0]):\n    plt.subplot(4, 4, ii + 1)\n    plt.plot(np.arange(-samp_length, 100), train.acoustic_data[f_idx - samp_length : f_idx + 100])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"98267e82692bd22bc3f8acb954015d554bf11e44"},"cell_type":"markdown","source":"## And look at some more data before events:\nWell this is not very informative /  representative.  "},{"metadata":{"_uuid":"71ebee5849ee9546d899b327c60a5ee46b6867d5","trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,10))\nfor ii, f_idx in zip(range(failures[0].shape[0]), failures[0]):\n    plt.subplot(4, 4, ii + 1)\n    plt.plot(np.arange(-1200000, 0), train.acoustic_data[f_idx - 1200000 : f_idx])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2e6b9a9c638739da8b942bc28fba62a5f44a77fb"},"cell_type":"markdown","source":"Kinda smooth plots before each event, not high value data. Just by eyeballing it seems like that there is a decay in high values, the closer we are to an earthquake. But well, this is not the best way to look at the data. \n\nAdditional point: We are not able to rely on too many data points before events to make our predictions. Rather (as the goal of the challenge is) we have to create a model to just predict the time before an event. \nCutting into samples could solve this maybe.\n\nWe also see that there are huge spikes in the data!"},{"metadata":{"_uuid":"550dc33c96e53c6484ae7cd8f60c108dd46d8fd3"},"cell_type":"markdown","source":"# Learn a bit more about time to failure"},{"metadata":{"_uuid":"34201298ddbbf4cb524eb47797f416f5bbe6db5a"},"cell_type":"markdown","source":"Our goal is to predict one value for time_to_failure for each training segment. A question that came to my mind, also due to the different sampling rates: How much does time_to_failure vary in random 150.000 (or 75.000) sample pieces. "},{"metadata":{"_uuid":"3632737debf3239a9875cb661480c252d6f718df","trusted":true},"cell_type":"code","source":"n_rnd_samples = 1000\nrnd_sample_idx = np.random.randint(low=samp_length, high=train.shape[0], size=n_rnd_samples)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"26652b28a3a448d11a6d0d71829c9f5f2d8ca487","trusted":true},"cell_type":"code","source":"variation = np.empty(n_rnd_samples)\n\nfor idx, samp in enumerate(rnd_sample_idx):\n    variation[idx] = np.max(train.time_to_failure[samp - samp_length : samp]) - np.min(train.time_to_failure[samp - samp_length : samp])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0103662a1b248640353f8e3d1270de81ac398c5c"},"cell_type":"markdown","source":"## Actually - not so much"},{"metadata":{"_uuid":"abf1872505c0092f6b9fdda6fa164ce7647cf730","trusted":true},"cell_type":"code","source":"b = plt.hist(variation[variation < 4])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b0516288e7322cf42c781de6f283f2709a84ca7d"},"cell_type":"markdown","source":"# Learn a bit more about acoustic data\n\n## Spikes\n\nWe have seen some huge spikes in the data. Are these common?"},{"metadata":{"_uuid":"b8bc31f498ca04a7235c326aac86346dc6b15921","trusted":true},"cell_type":"code","source":"b = plt.hist(np.clip(train.acoustic_data, -200, 200), bins=50)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2b7294ca5cd8da2e5f6022492453cc7e3e88dea3"},"cell_type":"markdown","source":"Apparently not"},{"metadata":{"_uuid":"960b5e0799a45ca8aecdb3a7de09e083455bf812"},"cell_type":"markdown","source":"## Can Spikes tell us something about time_to_failure?\nOnly using absolute values for spikes"},{"metadata":{"_uuid":"a49b42847644f7c5f2a27328234680f1b82246a9","trusted":true},"cell_type":"code","source":"from scipy.stats import spearmanr\nspike_size = 800\nspike_idx = np.abs(train.acoustic_data) > spike_size\n\nplt.figure(figsize=(15,5))\nspear = spearmanr(np.abs(train[spike_idx]['acoustic_data']), train[spike_idx]['time_to_failure'])\nplt.subplot(121)\nplt.scatter(np.abs(train[spike_idx]['acoustic_data']), train[spike_idx]['time_to_failure'])\nb = plt.title(f'Peaks and time to failure correlate with r ={spear[0]:.3f}')\nplt.subplot(122)\nb = plt.hist(train[spike_idx]['time_to_failure'])\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cf60869daf90659d9f0d619a49a8711c611d226d"},"cell_type":"markdown","source":"We can see that there seems to be some relationship between spikes and time to failure.\nHuge spikes seem to occur somewhat more readily before failure events!\n\nBut also caution: Some spikes occur long before!"},{"metadata":{"_uuid":"03df283ada4fe702774d2d5b0584db9e886827ac"},"cell_type":"markdown","source":"## Do some random sampling again: Look at the max-values in an epoch and time_to_failure"},{"metadata":{"_uuid":"28402dc011c06077d2f2bb0760f3c5727b60bccc","trusted":true},"cell_type":"code","source":"t_t_f = np.empty(n_rnd_samples)\nm_a_d = np.empty(n_rnd_samples)\n\nfor idx, samp in enumerate(rnd_sample_idx):\n    t_t_f[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    m_a_d[idx] = np.max(np.abs(train.acoustic_data[samp - samp_length : samp]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5058fcc52736a98d62b6f92cffcce599772ef934","trusted":true},"cell_type":"code","source":"spear = spearmanr(np.array(t_t_f), np.array(m_a_d))\n\nplt.scatter(np.array(t_t_f), np.array(m_a_d))\nb = plt.title(f'Peaks and time to failure correlate with r = {spear[0]:.3f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cc7f9a75dfd7500f984f3f90a6cc3e835176ce67"},"cell_type":"markdown","source":"We can actually see, that some really basic features (such as the max-value) has some predictive ability!"},{"metadata":{"_uuid":"b778a15e230960044941da4957d456ee3acc98fd"},"cell_type":"markdown","source":"# Extend the notion - can we use basic properties of the distribution?"},{"metadata":{"_uuid":"5e11683512b4af6192c141fd0febd7c08cad78e8","trusted":true},"cell_type":"code","source":"def basic_properties(data):\n    properties = np.zeros(7)\n    # For many of the basic properties we use the abs, that is the amplitude\n    properties[:3] = list(np.percentile(np.abs(data), [25, 50, 75]))\n    \n    for n_f, jj in enumerate([np.mean, np.std, np.max, np.min]):\n        properties[3 + n_f] = jj(data)\n    return properties\n    \nt_t_f = np.empty(n_rnd_samples)\nm_a_d = np.empty((n_rnd_samples, 7))\n\nfor idx, samp in enumerate(rnd_sample_idx):\n    t_t_f[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    m_a_d[idx, :] = basic_properties(np.abs(train.acoustic_data[samp - samp_length : samp]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c54beccdac695e9e3305fb2ec8ce547e137d8901","trusted":true},"cell_type":"code","source":"plt.figure(figsize=(25,20))\nfeat_names = ['25-percentile', '50-percentile', '75-percentile', 'mean', 'std', 'max', 'min']\nfor ii, fna in enumerate(feat_names):\n    ax = plt.subplot(3, 4, ii + 1)\n    sns.regplot(t_t_f, m_a_d[:, ii], ax=ax)\n    spear = spearmanr(np.array(t_t_f), np.array(m_a_d[:, ii]))\n    b = plt.title(f'{fna} and ttf correlate with r = {spear[0]:.3f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c78ba87153a0bf1fc8d82941e68472008008dd09"},"cell_type":"markdown","source":"Another value I just remembered that is often used in the auditory domanin ist the RSM (root-squared-mean). Which should correlated highly with values in the above. But I'll but it in here as well. "},{"metadata":{"trusted":true,"_uuid":"6a9f7386fe0eb8c49e6009aceecd0f383b8f3d46"},"cell_type":"code","source":"t_t_f = np.empty(n_rnd_samples)\nm_a_d = np.hstack([m_a_d, np.zeros((m_a_d.shape[0],1))])\n\nfor idx, samp in enumerate(rnd_sample_idx):\n    t_t_f[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    m_a_d[idx, 7] = np.sqrt(np.mean(train.acoustic_data[samp - samp_length : samp] ** 2))\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7ebb0065a1cde36d72a7e290ec986ddb2601dc65"},"cell_type":"code","source":"plt.figure(figsize=(15,5))\nax = plt.subplot(1, 2, 1)\nsns.regplot(t_t_f, m_a_d[:, 7], ax=ax)\nspear = spearmanr(np.array(t_t_f), np.array(m_a_d[:, 7]))\nb = plt.title(f'RMS and ttf correlate with r = {spear[0]:.3f}')\nax = plt.subplot(1, 2, 2)\nsns.regplot(m_a_d[:,3], m_a_d[:, 7], ax=ax)\nspear = spearmanr(m_a_d[:,3], np.array(m_a_d[:, 7]))\nb = plt.title(f'RMS and abs(mean) correlate with r = {spear[0]:.3f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35c4edaba46f59214d9b209ee41760564eff565a"},"cell_type":"markdown","source":"However, it might be interesting to check whether normalizing the data to have the same RMS will have an effect on predictive ability (generalization etc.)"},{"metadata":{"_uuid":"6b5c13120775cde40850fe8afe57f6b5066f054f"},"cell_type":"markdown","source":"## Punchline:\nSome basic features could actually do something already! We tried this already, see below Version 7, (LB=1.962)"},{"metadata":{"_uuid":"e4c863c910b8df7172f26f80c6b46cb31c2954c4"},"cell_type":"markdown","source":"# Addendum 1\nUsing amplitude, power and other simple features deal with the data at hand directly. Another important feature of the audio signal could be the change over time. That is calculating the derivative, but now using the simple difference over time.  \n**Important now**: The decimation step in the beginning is now affecting the data!"},{"metadata":{"trusted":true,"_uuid":"d3328b8ee9b092354ab202eda766809dbbcd6e74"},"cell_type":"code","source":"# This might use up a lot of memory... \ntrain['acoustic_diff'] = np.hstack([0, np.diff(train['acoustic_data'])])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a3385b0c11ade1edefd1ab581e9548596cc42747"},"cell_type":"code","source":"# simple description:\ntrain['acoustic_diff'].describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f9eb94731d64da925f0e94bcbdee07e246c14541"},"cell_type":"markdown","source":"## Side note\nI am not plotting the time series before the events, as those look pretty similar to the normal data. So there is not too much information to be gained by this. "},{"metadata":{"_uuid":"6bd4a637857942b1afe175974c077f0167c4fdcb"},"cell_type":"markdown","source":"## Redoing the basic features analysis:"},{"metadata":{"trusted":true,"_uuid":"16314930ff822d23c3e02399e34dc7a661dc6591"},"cell_type":"code","source":"\ndef RMS(data):\n    return np.sqrt(np.mean(data**2))\n\ndef basic_properties_up(data):\n    properties = np.zeros(8)\n    # For many of the basic properties we use the abs, that is the amplitude\n    properties[:3] = np.percentile(data, [25, 50, 75])\n    \n    for n_f, jj in enumerate([np.mean, np.std, np.max, np.min, RMS]):\n        properties[3 + n_f] = jj(data)\n    return properties\n    \nt_t_f = np.empty(n_rnd_samples)\nm_a_d = np.empty((n_rnd_samples, 8))\n\nfor idx, samp in enumerate(rnd_sample_idx):\n    t_t_f[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    m_a_d[idx, :] = basic_properties_up(train.acoustic_diff[samp - samp_length : samp])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8de0ab9fca3fb05c07c919b7266593ca921b17f8"},"cell_type":"code","source":"plt.figure(figsize=(25,20))\nfeat_names = ['25-percentile', '50-percentile', '75-percentile', 'mean', 'std', 'max', 'min', 'RMS']\nfor ii, fna in enumerate(feat_names):\n    ax = plt.subplot(3, 4, ii + 1)\n    sns.regplot(t_t_f, m_a_d[:, ii], ax=ax)\n    spear = spearmanr(np.array(t_t_f), np.array(m_a_d[:, ii]))\n    b = plt.title(f'{fna} and ttf correlate with r = {spear[0]:.3f}')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"38f8f17246fa77079db8a39082d2be0b05a6c6a0"},"cell_type":"markdown","source":"Some basic properties of change also correlate with time to failure. But without further analysis we cannot be sure, how unique those properties are! "},{"metadata":{"_uuid":"27b6a5cf2a1fc2ad8bb89c0eabc6f0aef73d3c0a"},"cell_type":"markdown","source":"# Addendum 2: time-frequency features"},{"metadata":{"_uuid":"8a92b25a3c3681355bf3d5913ea9ea8d7b596a5e"},"cell_type":"markdown","source":"Some common features that are used in the auditory domain are time frequency estimates. I am pretty much agnosticly applying those here. I have some ideas how to implement them more nicely, especially for visualization, but again, some evening hack... "},{"metadata":{"trusted":true,"_uuid":"20b046d48c8ccdedbb8cc32359c3b9f1a71f8f55"},"cell_type":"code","source":"import librosa # Seems to be a common library for all thing auditory","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"42548f36bff4f16b723814819c291ae84c8298a4"},"cell_type":"markdown","source":"A feature quite often used are mel frequencies. Do they work here? I don't know. Very often used in speech processing!\nWorking with defaults here. But let's see what happens in another domain, right before an event.\nBut we need a sampling rate here... "},{"metadata":{"trusted":true,"_uuid":"45aa6c391399a738b6b9bbfd2f6d80734009d4ba"},"cell_type":"code","source":"# Assuming that t_t_f is in seconds, yes, I'm lazy, still didn't look that up, shame on me:\nsr = round(1 / (train.iloc[0]['time_to_failure'] - train.iloc[1]['time_to_failure']))\nprint(sr)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ecc8053134f3f57dce7f101a15d3fa5df4282227"},"cell_type":"markdown","source":"Still not having much of an idea about earthquake data, that seems to be a quite high sr. "},{"metadata":{"trusted":true,"_uuid":"451bd44f46da2f665836c83463a4218f37878322"},"cell_type":"code","source":"(train.iloc[0]['time_to_failure'] - train.iloc[75000]['time_to_failure'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e245d5940a470c406738bb0ed9e05be664f5cac1"},"cell_type":"markdown","source":"And something is wrong. Let's just use the defaults and see how far we can get"},{"metadata":{"trusted":true,"_uuid":"72dc2b63d4978a901b348de1ead969502ef6ce8d"},"cell_type":"code","source":"def mfcc_wrap(data):\n    return librosa.power_to_db(librosa.feature.melspectrogram(data.values.astype('float32'), n_mels=25))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8394a620ec11d1b5bb49858fe71cc9d02ede3963"},"cell_type":"markdown","source":"# Hm... some look at db in a 75000 sample"},{"metadata":{"trusted":true,"_uuid":"f9ce3e7248fa2a27de7a3a066475acf423c35f40"},"cell_type":"code","source":"plt.figure(figsize=(15,10))\nfor ii, f_idx in zip(range(failures[0].shape[0]), failures[0]):\n    ax = plt.subplot(4, 4, ii + 1)\n    sns.heatmap(mfcc_wrap(train.acoustic_data[f_idx - samp_length : f_idx + 100]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5df484c549ab2d2dbe2dfb572edef3be5ec73481"},"cell_type":"code","source":"mfcc_wrap(train.acoustic_data[ : samp_length]).ravel().shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5dac4ef6d7ab4b9410e0b45b3b4250fb74953641"},"cell_type":"markdown","source":"Not really sure what to make off it... yet... this will happen another time..."},{"metadata":{"_uuid":"cb28ef3137be23a9e7404a5a65a98fd79e5cfba0"},"cell_type":"markdown","source":"# Let's try some really simple prediciton\nIn the prior version there was stuff done using fft features (which was probably quite wrong on many levels), it didn't do much too say the least^^\n\nBy looking at the other analysis, it was probably a relation between max values and time_to_failure"},{"metadata":{"_uuid":"cd95c5db838e531aed7a47f4e4aa5dc2c39be7a0","trusted":true},"cell_type":"code","source":"from sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_absolute_error","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c3daddaf99f729568fbdac8665235d0d3c718993","trusted":true},"cell_type":"code","source":"max_features = np.round(train.shape[0]/ samp_length).astype(int)\nprint(f'We can create {max_features} unique samples. Maybe this gets us somewhere')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b74a12dd472e32b087d69a751000501f66fd0e74","trusted":true},"cell_type":"code","source":"def basic_properties_time(data):\n    properties = np.zeros(7)\n    # For many of the basic properties we use the abs, that is the amplitude\n    properties[:3] = np.percentile(data, [25, 50, 75])\n    \n    for n_f, jj in enumerate([np.mean, np.std, np.max, RMS]):\n        properties[3 + n_f] = jj(data)\n    return properties","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8796525ed3b2a916982e96ba14c41f60be8e70a0"},"cell_type":"code","source":"def basic_properties_diff(data):\n    properties = np.zeros(6)\n    properties[:2] = np.percentile(data, [25, 75])\n    \n    for n_f, jj in enumerate([np.min, np.std, np.max, RMS]):\n        properties[2 + n_f] = jj(data)\n    return properties","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0d2cf4724c7df8cad69245ba4a02037b9c123c6b"},"cell_type":"markdown","source":"## Some things that might help\nWe can create 4194 unique samples. But I think random sampling from the data to create our test and validation sets might be the way to got. Let's see how far we can get! "},{"metadata":{"trusted":true,"_uuid":"c2ab36d9a342069032d43676c15821222617e8cd"},"cell_type":"code","source":"train =  train.drop(['acoustic_diff'], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25eb5f85808d45560b0f3cfe4a07824642acfe3c","trusted":true},"cell_type":"code","source":"n_samples = 15000\nX_time = np.zeros((n_samples, 7))\nX_diff = np.zeros((n_samples, 6))\nX_mfcc = np.zeros((n_samples, 3675))\ny = np.zeros((n_samples))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cd57150367b53d5b1aa603dbac4a55687790a3f0","trusted":true},"cell_type":"code","source":"sample_idx = np.random.randint(low=samp_length, high=train.shape[0], size=n_samples)\n\nfor idx, samp in enumerate(sample_idx):\n    y[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    X_time[idx, :] = basic_properties_time(np.abs(train.acoustic_data[samp - samp_length : samp]))\n    X_diff[idx, :] = basic_properties_diff(np.diff(train.acoustic_data[samp - samp_length : samp]))\n    X_mfcc[idx, :] = mfcc_wrap(train.acoustic_data[samp - samp_length : samp]).ravel()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7032d0572f03d6e8ee57d29a92dc8397d8f0ce78","trusted":true},"cell_type":"code","source":"pred_mean = np.ones(y.shape) * np.mean(X_time[:, 3])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fbc05d0474f0a7e4912f304965c7ce5afd617015"},"cell_type":"markdown","source":"# Predictions:"},{"metadata":{"_uuid":"4fc19bd05b75a615264b61d1a0da600497904d8a"},"cell_type":"markdown","source":"## Base line - what we definitely want to beat!"},{"metadata":{"_uuid":"ef0bba09342c30fc5d57ab42bae73bf2a2fb8fab","trusted":true},"cell_type":"code","source":"print(mean_absolute_error(y, pred_mean))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8c89e849ca0190a996e2973291fcde6d27446e07"},"cell_type":"markdown","source":"## Trees have been promising before"},{"metadata":{"_uuid":"ea597888921d31a671caeade5040a4e6348461c5","trusted":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d19b2267658554a9890c5283d4d00e1f538cc7a2","trusted":true},"cell_type":"code","source":"# Make pred_tree greater than 0\npred_time = cross_val_predict(RandomForestRegressor(n_estimators=100), X_time, y, cv=3)\npred_time[pred_time < 0] = 0","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"83253b4a8b1b4f115177f0691549c48edd0c5a1e","trusted":true},"cell_type":"code","source":"print(mean_absolute_error(y, pred_time))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"425768296f43b1ffc6dcf31c729959069f8bb92c"},"cell_type":"markdown","source":"At least better than baseline... "},{"metadata":{"_uuid":"943d8ed7e958032354e5dea2997dbbd307c3abb9"},"cell_type":"markdown","source":"## Time derivative features:"},{"metadata":{"trusted":true,"_uuid":"5c97389adafb4585c7d76cb6f33e3ff217dfc428"},"cell_type":"code","source":"pred_diff = cross_val_predict(RandomForestRegressor(n_estimators=100), X_diff, y, cv=3)\npred_diff[pred_diff < 0] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"22ff745f2365c7c2e236a6e6ad935fa39343ffd9"},"cell_type":"code","source":"print(mean_absolute_error(y, pred_diff))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b9b55a6a6513c77a84da234d99df72099cc3f39e"},"cell_type":"markdown","source":"A bit better still ..."},{"metadata":{"trusted":true,"_uuid":"c4ad494a3310ae787c2810525bf54fcb9c2304c5"},"cell_type":"code","source":"# simple average:\nprint(mean_absolute_error(y, (pred_time + pred_diff)/2))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"306a6842405d4a6276610823c7146f254ba09194"},"cell_type":"markdown","source":"## Time and derivatives combined"},{"metadata":{"trusted":true,"_uuid":"9e960105d3dce076fd5be34aa127fa2c04f317af"},"cell_type":"code","source":"pred_comb = cross_val_predict(RandomForestRegressor(n_estimators=100), np.hstack([X_time, X_diff]), y, cv=3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3288f0984bcf693405e4bf34a5a0c739ff72efcb"},"cell_type":"code","source":"pred_comb[pred_comb < 0] = 0\nprint(mean_absolute_error(y, pred_comb))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"32ff1256c3fc0c3de53f982ef35b6e63cfddcb98"},"cell_type":"markdown","source":"#### Bottom line:\nTime derivative features seem to contain at least some more information"},{"metadata":{"_uuid":"fb644951f33c32d872007c892cb3726e3a60be7c"},"cell_type":"markdown","source":"## A quick look a the mel frequencies"},{"metadata":{"trusted":true,"_uuid":"01b7c9e1b0bd12bcb39b2a84b5721fba39013308"},"cell_type":"code","source":"pred_mel = cross_val_predict(RandomForestRegressor(n_estimators=10, n_jobs=3), X_mfcc, y, cv=3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c9d16dd5e3e51b8f8ac79faebd2121f779c610fc"},"cell_type":"code","source":"pred_mel[pred_mel < 0] = 0\nprint(mean_absolute_error(y, pred_mel))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f216c51c8a8be9a26683e68594d69c5d51ab3df"},"cell_type":"markdown","source":"Something more has to be done here! ... Probably understanding wise... "},{"metadata":{"_uuid":"38c59ff59ef991e665a34b25e21dba0602d03cec"},"cell_type":"markdown","source":"## And the average"},{"metadata":{"trusted":true,"_uuid":"ef5de2361a1be79dbcc305bb5792e240d87d28b5"},"cell_type":"code","source":"print(mean_absolute_error(y, (pred_mel + pred_comb)/2))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0cd2b2d5af331c4db88c31e842adbcca4ec99baf"},"cell_type":"markdown","source":"# A last test - does normalization of batches help?"},{"metadata":{"trusted":true,"_uuid":"03a0de73cbeafab433d931c3782a5f78b9f05763"},"cell_type":"code","source":"data_rms = RMS(train.acoustic_data)\nprint(data_rms)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"613a4708b17b471291c17812ca0503860d90d771"},"cell_type":"markdown","source":"## We are now recreating the data set and look at a first basic preprocessing step"},{"metadata":{"_uuid":"0ecc93b0a7151401393d0a4f636db3e4f31f4838"},"cell_type":"markdown","source":"Basically we are rescaling every batch in the training set by the factor data_rms / set_rms"},{"metadata":{"trusted":true,"_uuid":"6576aaad7ea5248df3b58febd616af71b943d5d5"},"cell_type":"code","source":"n_samples = 15000\nX_time = np.zeros((n_samples, 7))\nX_diff = np.zeros((n_samples, 6))\nX_mfcc = np.zeros((n_samples, 3675))\ny = np.zeros((n_samples))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6ab38ac59cad6418e8d766a7d2698ce3c8c1a142"},"cell_type":"code","source":"sample_idx = np.random.randint(low=samp_length, high=train.shape[0], size=n_samples)\n\nfor idx, samp in enumerate(sample_idx):\n    y[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    temp_data = train.acoustic_data[samp - samp_length : samp]\n    temp_rms = RMS(temp_data)\n    temp_data = temp_data * (data_rms/temp_rms)\n    \n    X_time[idx, :] = basic_properties_time(np.abs(temp_data))\n    X_time[idx, -1] = temp_rms # keeping the old rms as a feature basic_properties_time(np.abs(temp_data))\n    X_diff[idx, :] = basic_properties_diff(np.diff(temp_data))\n    X_mfcc[idx, :] = mfcc_wrap(temp_data).ravel()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3a9f79cdcf1bcf040a34881095cd06745133180b"},"cell_type":"markdown","source":"## The rescaled error values:"},{"metadata":{"trusted":true,"_uuid":"47d5f7b84b3c115c48a0828bf6e0d499ad3bb218"},"cell_type":"code","source":"pred_time = cross_val_predict(RandomForestRegressor(n_estimators=100), X_time, y, cv=3)\npred_time[pred_time < 0] = 0\nprint(mean_absolute_error(y, pred_time))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3fc841df12dd36a1e035863ef49352f1382e9626"},"cell_type":"code","source":"pred_diff = cross_val_predict(RandomForestRegressor(n_estimators=100), X_diff, y, cv=3)\npred_diff[pred_diff < 0] = 0\nprint(mean_absolute_error(y, pred_diff))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"721931367e9cd5e29765332a2320ba912ec22671"},"cell_type":"code","source":"pred_comb = cross_val_predict(RandomForestRegressor(n_estimators=100), np.hstack([X_time, X_diff]), y, cv=3)\npred_comb[pred_comb < 0] = 0\nprint(mean_absolute_error(y, pred_comb))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68eb181ba3fb2b58e7936489f5b03b29387d85f0"},"cell_type":"markdown","source":"### Bottom line: Rescaling didn't help us "},{"metadata":{"_uuid":"d2827e71af26a8243fd242a2f870685d72f6671e"},"cell_type":"markdown","source":"# At last: Creating a simple submission"},{"metadata":{"trusted":true,"_uuid":"5d14638bbd0dfd47706ca5f0ea603550426031d2"},"cell_type":"code","source":"n_samples = 15000\nX_time = np.zeros((n_samples, 7))\nX_diff = np.zeros((n_samples, 6))\nX_mfcc = np.zeros((n_samples, 3675))\ny = np.zeros((n_samples))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d6a6b5940ff9cd2a63b2570e2484460dc2750679"},"cell_type":"code","source":"sample_idx = np.random.randint(low=samp_length, high=train.shape[0], size=n_samples)\n\nfor idx, samp in enumerate(sample_idx):\n    y[idx] = np.median(train.time_to_failure[samp - samp_length : samp])\n    X_time[idx, :] = basic_properties_time(np.abs(train.acoustic_data[samp - samp_length : samp]))\n    X_diff[idx, :] = basic_properties_diff(np.diff(train.acoustic_data[samp - samp_length : samp]))\n    X_mfcc[idx, :] = mfcc_wrap(train.acoustic_data[samp - samp_length : samp]).ravel()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"940882c2c1a51b90450b9adc689e60dc8faa2313"},"cell_type":"markdown","source":"# And off to the leader board:"},{"metadata":{"_uuid":"1da00a6eae8602a180398f394518f665abb9cf52","trusted":true},"cell_type":"code","source":"model_basic = RandomForestRegressor(n_estimators=100)\nmodel_basic.fit(np.hstack([X_time, X_diff]), y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8330b12dc7f74b10103303f4945b2faee6bff29c"},"cell_type":"markdown","source":"### We are using very basic features, hoping that the distribution is the same, should not mean too much (I hope)"},{"metadata":{"_uuid":"bab429806c6691cab390d955472a621c3aa8560b","trusted":true},"cell_type":"code","source":"submit_df = test_ids.copy()\nfor seg_id in test_ids.seg_id:\n    temp = pd.read_csv(f'{PATH}test/{seg_id}.csv')\n    temp = temp.iloc[::2, :] # we've been using decimated data the whole time...\n    x_pred = np.zeros((1, 13))\n    x_pred[0, :7] = basic_properties_time(np.abs(temp.acoustic_data))\n    x_pred[0, 7:] = basic_properties_diff(np.diff(temp.acoustic_data))\n    X_pred = model_basic.predict(x_pred)\n    \n    submit_df.loc[test_ids.seg_id == seg_id, 'time_to_failure'] = np.max([0, X_pred])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"92327870303a673135879f1fa46a72b9ed6253c6","trusted":true},"cell_type":"code","source":"submit_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","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}