{"cells":[{"metadata":{"_uuid":"93e36bbb6826818b439ae57db1f438168c381646"},"cell_type":"markdown","source":"## What is this competition all about?\n\n* Given seismic signals we are asked to predict the time until the onset of laboratory earthquakes.\n* The training data is a single sequence of signal and seems to come from one experiment alone.\n* In contrast the test data consists of several different sequences, called segments, that may correspond to different experiments. The regular pattern we might find in the train set does not match those of the test segments. \n* For each test data segment with its corresponding seg_id we are asked to predict it's single time until the lab earthquake takes place."},{"metadata":{"_uuid":"d3457484c9d6ba5136f23749b1cd94142df23467"},"cell_type":"markdown","source":"## Loading packages"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":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)\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport seaborn as sns\nsns.set()\n\nfrom IPython.display import HTML\n\nfrom os import listdir\nprint(listdir(\"../input\"))\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=DeprecationWarning)\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", category=FutureWarning)\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d974922e32763d325dcdb936a619ea4733134155"},"cell_type":"markdown","source":"## What is an earthquake in the lab?\n\nCurrently I don't know how an earthquake in the laboratory works and as I like to know I googled around and found this nice video that shows how such a lab looks like. If you like, feel free to take a look at it. I'm still on my journey to understand the problem. "},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"745e8edbf4909f4785f99d237970ec01da376d63"},"cell_type":"code","source":"HTML('<iframe width=\"800\" height=\"400\" src=\"https://www.youtube.com/embed/m_dBwwDJ4uo\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"718d9267ca54042867d0e487b640bb7c6070f3ef"},"cell_type":"markdown","source":"In the end we can see that the probes that are used are put under some kind of **normal pressure but there is a shear stress working on it as well**. Then, after some time, the probe splits. If you take a look at the additional material given, you can see that we have: \n\n### 3 kind of plates\n\n* 2 plates left and right that are under normal pressure: Forces are acting with 90 degree on the plate, pushing the two plates together. \n* In the middle we find a third plate which is separated by some granular material. This plate moves downwards with constant velocity. \n\nI'm not sure if I understand this right, but it seems that this granular material is the \"rock\" that can split and load again to produce this kind of lab earthquakes in repetitive cycles. Even though the train set contains continuous data it contains several such splits (earthquakes)."},{"metadata":{"_uuid":"2eb60d638f110958f6420b79a63dc4083503843f"},"cell_type":"markdown","source":"## Let's get familiar with the data!\n\n### Training data\n\nThe total size of the train data is almost 9 GB and we don't want to wait too long just for a first impression, let's load only some rows: "},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/train.csv\", nrows=10000000,\n                    dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntrain.head(5)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9cd3a4b5a33e3380a21f5c773f8433f871c9599e"},"cell_type":"markdown","source":"We can see two columns: Acoustic data and time_to_failure. The further is the seismic singal and the latter corresponds to the time until the laboratory earthquake takes place. Ok, personally I like to rename the columns as typing \"acoustic\" every time is likely for me to produce errors:"},{"metadata":{"trusted":true,"_uuid":"fd1f72651a03eb7f1cd77df583f668f56d44aedc"},"cell_type":"code","source":"train.rename({\"acoustic_data\": \"signal\", \"time_to_failure\": \"quaketime\"}, axis=\"columns\", inplace=True)\ntrain.head(5)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"763c28eb26d64a9359a8293f27b8718f3069c382"},"cell_type":"markdown","source":"We can see that the quaketime of these first rows seems to be always the same. But is this really true?"},{"metadata":{"trusted":true,"_uuid":"8922228579b246fb69f4a7b311731bcc276acb7b"},"cell_type":"code","source":"for n in range(5):\n    print(train.quaketime.values[n])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"70e6e8faac19376450d1006035b33036995fa5f1"},"cell_type":"markdown","source":"Aha! We can see that they are not the same and that pandas has rounded them off. And we can see that the time seems to decrease. Let's plot the time to get more familiar with this pattern:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"9a593f59bf16cdd00d40f0492e90a004598022fa"},"cell_type":"code","source":"fig, ax = plt.subplots(2,1, figsize=(20,12))\nax[0].plot(train.index.values, train.quaketime.values, c=\"darkred\")\nax[0].set_title(\"Quaketime of 10 Mio rows\")\nax[0].set_xlabel(\"Index\")\nax[0].set_ylabel(\"Quaketime in ms\");\nax[1].plot(train.index.values, train.signal.values, c=\"mediumseagreen\")\nax[1].set_title(\"Signal of 10 Mio rows\")\nax[1].set_xlabel(\"Index\")\nax[1].set_ylabel(\"Acoustic Signal\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3b844108295488f0d636513afe19378727378a59"},"cell_type":"markdown","source":"### Take-Away\n\n* We can see only one time in 10 Mio rows when quaketime goes to 0. This is a timepoint where an earthquake in the lab occurs. \n* There are many small oscillations until a heavy peak of the signal occurs. Then it takes some time with smaller oscillations and the earthquake occurs.\n\n\nIf we take a look at the first 50000 indizes we can see that there is a second pattern of quaketime that may has something to do with the resolution of the experimental equipment:"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"33104fb28f462352eb6da660484088a1d21e24ad"},"cell_type":"code","source":"fig, ax = plt.subplots(3,1,figsize=(20,18))\nax[0].plot(train.index.values[0:50000], train.quaketime.values[0:50000], c=\"Red\")\nax[0].set_xlabel(\"Index\")\nax[0].set_ylabel(\"Time to quake\")\nax[0].set_title(\"How does the second quaketime pattern look like?\")\nax[1].plot(train.index.values[0:49999], np.diff(train.quaketime.values[0:50000]))\nax[1].set_xlabel(\"Index\")\nax[1].set_ylabel(\"Difference between quaketimes\")\nax[1].set_title(\"Are the jumps always the same?\")\nax[2].plot(train.index.values[0:4000], train.quaketime.values[0:4000])\nax[2].set_xlabel(\"Index from 0 to 4000\")\nax[2].set_ylabel(\"Quaketime\")\nax[2].set_title(\"How does the quaketime changes within the first block?\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"efb67f959e7271a458e18e51bd28564e9f7b4736"},"cell_type":"markdown","source":"### Take-Away\n\nVery interesting!\n\n* The first plot shows that the quaketime seems to stay almost constant up to index 4000. Then a steep decrease occurs. Afterwards this kind of pattern is repeated.\n* The second plot reveals that the second jump of the quaketime is larger than the first.\n* The third plot shows that the quaketime within such a \"constant\" block is not really constant but linear decreasing even though with very small numbers."},{"metadata":{"_uuid":"76d33aef4c795ced5d369aa8a1163f2180cbb70f"},"cell_type":"markdown","source":"### First conclusion\n\nThank you to @pete **who pointet out that there are 16 earthquakes in the train set and that they really occur when quaketime goes to 0**. It's still interesting why we have **several blocks and jumps of quaketime on low resolution** where the time decreases linear with a very small stepsize.\n\nCurrently I'm not sure why these jumps between our target quaketime occur between values with only small scaled differences.  **What do you think? Let's me know if you like in the commets ;-) **\n\nOk, I feel more familiar with the train data right now. Let's turn to the test data before switching to more explorations."},{"metadata":{"_uuid":"25a0e765f82730bf5fd9e72476d091bcf22a4944"},"cell_type":"markdown","source":"### Test data\n\nWe can find multiple segments of sequences in the test folder. Let's peek at their names:"},{"metadata":{"trusted":true,"_uuid":"f5f067e8975e93384f3e029c0e2fd11290ab1864"},"cell_type":"code","source":"test_path = \"../input/test/\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bdebb1e826a0b554feec3c05b26d7a29fa86e34a"},"cell_type":"code","source":"test_files = listdir(\"../input/test\")\nprint(test_files[0:5])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f1572948e1b35f70374520339466c3cb5565e7a9"},"cell_type":"markdown","source":"How many segments do we have?"},{"metadata":{"trusted":true,"_uuid":"2a819bb0c0c1b4ab734a768243739ba9a622c601"},"cell_type":"code","source":"len(test_files)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1b6cf8b77cb61c61b703b3e067fdadc91fc83486"},"cell_type":"markdown","source":"Does this match with the number of seg_ids in the sample submission?"},{"metadata":{"trusted":true,"_uuid":"38da3737512aaae1d43ec59fa14a8afe74a0eae3"},"cell_type":"code","source":"sample_submission = pd.read_csv(\"../input/sample_submission.csv\")\nsample_submission.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5b734eb0e0d3eac95bc7aa426f3b8c3f82c86827"},"cell_type":"code","source":"len(sample_submission.seg_id.values)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"73226056634bbdb29f4ea22be62635956aff23e0"},"cell_type":"markdown","source":"Ok. How does the signal of the test data look like?"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"0fec497da40633de95566b5a0c300ba552bea45d"},"cell_type":"code","source":"fig, ax = plt.subplots(4,1, figsize=(20,25))\n\nfor n in range(4):\n    seg = pd.read_csv(test_path  + test_files[n])\n    ax[n].plot(seg.acoustic_data.values, c=\"mediumseagreen\")\n    ax[n].set_xlabel(\"Index\")\n    ax[n].set_ylabel(\"Signal\")\n    ax[n].set_ylim([-300, 300])\n    ax[n].set_title(\"Test {}\".format(test_files[n]));","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7056a02a965a7f63bf75337c4bb0a37ac378b58d"},"cell_type":"markdown","source":"### Take-Away\n\n* These test segment examples differ a lot in the occurences of small peaks that seem to be similar to those in the train data before and after the heavy signal peak that occured some time before the lab earthquake took place. \n* They probably came from the same experiment but do neither form a continuous signal nor directly follow after the train data. "},{"metadata":{"_uuid":"ea46ab4cda6301fe610612b6ac4269c3ed1967cb"},"cell_type":"markdown","source":"## A question collection\n\nEven though there is only one signal feature and one target column, there is so much we need to understand and explore. Before we add more and more features, it's probably better to work with what is given so far. I'm afraid of feeling puzzled too fast. ;-) And to prevent that feeling even further, I like to collect the questions that draw circles in my mind:\n\n### Questions for the train set\n\n* Why do we have this low resolution jumps and why are they different? Is there some periodicity that may correlate with the signal? Would it be helpful to reconstruct that for test segments?\n* Why do we only have 16 earthquakes and such a high resolution of signal inbetween?\n\n\n### Questions for the test set\n\n* Are all segments of the test set of the same length?\n* Are they similar in their distributions or in the strength and time period between strong peaks?\n* Can we find some groups of similar test set segments?\n"},{"metadata":{"trusted":true,"_uuid":"376786419dddcb03ba51be8131070d854fae408e"},"cell_type":"markdown","source":"## A jump into train explorations"},{"metadata":{"trusted":true,"_uuid":"23c4add75254752bf097fb1af8b67a2c0ae83eb5"},"cell_type":"code","source":"train.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d7a2b169f4fdb29bfe32bbe4e224eea0225316ea"},"cell_type":"markdown","source":"* We can see that the mean is shifted towards higher values due to the earthquake. In addition we can see that the 25% up to 75% quartils are looking very discrete.\n* Looking at the quaketime we can't say much about it."},{"metadata":{"_uuid":"5c41dedd5253b78996c424552b1f3cc9125e6f36"},"cell_type":"markdown","source":"### The train signal distribution"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"38973e6d701a80b5f88242c0d89e78c4408ebf57"},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(20,5))\nsns.distplot(train.signal.values, ax=ax[0], color=\"Red\", bins=100, kde=False)\nax[0].set_xlabel(\"Signal\")\nax[0].set_ylabel(\"Density\")\nax[0].set_title(\"Signal distribution\")\n\nlow = train.signal.mean() - 3 * train.signal.std()\nhigh = train.signal.mean() + 3 * train.signal.std() \nsns.distplot(train.loc[(train.signal >= low) & (train.signal <= high), \"signal\"].values,\n             ax=ax[1],\n             color=\"Orange\",\n             bins=150, kde=False)\nax[1].set_xlabel(\"Signal\")\nax[1].set_ylabel(\"Density\")\nax[1].set_title(\"Signal distribution without peaks\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b40ef6537cb8917c4bf95d20a46c9edc318646ea"},"cell_type":"markdown","source":"### Take-Away\n\n* In the first plot the we can see that **the earthquake causes very high outlier values**. This way we can't say much about the signal distribution close to zero. \n* The second plot shows the signal at the median and mean around 4. We can see that it **looks very gaussian and balanced. It seems that the signal is somehow discrete.**"},{"metadata":{"_uuid":"de5f043d56f72f54c361ea5d4f82d05a0f3cac18"},"cell_type":"markdown","source":"### The stepsize\n\nBy computing the difference between to quaketimes we obtain some kind of stepsize that is probably equal within blocks and can show us the jump strength between different quaketimes. The following code computes differences first and drops the last row of train such that we can add the stepsize to the data. I think we won't loose fruitful information this way."},{"metadata":{"_kg_hide-input":false,"trusted":true,"_uuid":"34177c0b9415954f1ec609be0697389777d1ce6d"},"cell_type":"code","source":"stepsize = np.diff(train.quaketime)\ntrain = train.drop(train.index[len(train)-1])\ntrain[\"stepsize\"] = stepsize\ntrain.head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"deddf4f8c5f02057b7070ef0cf4f98e4759f0aa9"},"cell_type":"code","source":"train.stepsize = train.stepsize.apply(lambda l: np.round(l, 10))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1d1278916c8e5d07a0af08c3987a7cb299f4ec50"},"cell_type":"code","source":"stepsize_counts = train.stepsize.value_counts()\nstepsize_counts","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dba857ec8dd3d23675986530569716931d66c0f5"},"cell_type":"markdown","source":"Still strange! **The stepsize within two quaketimes is often given by either -1, -1,1 or -2 ns**. In addition we have **stepsizes close to -1,0955 ms or -0,9955 ms** :-) Do you see it? It's something centered at -1 ms. As we have one earthquake in train the **heavy stepsize of 11.5 is the step between the earthquake and a new cycle** that goes on until the next earthquake occurs."},{"metadata":{"_uuid":"c0f4ca8b8c97de8987d5efdad2062b416f9f1533"},"cell_type":"markdown","source":"## Setting up a validation strategy\n\nOk, I think that is part is difiicult. In my opinion it does not make sense to split the rows of the train data randomly to obtain validation data. This is temporal data and in the test segments we probably don't have any future values given. In the description you can see that:\n\n> The training data is a single, continuous segment of experimental data. The test data consists of a folder containing many small segments. The data within each test file is continuous, but the test files do not represent a continuous segment of the experiment; thus, the predictions cannot be assumed to follow the same regular pattern seen in the training file.\n\nHmm. What does that mean for us? Let's collect some scenarios:\n\n### Scenario A\n\n* The **cycles in the train data are independent of each other**. \n* If there was already an earthquake **does not change the material** in its fault behaviors. \n* Let's suppose the **test data uses the same experimental setup**. \n\nIf this is true it's not important that we need a second (or multiple) experiment to generate test data. We could only do one experiment and cut out segments to produce different data snippets. As this does not fit well to the description of the data **this scenario is not likely**.\n\n### Scenario B\n\n* The cycles depend on each other. There is temporal correlation of future signals with past ones. \n* The **material changes somehow from cycle to cylce** and perhaps future cycles are clearly different from past ones.\n* The **test data** uses the **same experimental setup** but to prevent leakage there were **at least one new experiment** done to produce segments that may depend on each other, e.g some of them may have temporal correlations. \n* To make it more tricky there are **only some snippets given in the test data whereas others are simply dropped**. \n\nIn this case we need a model that is able to capture the temporal dependence of cycles and it should be able to make nice predictions for different cycles. As we need to make predictions for several different cycle-phases it could be fruitful to use some kind of rolling window validation. \n\n### Scenario C\n\n* The **cylces depend on each other**. There is temporal correlation of future signals with past ones.\n* The state of the **material changes from cycle to cycle**.\n* The **test data uses experimental setups that are different from the train data generation** process. In the worst case we have multiple different experiments that produce test data. \n* To make it more tricky there are **only some snippets given in the test data whereas others are simply dropped**. \n\nThis is my personal worst case and this would probably cause high differences between scores of validation and test data. \n\n### How to split now?\n\nTo start with this competition **I prefer scenario B** and I haven't done it before but there is a scikit-learn implementation that may help us now. But before we keep going on... there is one more topic that should be considered: **The test set does not contain any earthquakes. Hence it's perhaps not a good idea to include them during training and validation**. In training they would cause extreme target outliers of quaketime. Our model may always try to match its predictions with these targets and this would hinder learning for the relevent parts, for those kind of predictions we need for the test set. "},{"metadata":{"trusted":true,"_uuid":"f1e541fb7ce6c00eb449f12a917b74020bd89369"},"cell_type":"code","source":"from sklearn.model_selection import TimeSeriesSplit\n\ncv = TimeSeriesSplit(n_splits=5)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25302279ada71973304868b8f8af866d21e39130"},"cell_type":"markdown","source":"## Rolling features\n\nI sometimes get stuck in too much details that may not be neccessary or are more fun to explore by yourself. As this is just a starter, I like to continue with some ideas and visualisations that have been used in initial work of Bertrand Rouet-Leduc and the LANL-group. One of the ideas was to use features extracted by a rolling window approach. Let's do the same and make some visualisations what goes on with these features until the first lab earthquake occurs.\n\n### Window size\n\nI don't know in adcance which kind of window size would be an appropriate choice and I think it's an hyperparameter we should try to optimize. But to start, let's try out some different sizes and the mean and standard deviation to select one that may be sufficient to play around:"},{"metadata":{"trusted":true,"_uuid":"066ee08ef873bfccbe743f7d701e63f5d5f256e2"},"cell_type":"code","source":"window_sizes = [10, 50, 100, 1000]\nfor window in window_sizes:\n    train[\"rolling_mean_\" + str(window)] = train.signal.rolling(window=window).mean()\n    train[\"rolling_std_\" + str(window)] = train.signal.rolling(window=window).std()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"1cc2e8b6affafe38880b16d3409f667faf6e79d4"},"cell_type":"code","source":"fig, ax = plt.subplots(len(window_sizes),1,figsize=(20,6*len(window_sizes)))\n\nn = 0\nfor col in train.columns.values:\n    if \"rolling_\" in col:\n        if \"mean\" in col:\n            mean_df = train.iloc[4435000:4445000][col]\n            ax[n].plot(mean_df, label=col, color=\"mediumseagreen\")\n        if \"std\" in col:\n            std = train.iloc[4435000:4445000][col].values\n            ax[n].fill_between(mean_df.index.values,\n                               mean_df.values-std, mean_df.values+std,\n                               facecolor='lightgreen',\n                               alpha = 0.5, label=col)\n            ax[n].legend()\n            n+=1\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"20e7d5bd68798e0ab8c05929537d8668ba7beb37"},"cell_type":"markdown","source":"A window size of 50 looks good enough to cover most fluctuations in sufficient detail without averaging out important signals. \n\n### Statistical features\n\nNow, that we have found a first choice for the window size, let's compute some basic rolling statistical features:\n\n* mean \n* standard deviation\n* 25% quartile\n* 50% quartile (median)\n* 75% quartile\n* interquartile range (75%-25%)\n* min\n* max\n* skewness \n* kurtosis\n\nJust for curiosity let's keep the mean and std for the window sizes we have excluded from above. Perhaps we can see later that these features were more important than the 50-size-window-features. "},{"metadata":{"trusted":true,"_uuid":"13bd2d78efd982e68e5db82b630ad0da2aca78aa"},"cell_type":"code","source":"train[\"rolling_q25\"] = train.signal.rolling(window=50).quantile(0.25)\ntrain[\"rolling_q75\"] = train.signal.rolling(window=50).quantile(0.75)\ntrain[\"rolling_q50\"] = train.signal.rolling(window=50).quantile(0.5)\ntrain[\"rolling_iqr\"] = train.rolling_q75 - train.rolling_q25\ntrain[\"rolling_min\"] = train.signal.rolling(window=50).min()\ntrain[\"rolling_max\"] = train.signal.rolling(window=50).max()\ntrain[\"rolling_skew\"] = train.signal.rolling(window=50).skew()\ntrain[\"rolling_kurt\"] = train.signal.rolling(window=50).kurt()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"95ac0cee2c6005bda0262a4a92fecc38421f511c"},"cell_type":"markdown","source":"## The end ;-) \n\nI hope my kernel has supported you in getting started with the data and perhaps it has also provided some ideas to try out. If you like it, you can make me happy with a comment and/or an upvote. Thank you! :-)"}],"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}