{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import pyarrow.parquet as pq # convert parguet formatted files\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy.signal import *\nimport statsmodels.api as sm\nfrom scipy import fftpack # Fast Fourier Transform functions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f089822757c0ac9111d1ac6fc1611b247d805b19"},"cell_type":"code","source":"# plot settings\nrand_seed = 135\nnp.random.seed(rand_seed)\nxsize = 12.0\nysize = 8.0\n\nfrom pylab import plot, show, savefig, xlim, figure, \\\n                hold, ylim, legend, boxplot, setp, axes","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35e1653752054258e33afe2b5b16fbec0270315b"},"cell_type":"markdown","source":"# Introduction\nPowerlines are measured using voltage, but you don't just measure voltage like you measure height. It's not just a single number. It goes up and down making waves. Each wave is a cycle and you can make judgements of each cycle by taking a multitude of measurements over time. For this competition, the underlying electric grid operates at 50 Hz (AKA 50 cycles per 1000 milliseconds). In other words, if we want to measure 50 cycles of voltage on this electric grid we would take measurements for 1000 milliseconds. The measurements in the competition were only performed for 20 milliseconds so we will only see one cycle. See basic math below:"},{"metadata":{"_uuid":"7efb98e52bd6eaa53dd30f4718061bf7df0a819d"},"cell_type":"markdown","source":"(50 cycles/ 1000 milliseconds) x (20 milliseconds) = 1 cycle"},{"metadata":{"_uuid":"2093996fdfe72270d3024b2c296ed872084645de"},"cell_type":"markdown","source":"In other words, the 800,000 measurements make up a time-course showing the voltage over time. "},{"metadata":{"_uuid":"56cba136d866948238298e3c450aafbc4395727a"},"cell_type":"markdown","source":"Now that we have that figured out, what is a 3-phase power scheme and what does that mean for us in this kaggle competition? To better understand, let's look at the data... "},{"metadata":{"trusted":true,"_uuid":"939399755169feaf275f9de986a4b338afe4794f"},"cell_type":"code","source":"%%time\n\ntrain_meta_df = pd.read_csv(\"../input/metadata_train.csv\")\ntrain_df = pq.read_pandas(\"../input/train.parquet\").to_pandas()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d126a8d6c94517fc019223beb0b36bfdc0e1a001"},"cell_type":"code","source":"train_meta_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"79a916a751791338c0d9b9631be0805bdb9e1a2a"},"cell_type":"code","source":"train_meta_df.head(n=9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eab72034b261532aa173b88faa8680711f2f47c2"},"cell_type":"code","source":"train_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3d6f8a41d70d2ff7d296f21bce12a6ebe84a9840"},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"580bda1c767029b01ed0249e112a9c22612fc3e1"},"cell_type":"markdown","source":"I just have book-keeping variables to help with transforming the voltage signals later..."},{"metadata":{"trusted":true,"_uuid":"29bd73c51976fd2204b178552b69a8e67dab9d82"},"cell_type":"code","source":"# sampling rate\nnum_samples = train_df.shape[0] # 800,000 samples per signal\nperiod = 0.02 # over a 20ms period\nfs = num_samples / period # 40MHz sampling rate\n\n# time array support\nt = np.array([i / fs for i in range(num_samples)])\n\n# frequency vector from FFT\nfreqs = fftpack.fftfreq(num_samples, d=1/fs)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"66f66c0356dbfd7b63de1bafe1acb693c1528411"},"cell_type":"markdown","source":"We have two training set files:"},{"metadata":{"_uuid":"69ad9d4465052f91b4c521cd634bd658348ac1e8"},"cell_type":"markdown","source":"* `metadata_train.csv` which contains four columns: \n * `signal_id`: a unique identifier so its meaningless\n  * `id_measurement`: ID code for each powerline. There should be 3 of each number which represents each phase in the 3-phase power scheme (AKA the trio)\n  * `phase`: the phase ID within the trio (0,1,2). \n  * `target`: 0 fixed, 1 broken"},{"metadata":{"_uuid":"1d46c75129f89aa42da543360e2dee79781e61e6"},"cell_type":"markdown","source":"* `train.parquet` which contains the 800,000 rows for the 800,000 measurements for each respective `signal_id`"},{"metadata":{"_uuid":"671f1af2cb5944d7df43ffa85f5f338bd8f59249"},"cell_type":"markdown","source":"Since each powerline has 3 rows of data (1 for each phase) in the `metadata_train.csv` and the `train.parquet` file that means the `target` is the same for each `signal_id` that has the same `id_measurement` (see plot below which shows the target variable for each phase seperately).  It also means we have 800,000 x 3 measurements per powerline."},{"metadata":{"trusted":true,"_uuid":"75715ed42cf3e75bbdf6135ce2d7156fc532ff03"},"cell_type":"code","source":"%%time\n\nfig, ax = plt.subplots()\nfig.set_size_inches(xsize, ysize)\n\nax =sns.countplot(x=\"phase\", hue=\"target\", data=train_meta_df, ax=ax)\nax.set_title(\"Distributions of `Target` variable for each phase is equal\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5218d63853edc37f04a25952f4463806953fe73d"},"cell_type":"markdown","source":"We can also note from the plot above that there is a big target class imbalancement issue which could be dealt with via subsampling or loss functions, but I digress."},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"47eb28678ef8f0cfd4be6dc89055eb275b8ea098"},"cell_type":"code","source":"%%time\n\nplt.figure(figsize=(15, 10))\nplt.title(\"ID measurement:0, Target:0\",\n         fontdict={'fontsize':36})\nplt.plot(train_df[\"0\"].values, marker=\"o\", label='Phase 0')\nplt.plot(train_df[\"1\"].values, marker=\"o\", label='Phase 1')\nplt.plot(train_df[\"2\"].values, marker=\"o\", label='Phase 2')\nplt.ylim(-50,50)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"616c6a68a0b2fbe614e5635246bd1b120af8bbfb"},"cell_type":"code","source":"%%time\n\nplt.figure(figsize=(15, 10))\nplt.title(\"ID measurement:1, Target:1\",\n         fontdict={'fontsize':36})\nplt.plot(train_df[\"3\"].values, marker=\"o\", label='Phase 0')\nplt.plot(train_df[\"4\"].values, marker=\"o\", label='Phase 1')\nplt.plot(train_df[\"5\"].values, marker=\"o\", label='Phase 2')\nplt.ylim(-50,50)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7440926c706b7cfcc9feb0612f7a1f4d7f55266e"},"cell_type":"code","source":"%%time\n\nplt.figure(figsize=(15, 10))\nplt.title(\"ID measurement:2, Target:0\",\n         fontdict={'fontsize':36})\nplt.plot(train_df[\"6\"].values, marker=\"o\", label='Phase 0')\nplt.plot(train_df[\"7\"].values, marker=\"o\", label='Phase 1')\nplt.plot(train_df[\"8\"].values, marker=\"o\", label='Phase 2')\nplt.ylim(-50,50)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a6abee81e6be440687632f53f4d2756bb19fcf95"},"cell_type":"markdown","source":"So based on our tiny sampling of these three powerlines it looks like amplitude could be a significant factor in predicting whether these powerlines are faulty."},{"metadata":{"_uuid":"862776543278c4a2e048073b8b1615d6bbb2a3e2"},"cell_type":"markdown","source":"According to the wiki on the 3-phase power scheme:\n> In a symmetric three-phase power supply system, three conductors each carry an alternating current of the same frequency and voltage amplitude relative to a common reference but with a phase difference of one third of a cycle between each."},{"metadata":{"_uuid":"0907769e04d55b082e0eb36d6adec5ceead85e7c"},"cell_type":"markdown","source":"Therefore, as competitors we need to look into smoothing techniques to measure the amplitude of these waves. We should also measure the variance of our estimated amplitude for each powerline across the 3 phases. "},{"metadata":{"_uuid":"6a13db3986e9b695ddf9220b9ce91300d573273a"},"cell_type":"markdown","source":"# Feature Engineering\nSince we don't want to do machine learning on the raw numbers from the parquet files (which could take forever and may not even be useful) I want to create features I can add to the meta tables to train on instead. We will also needs to add these features to the test set."},{"metadata":{"_uuid":"b38f36ba814e062a2d5b4b40e99f497c0d2e96f8"},"cell_type":"markdown","source":"Note that to make the code run faster I subset the training data for when I am editing this kernel but when I commit the code I do not do this..."},{"metadata":{"trusted":true,"_uuid":"f80e7d4238a3e194ec3354d1c9cfb9d30954b651"},"cell_type":"code","source":"# uncomment to subset the data (Note the graphs look really different when you do this)\n#train_subset_df = train_df.iloc[:,range(0,99)]\n#train_subset_meta_df = train_meta_df.iloc[range(0,99),:]\n\n# uncomment to use the full dataset\ntrain_subset_df = train_df\ntrain_subset_meta_df = train_meta_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"220351509b26817918c12a1c2439c0873ceb7bb8"},"cell_type":"markdown","source":"## Mean, Median, and Standard Deviation of Measurements\nMultiple kernels have done this and found these numbers to be slightly useful so let's see..."},{"metadata":{"trusted":true,"_uuid":"4603cd056d422ff3aa465cc568b612a54fab7ac8"},"cell_type":"code","source":"%%time\n\nmean_list = train_subset_df.apply(np.mean)\nmedian_list = train_subset_df.apply(np.median)\nstd_list = train_subset_df.apply(np.std)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2ee2e55553070acf27fb100b60ce5d5c3c337b1e"},"cell_type":"code","source":"mean_signal_df = mean_list.to_frame()\nmean_signal_df = mean_signal_df.reset_index()\nmean_signal_df = mean_signal_df.drop(\"index\",axis=1)\ntrain_subset_meta_df =train_subset_meta_df.merge(mean_signal_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"mean\"})\n\nmedian_signal_df = median_list.to_frame()\ntrain_subset_meta_df =train_subset_meta_df.merge(median_signal_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"median\"})\n\nstd_signal_df = std_list.to_frame()\ntrain_subset_meta_df =train_subset_meta_df.merge(std_signal_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"std_dev\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2671384a71acd307df292780aebe6d267c7ea877"},"cell_type":"code","source":"train_subset_meta_df.head(n=9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"81613f8aab6e10445453390a8ce3642bca939fd1"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Mean Across Target Variables\")\nax = sns.boxplot(x=\"target\", y=\"mean\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d01809f5bf54f25a0a34203fa070c207a98412cd"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Median Across Target Variables\")\nax = sns.boxplot(x=\"target\", y=\"median\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a3b4fff1fcab7822257c4128c80c1c55786fb92"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Stardard Deviation Across Target Variables\")\nax = sns.boxplot(x=\"target\", y=\"std_dev\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3c9cdf725db4aafef78293636dfbbd2d4300b155"},"cell_type":"markdown","source":"## Amplitude of Rolling Series\nLets smooth the plots of the first 3 power lines to see if we can see any differences between these waves in target=0 verses target=1"},{"metadata":{"trusted":true,"_uuid":"f428e22dfc648058f47a9a49110757fa7e818943"},"cell_type":"code","source":"ts1 = train_df[\"0\"]\nts2 = train_df[\"1\"]\nts3 = train_df[\"2\"]\n\nplt.figure(figsize=(16,6))\nplt.title(\"ID measurement:0, Target:0\",\n         fontdict={'fontsize':36})\nplt.plot(ts1.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts1.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts2.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts2.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts3.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts3.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7eaaa902e0780cc2ba3c47ba66bca93cb95dab6f"},"cell_type":"code","source":"ts1 = train_df[\"3\"]\nts2 = train_df[\"4\"]\nts3 = train_df[\"5\"]\n\nplt.figure(figsize=(16,6))\nplt.title(\"ID measurement:1, Target:1\",\n         fontdict={'fontsize':36})\nplt.plot(ts1.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts1.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts2.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts2.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts3.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts3.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"acae79ae72c662983f15b6e596ed79ecfb392330"},"cell_type":"code","source":"ts1 = train_df[\"6\"]\nts2 = train_df[\"7\"]\nts3 = train_df[\"8\"]\n\nplt.figure(figsize=(16,6))\nplt.title(\"ID measurement:2, Target:0\",\n         fontdict={'fontsize':36})\nplt.plot(ts1.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts1.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts2.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts2.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.plot(ts3.rolling(window=100000,center=False).mean(),label='Rolling Mean');\nplt.plot(ts3.rolling(window=100000,center=False).std(),label='Rolling sd');\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b81a8ed6224695629e11e216c59073077cdb1730"},"cell_type":"markdown","source":"It doesn't look like we can see any difference, but this is a small sample size to work with. Let's look at the amplitude across each target group. To calculate the amplitude, I smooth the powerline signals to create a single wave then I subtract the lowest and highest point. "},{"metadata":{"trusted":true,"_uuid":"ad1fcf8ee21f8250ef8896a81984326d22708ecb"},"cell_type":"code","source":"%%time\n\ndef calc_rolling_amp(row, window=100000):\n    return np.max(row.rolling(window,center=False).mean()) - np.min(row.rolling(window=100000,center=False).mean())\n\nrolling100k_amp = train_subset_df.apply(calc_rolling_amp)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b147b70901e632dc82842204fb91604248b006ba"},"cell_type":"code","source":"rolling100k_amp_df = rolling100k_amp.to_frame()\ntrain_subset_meta_df =train_subset_meta_df.merge(rolling100k_amp_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"rolling100k_amp\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"da8edfa6f878e54ff32096bffe976fba6447db96"},"cell_type":"code","source":"train_subset_meta_df.head(n=9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f99f1dae51bc64a85c356952df7bbef2a7d27f54"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Amplitude Across Target Variables\")\nax = sns.boxplot(x=\"target\", y=\"rolling100k_amp\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"174b574d5aba87f227c646220a0ec207a5b6fb90"},"cell_type":"markdown","source":"It looks like there is not much of a difference unfortunately..."},{"metadata":{"_uuid":"dfa3caee12daf02d18f8eea669bb4259f3146b9a"},"cell_type":"markdown","source":"\n## Measuring Amount of Noisy Points\n### Number of points 1SD from the mean\nThe next feature I am intersted in looking at is the number of data points in each signal that is greater than 1 SD from the mean."},{"metadata":{"trusted":true,"_uuid":"9821c9dbc7f0f4f6e8677bf51e37a36b7a704016"},"cell_type":"code","source":"def count1SDfromTheMean(row):\n    max_1sd = np.mean(row) + np.std(row)\n    min_1sd = np.mean(row) - np.std(row)\n    noise_points = [x for x in row if (x > max_1sd) or (x < min_1sd)]\n    return (len(noise_points))\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"66743e0d981feefbc1a8f1170790aa74a30b65a2","scrolled":true},"cell_type":"code","source":"%%time\n\ncount1SDfromTheMean_list = train_subset_df.apply(count1SDfromTheMean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8c0d2ade47acab4bb9d91675a49021afb009bf50"},"cell_type":"code","source":"#%%time\n#count1SDfromTheMean_list = train_df.apply(count1SDfromTheMean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"88bf2cafb605e4ff7d1a3ed3bb0f4114e986c28f"},"cell_type":"code","source":"count1SDfromTheMean_df = count1SDfromTheMean_list.to_frame()\ntrain_subset_meta_df =train_subset_meta_df.merge(count1SDfromTheMean_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"count1SDfromTheMean\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"09e1e976526e500b8bfdee67664baa8c624a8a0d"},"cell_type":"code","source":"train_subset_meta_df.head(n=9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"281f3b4295b5177ce420f1930784e3e3a9e8089e"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Noise Count Across Target Variables\")\nax =sns.boxplot(x=\"target\", y=\"count1SDfromTheMean\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d51685306dc0cd68813c772df292e02676dfd3f0"},"cell_type":"markdown","source":"There is not that much of a difference and spoiler alert the next parameter is probably better and gives basically the same information. I am going to drop this column."},{"metadata":{"trusted":true,"_uuid":"e3464360bcb4c6a48cfd3c6dd79b76b75655919f"},"cell_type":"code","source":"# drop \"count1SDfromTheMean\" from the training set\ntrain_subset_meta_df = train_subset_meta_df.drop([\"count1SDfromTheMean\"], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e936939b83beb49751c35113fc01a7a3e5eedc4d"},"cell_type":"markdown","source":"\n### Number of points 2SD from the mean\nThe next feature I am intersted in looking at is the number of data points in each signal that is greater than 2 SD from the mean."},{"metadata":{"trusted":true,"_uuid":"1d9cf9a82c3e36c054e919aa94c4f5f67b97c239"},"cell_type":"code","source":"def count2SDfromTheMean(row):\n    max_1sd = np.mean(row) + (2 * np.std(row))\n    min_1sd = np.mean(row) - (2 * np.std(row))\n    noise_points = [x for x in row if (x > max_1sd) or (x < min_1sd)]\n    return (len(noise_points))\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"01f3f632cd2361120c53354227709228ff17ae5a"},"cell_type":"code","source":"%%time\n\ncount2SDfromTheMean_list = train_subset_df.apply(count2SDfromTheMean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9894399c98bff82982e0b03159e080ba528d613b"},"cell_type":"code","source":"count2SDfromTheMean_df = count2SDfromTheMean_list.to_frame()\ntrain_subset_meta_df =train_subset_meta_df.merge(count2SDfromTheMean_df,\"inner\", \n                    left_index=True,right_index=True)\ntrain_subset_meta_df = train_subset_meta_df.rename(index=str, columns={0:\"count2SDfromTheMean\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"86c432ffe451d6047dd9c59d55839da49ce49580"},"cell_type":"code","source":"train_subset_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2ceb323b104b74745d2ab5a0e0a91e5ec59848cd"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Noise Count Across Target Variables\")\nax =sns.boxplot(x=\"target\", y=\"count2SDfromTheMean\", data=train_subset_meta_df)\nplt.ylim(0,250)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"795be5498e1f76250e55172c3fb7e22d57cfd694"},"cell_type":"markdown","source":"This looks like it could be a really useful feature!"},{"metadata":{"_uuid":"be3d0f14f9c4236b5cededf2d8701d43df4c2e1b"},"cell_type":"markdown","source":"1. ## Sum \"extremeness\" of each bin\nI want to try and bin the 800,000 data points and see if noise in a specfic location is significantly different in the different targets, but to do this I need to transform the waves so that they can be referenced equally. See kernel for explanation of how this transformation works:  https://www.kaggle.com/fernandoramacciotti/sync-waves-with-fft-coeffs"},{"metadata":{"trusted":true,"_uuid":"1f3d880e8b11908dc869eb0b95f6090e305db46a"},"cell_type":"code","source":"n_signals_to_load = 3\nsignals = pq.read_pandas(\n    '../input/train.parquet', \n    columns=[str(i) for i in range(n_signals_to_load)]).to_pandas()\nsignals.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f1a8e6e40a4cb0501288c5bc59ef62865c50316e"},"cell_type":"code","source":"# get fft coeffs\ndef get_fft_coeffs(sig):\n    return fftpack.fft(sig)\n\n# get coeff with highest norm\ndef get_highest_coeff(fft_coeffs, freqs, verbose=True):\n    coeff_norms = np.abs(fft_coeffs) # get norms (fft coeffs are complex)\n    max_idx = np.argmax(coeff_norms)\n    max_coeff = fft_coeffs[max_idx] # get max coeff\n    max_freq = freqs[max_idx] # assess which is the dominant frequency\n    max_amp = (coeff_norms[max_idx] / num_samples) * 2 # times 2 because there are mirrored freqs\n    if verbose:\n        print('Dominant frequency is {:,.1f}Hz with amplitude of {:,.1f}\\n'.format(max_freq, max_amp))\n    \n    return max_coeff, max_amp, max_freq\n\n# get max coeff phase\ndef get_max_coeff_phase(max_coeff):\n    return np.angle(max_coeff)\n\n# construct the instant angular phase vector indexed by pi, i.e. ranges from 0 to 2\ndef get_instant_w(time_vector, f0, phase_shift):\n    w_vector = 2 * np.pi * time_vector * f0 + phase_shift\n    w_vector_norm = np.mod(w_vector / (2 * np.pi), 1) * 2 # range between cycle of 0-2 \n    return w_vector, w_vector_norm\n\n# find index of chosen phase to align\ndef get_align_idx(w_vector_norm, align_value=0.5):\n    candidates = np.where(np.isclose(w_vector_norm, align_value))\n    # since we are in discrete time, threre could be many values close to the desired one\n    # so let's take the one in the middle\n    return int(np.median(candidates))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d303242a07245257edcf39991a374fcc5e1ab211"},"cell_type":"code","source":"# align waves with np.roll()\nalign_phase = 0.5 # w_i = pi/2\n\nfig = plt.figure(figsize=(12, 9))\nplot_number = 0\n\nfor signal_id in signals.columns:\n    # get samples\n    sig = signals[signal_id]\n    \n    # fft\n    fft_coeffs = get_fft_coeffs(sig)\n    \n    # asses dominant frequency\n    max_coeff, amp, f0 = get_highest_coeff(fft_coeffs, freqs, verbose=True)\n    \n    # phase shift\n    ps = get_max_coeff_phase(max_coeff)\n    \n    # get angular phase vector\n    w, w_norm = get_instant_w(t, f0, ps)\n    \n    # generate dominant signal at f0\n    dominant_wave = amp * np.cos(w)\n    \n    # idx to roll\n    origin = get_align_idx(w_norm, align_value=align_phase)\n    \n    # roll signal and dominant wave\n    sig_rolled = np.roll(sig, num_samples - origin)\n    dominant_wave_rolled = np.roll(dominant_wave, num_samples - origin)\n    \n    # plot signals\n    plot_number += 1\n    ax = fig.add_subplot(3, 1, plot_number)\n    \n    ax.plot(t * 1000, sig_rolled, label='Rolled Original') # original signal\n    ax.plot(t * 1000, dominant_wave_rolled, color='red', label='Rolled Wave at {:.0f}Hz'.format(f0)) # wave at f0\n    ax.legend()\n    ax.set_xlabel('time (ms)')\n    ax.set_ylabel('Amplitude')\n    ax.set_title('Signal {} rolled'.format(signal_id))\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b34831ca5b117303eeb3e8fde346fecac1e03dff"},"cell_type":"markdown","source":"[](http://)What if you get the sum of the difference between the rolled original and rolledwave at 50 Hz in 125 bins?"},{"metadata":{"trusted":true,"_uuid":"49676d93fed833254e91fe7e90943a2db99d6580"},"cell_type":"code","source":"# align waves with np.roll()\nalign_phase = 0.5 # w_i = pi/2\nnbins=125\ntrain_subset_meta_df.index = pd.RangeIndex(start=0, stop=len(train_subset_meta_df), step=1)\n\n\ndiff_df= pd.DataFrame()\nfor signal_id in train_subset_df.columns:\n    # get samples\n    sig = train_subset_df[signal_id]\n    # fft\n    fft_coeffs = get_fft_coeffs(sig)\n    \n    # asses dominant frequency\n    max_coeff, amp, f0 = get_highest_coeff(fft_coeffs, freqs, verbose=False)\n    \n    # phase shift\n    ps = get_max_coeff_phase(max_coeff)\n    \n    # get angular phase vector\n    w, w_norm = get_instant_w(t, f0, ps)\n    \n    # generate dominant signal at f0\n    dominant_wave = amp * np.cos(w)\n    \n    # idx to roll\n    origin = get_align_idx(w_norm, align_value=align_phase)\n    \n    # roll signal and dominant wave\n    sig_rolled = np.roll(sig, num_samples - origin)\n    dominant_wave_rolled = np.roll(dominant_wave, num_samples - origin)\n    \n    diff_bw_signAndDom = np.abs(dominant_wave_rolled-sig_rolled)\n    sum_signals=[]\n    numSignalsInBin=int(num_samples/nbins)\n    #print(num_samples)\n    #print(numSignalsInBin)\n    for i in range(0,num_samples,numSignalsInBin):\n        bin_sum = np.sum(diff_bw_signAndDom[i:i+numSignalsInBin])\n        sum_signals.append(bin_sum)\n    diff_df = diff_df.append(pd.Series(sum_signals), ignore_index=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"624b97ca50c977a3be59f721393e13505168a05c"},"cell_type":"code","source":"colnamesOfDiffTable=[\"rolldiff\"+str(x) for x in list(diff_df.columns)]\ndiff_df.columns = colnamesOfDiffTable\ntrain_subset_meta_df =train_subset_meta_df.merge(diff_df,\"inner\", \n                    left_index=True,right_index=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5e50371f65f776aecdb25bb12ed58dd09cbf4f67"},"cell_type":"code","source":"train_subset_meta_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"56452b81bbdc15b8b4f52382903b6ebbf99af0f7"},"cell_type":"code","source":"# function for setting the colors of the box plots pairs\ndef setBoxColorsOfTargets(bp):\n    setp(bp['boxes'][0], color='blue')\n    setp(bp['caps'][0], color='blue')\n    setp(bp['caps'][1], color='blue')\n    setp(bp['whiskers'][0], color='blue')\n    setp(bp['whiskers'][1], color='blue')\n    setp(bp['fliers'][0], color='blue')\n    setp(bp['fliers'][1], color='blue')\n    setp(bp['medians'][0], color='blue')\n\n    setp(bp['boxes'][1], color='orange')\n    setp(bp['caps'][2], color='orange')\n    setp(bp['caps'][3], color='orange')\n    setp(bp['whiskers'][2], color='orange')\n    setp(bp['whiskers'][3], color='orange')\n    setp(bp['fliers'][2], color='orange')\n    setp(bp['fliers'][3], color='orange')\n    setp(bp['medians'][1], color='orange')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d38d7e21cfb05f6a1a356c33c695d35e2f36bac3"},"cell_type":"code","source":"np.arange(0, 25, step=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e2639703882b9e416b4c0349be63824142147c82"},"cell_type":"code","source":"for i in range(0,len(colnamesOfDiffTable),25):\n    rolldiff_train_df = train_subset_meta_df.loc[:,[\"target\"]+colnamesOfDiffTable[i:i+25]].copy()\n    rolldiff_train_df = rolldiff_train_df.set_index(\"target\")\n    rolldiff_train_df = rolldiff_train_df.stack()\n    rolldiff_train_df = rolldiff_train_df.reset_index()\n    rolldiff_train_df.columns=[\"target\",\"rolldiff\",\"value\"]\n\n    plt.figure(figsize=(6,8))\n    sns.set(style=\"whitegrid\")\n    plt.title(\"Sum of Extremeness across Bins in Target Variables\")\n    ax = sns.boxplot(x=\"rolldiff\",y=\"value\", hue=\"target\", data=rolldiff_train_df)\n    plt.xticks(np.arange(0, 25, step=1),np.arange(i, i+25, step=1))\n    plt.xlabel(\"bin\")\n    plt.show()\n    plt.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6ee80e90d21dc491579729a3a763a1a161c6062c"},"cell_type":"code","source":"plt.figure(figsize=(6,8))\nsns.set(style=\"whitegrid\")\nplt.title(\"Sum of Extremeness across Bins in Target Variables\")\nax = sns.boxplot(x=\"target\", y=\"std_dev\", data=train_subset_meta_df)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d4ef724a8a44bae78e892b93d2c126f70894fecc"},"cell_type":"markdown","source":" Maybe next we can actually try and predict using these features..."},{"metadata":{"trusted":true,"_uuid":"a13d58ab983a3d978ef4f2c05b553df6abafaba3"},"cell_type":"code","source":"train_subset_meta_df.to_csv('metadata_train_V2.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"298a20349263f0bdbcd7d7ea9690a3af18e2d716"},"cell_type":"markdown","source":"# Resources\n* https://www.kaggle.com/timothycwillard/vsb-power-line-faults-eda-feature-engineering\n* https://www.kaggle.com/theoviel/fast-fourier-transform-denoising"},{"metadata":{"trusted":true,"_uuid":"9ef1f9c8385242274127c24c0f252d653875b17b"},"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}