{"cells":[{"metadata":{"trusted":true,"_uuid":"5226a95ba6e6ea70ec33bfe1334afdbfdb3d62e3"},"cell_type":"code","source":"\"\"\"\nLargely followed methodology outlined in this paper:\n\"A Complex Classification Approach of Partial Discharges from Covered Conductors in Real Environment\"\n2017\nS. Mišák, J. Fulneček, T. Vantuch, T. Buriánek and T. Ježowicz\nPreprint version available here:\nhttps://www.dropbox.com/s/2ltuvpw1b1ms2uu/A%20Complex%20Classification%20Approach%20of%20Partial%20Discharges%20from%20Covered%20Conductors%20in%20Real%20Environment%20%28preprint%29.pdf?dl=0\n\nMore details of the methodology from the paper above are provided in this thesis:\n\"Analysis of Time Series Data\"\nTomáš Vantuch\n2018\nAvailable here:\nhttp://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n\nBegan with kernel from:\nhttps://www.kaggle.com/jackvial/dwt-signal-denoising\nto start with Discrete Wavelet Transform denoising/peak identification\n\"\"\"\n\n# Import Libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nimport pyarrow.parquet as pq\nimport gc\nimport pywt\nfrom statsmodels.robust import mad\nimport scipy\nfrom scipy import signal\nfrom scipy.signal import butter, find_peaks\nfrom scipy import fftpack # Fast Fourier Transform functions\nimport os.path\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import confusion_matrix, make_scorer, matthews_corrcoef\nfrom sklearn.model_selection import cross_val_score\nimport warnings\n\n# Suppress pandas future warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\n# Define data directory\ndata_dir = '../input'\n\n# Print scipy version - 1.1.0 at time of creation\nprint(scipy.__version__)","execution_count":2,"outputs":[{"output_type":"stream","text":"1.1.0\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"# Define Functions\n\n**Calculate Mean Absolute Deviation** \nmaddest()\n\n**Synchronise All Waveforms** \nsync_phase()\n\n**High Pass Filter** \nhigh_pass_filter()\n\n**Discrete Wavelet Transform Denoising** \ndenoise_signal()\n\n**Cancel False Peaks** \ncancel_false_peaks()\n\n**Extract Peak and Valley Features:** \npv_features()"},{"metadata":{"trusted":true},"cell_type":"code","source":"# 800,000 data points taken over 20 ms\n# Grid operates at 50hz, 0.02 * 50 = 1, so 800k samples in 20 milliseconds will capture one complete cycle\nn_samples = 800000\n\n# Sample duration is 20 miliseconds\nsample_duration = 0.02\n\n# Sample rate is the number of samples in one second\n# Sample rate will be 40mhz\nsample_rate = n_samples * (1 / sample_duration)\n\n# time array support\nt = np.array([i / sample_rate for i in range(n_samples)])\n\n# frequency vector for FFT\nfreqs = fftpack.fftfreq(n_samples, d=1/sample_rate)\n\n# Mean Absolute Deviation\ndef maddest(d, axis=None):\n return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\n# Synchronise all waveforms\n# Adapted from https://www.kaggle.com/fernandoramacciotti/sync-waves-with-fft-coeffs\ndef sync_phase(x, align_value=0.5):\n \n # fft\n fft_coeffs = fftpack.fft(x)\n \n # asses dominant frequency\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] / n_samples) * 2 # times 2 because there are mirrored freqs\n if max_freq != 50:\n # print('Dominant frequency is {:,.1f}Hz with amplitude of {:,.1f}\\n'.format(max_freq, max_amp))\n # print('Signal ID is {:,.1f} and its target is {:,.1f}\\n'.format(metadata_train.signal_id[i], metadata_train.target[i]))\n return x, max_freq\n else:\n \n # phase shift\n phase_shift = np.angle(max_coeff)\n \n # get angular phase vector\n w = 2 * np.pi * t * max_freq + phase_shift\n w_norm = np.mod(w / (2 * np.pi), 1) * 2 # range between cycle of 0-2\n \n # idx to roll\n candidates = np.where(np.isclose(w_norm, align_value))\n # since we are in discrete time, threre could be many values close to the desired one\n # so take the median\n origin = int(np.median(candidates))\n \n # roll/sync signal\n sig_rolled = np.roll(x, n_samples - origin)\n return sig_rolled, max_freq\n\n# High Pass Filter\n# Adapted from https://github.com/randxie/Kaggle-VSB-Baseline/blob/master/src/utils/util_signal.py\ndef high_pass_filter(x, low_cutoff=1000, sample_rate=sample_rate): # 1000 is default, correctly set to 10^4 below \n \n # nyquist frequency (half the sample rate)\n nyquist = 0.5 * sample_rate\n norm_low_cutoff = low_cutoff / nyquist\n \n # Fault pattern usually exists in high frequency band. \n # According to literature, the pattern is visible above 10^4 Hz.\n # For digital filters, Wn is normalized from 0 to 1, \n # where 1 is the Nyquist frequency, pi radians/sample (Wn is thus in half-cycles / sample).\n sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n filtered_sig = signal.sosfilt(sos, x)\n\n return filtered_sig\n\n# Discrete Wavelet Transform Denoising\n# 1. Adapted from waveletSmooth function found here:\n# http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n# 2. Threshold equation and using hard mode in threshold as mentioned\n# in section '3.2 denoising based on optimized singular values' from the Vantuch thesis \ndef denoise_signal(x, wavelet='db4', level=1):\n\n # Decompose to get the wavelet coefficients\n coeff = pywt.wavedec(x, wavelet, mode=\"per\" )\n \n # Calculate sigma for threshold as defined in the Vantuch thesis, using Mean Absolute Deviation\n sigma = (1/0.6745) * maddest( coeff[-level] ) # eqn 3.8 \n\n # Calculte the univeral threshold\n uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) ) # eqn 3.9\n coeff[1:] = (pywt.threshold( i, value=uthresh, mode='hard' ) for i in coeff[1:] )\n \n # Reconstruct the signal using the thresholded coefficients\n return pywt.waverec(coeff, wavelet, mode='per')\n\n# Cancel False Peaks\n# Such as those from corona discharges\n# Adapted from https://www.kaggle.com/jeffreyegan/vsb-power-line-fault-detection-approach\ndef cancel_false_peaks(signal, peak_indexes, min_height_fp=50):\n\n false_peak_indexes = []\n min_height_fp = min_height_fp\n max_sym_distance = 20 \n max_pulse_train = 750 \n max_height_ratio = 1.25 \n min_height_ratio = 0.25\n for pk in range(len(peak_indexes)-1):\n if not peak_indexes[pk] in false_peak_indexes:\n if (signal[peak_indexes[pk]] > min_height_fp and signal[peak_indexes[pk+1]] < 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:\n # if 1 > abs(signal[peak_indexes[pk+1]])/abs(signal[peak_indexes[pk]]) > max_height_ratio:\n if max_height_ratio > abs(signal[peak_indexes[pk+1]])/abs(signal[peak_indexes[pk]]) > min_height_ratio:\n for x in range(len(peak_indexes)):\n if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train:\n false_peak_indexes.append(peak_indexes[x]) \n \n if (signal[peak_indexes[pk]] < -min_height_fp and signal[peak_indexes[pk+1]] > 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:\n if max_height_ratio > abs(signal[peak_indexes[pk+1]])/abs(signal[peak_indexes[pk]]) > min_height_ratio:\n for x in range(len(peak_indexes)):\n if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train:\n false_peak_indexes.append(peak_indexes[x])\n \n true_peak_indexes = list(set(peak_indexes) - set(false_peak_indexes))\n \n return np.array(true_peak_indexes, dtype=np.int32), np.array(false_peak_indexes, dtype=np.int32)\n\n\n# Extract Peak and Valley Features\ndef pv_features(x_dn, pv_true, peaks, valleys, rel_height=0.2):\n\n # p=peaks, v=valleys, h=height, w=width\n true_peaks = np.array(list(set(peaks) & set(pv_true)), dtype=np.int32) # intersection of true peaks/valleys and peaks\n if true_peaks.shape[0] == 0:\n p_n = 0 \n ph_mean = 0\n ph_max = 0\n pw_mean = 0\n pw_max = 0\n else:\n p_n = true_peaks.shape[0]\n peak_heights, _ , _ = signal.peak_prominences(x_dn, true_peaks)\n peak_heights = np.vstack((x_dn[true_peaks], peak_heights)).min(axis=0)\n ph_mean = peak_heights.mean()\n ph_max = peak_heights.max()\n peak_widths, peak_width_heights, peak_left_ips , peak_right_ips = signal.peak_widths(x_dn, true_peaks, rel_height=rel_height)\n pw_mean = peak_widths.mean()\n pw_max = peak_widths.max()\n \n true_valleys = np.array(list(set(valleys) & set(pv_true)), dtype=np.int32) # intersection of true peaks/valleys and valleys\n if true_valleys.shape[0] == 0:\n v_n = 0 \n vh_mean = 0\n vh_max = 0\n vw_mean = 0\n vw_max = 0\n else:\n v_n = true_valleys.shape[0]\n valley_heights, _ , _ = signal.peak_prominences(-x_dn, true_valleys)\n valley_heights = np.vstack((-x_dn[true_valleys], valley_heights)).min(axis=0)\n vh_mean = valley_heights.mean()\n vh_max = valley_heights.max()\n valley_widths, valley_width_heights, valley_left_ips , valley_right_ips = signal.peak_widths(-x_dn, true_valleys, rel_height=rel_height)\n vw_mean = valley_widths.mean()\n vw_max = valley_widths.max()\n \n return np.array([p_n, ph_mean, ph_max, pw_mean, pw_max, v_n, vh_mean, vh_max, vw_mean, vw_max])\n ","execution_count":3,"outputs":[]},{"metadata":{"_uuid":"d9a8df3c33016b94ae63d223de4cb652a5fe9dd9"},"cell_type":"markdown","source":"### Import Train Metadata"},{"metadata":{"trusted":true,"_uuid":"5cfda79d0648204193833b9a5aa1d4e8f4c90666","scrolled":true},"cell_type":"code","source":"metadata_train = pd.read_csv(data_dir + '/metadata_train.csv')\nmetadata_train.head()","execution_count":4,"outputs":[{"output_type":"execute_result","execution_count":4,"data":{"text/plain":" signal_id id_measurement phase target\n0 0 0 0 0\n1 1 0 1 0\n2 2 0 2 0\n3 3 1 0 1\n4 4 1 1 1","text/html":"
\n\n
\n \n
\n
\n
signal_id
\n
id_measurement
\n
phase
\n
target
\n
\n \n \n
\n
0
\n
0
\n
0
\n
0
\n
0
\n
\n
\n
1
\n
1
\n
0
\n
1
\n
0
\n
\n
\n
2
\n
2
\n
0
\n
2
\n
0
\n
\n
\n
3
\n
3
\n
1
\n
0
\n
1
\n
\n
\n
4
\n
4
\n
1
\n
1
\n
1
\n
\n \n
\n
"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"**Check number of phases on which a fault occurs** \nThe majority of faults (target=1) occur on all 3 phases"},{"metadata":{"trusted":true,"_uuid":"4e1a09820640ce2c05f46615d22b5e9484c673e7"},"cell_type":"code","source":"metadata_train.groupby([\"id_measurement\"]).sum()['target'].value_counts()","execution_count":5,"outputs":[{"output_type":"execute_result","execution_count":5,"data":{"text/plain":"0 2710\n3 156\n1 19\n2 19\nName: target, dtype: int64"},"metadata":{}}]},{"metadata":{"_uuid":"f88a87ebd7e6f322caad76783f54055e99b23cbd"},"cell_type":"markdown","source":"### Import Small Subset of Train Data for Exploration"},{"metadata":{"trusted":true,"_uuid":"4524f1ddb1414f8072a7f3d07189957ad6518532","scrolled":false},"cell_type":"code","source":"subset_train = pq.read_pandas(data_dir + '/train.parquet', columns=[str(i) for i in range(21)]).to_pandas()\nsubset_train.shape","execution_count":6,"outputs":[{"output_type":"execute_result","execution_count":6,"data":{"text/plain":"(800000, 21)"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"### Plot the Target Variable Counts by Phase"},{"metadata":{"trusted":true,"_uuid":"1212477ae2d75474717fe2e0e36039ce146466cd"},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))\nsns.countplot(x=\"target\", data=metadata_train, ax=ax1)\nsns.countplot(x=\"target\", data=metadata_train, hue=\"phase\", ax=ax2)","execution_count":7,"outputs":[{"output_type":"execute_result","execution_count":7,"data":{"text/plain":""},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"
","image/png":"iVBORw0KGgoAAAANSUhEUgAAA00AAAEKCAYAAADOwb7RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu0XXV16PHvNOGlggQ45GJObLCmSFAbyOFh4XpVSgi5vYSi8vCRCLHp1SB47agGRluQh4OOoggoOFKJEBViRCnBy4CmQapyRTiBFAKUJoVAzhmBHJIY8cErzvvH/iXdxJzNzmGv8/x+xthjrzXXb601l574c671278VmYkkSZIkacdeN9AJSJIkSdJgZtEkSZIkSQ1YNEmSJElSAxZNkiRJktSARZMkSZIkNWDRJEmSJEkNWDRJkiRJUgMWTZIkSZLUgEWTJEmSJDUweqATqMJ+++2XEyZMGOg0JGlEW758+bOZ2TbQeQxG9lOSNDg021cNy6JpwoQJdHZ2DnQakjSiRcSTA53DYGU/JUmDQ7N9lcPzJEmSJKkBiyZJkiRJasCiSZIkSZIaGJa/aZIkSZJUvZdeeomuri6ef/75gU6lod1335329nZ22WWXPu1v0SRJkiSpT7q6uthzzz2ZMGECETHQ6exQZrJhwwa6uro48MAD+3QMh+dJkiRJ6pPnn3+efffdd9AWTAARwb777vuanoZZNEmSJEnqs8FcMG31WnO0aJIkSZKkBiyaJEmSJPWbCRMm8Oyzzw50GjvFiSAamPLXCwc6BQ1By/9h5kCnIElDwlMXvnOgU+izt/zdQwOdgqR+VGnRFBH/B/gEkMBDwBnAAcAiYF9gOfCxzHwxInYDFgJTgA3AqZm5phznXGA2sAU4OzPvqDJvSZKGgqF+c+/mPQc6g747+qqjBzqFPrv703cPdAoaIdasWcO0adOYMmUK999/P4cccggLF9b+d+uqq67i1ltv5aWXXuJ73/seb3/727n33ns555xzeP7559ljjz345je/yUEHHcTDDz/MGWecwYsvvsjvfvc7vv/97zNx4kS+/e1vc+WVV/Liiy9y5JFHcvXVVzNq1KhKrqWy4XkRMQ44G+jIzHcAo4DTgL8HLs/MtwGbqBVDlO9NJX55aUdETCr7HQJMA66OiGr+05AkSZLUMo899hif+tSnePTRR9lrr724+uqrAdhvv/24//77+eQnP8lll10GwNvf/nZ+8pOf8MADD3DhhRdy3nnnAfD1r3+dc845hxUrVtDZ2Ul7ezuPPvoo3/3ud7n77rtZsWIFo0aN4jvf+U5l11H18LzRwB4R8RLwemAd8H7gw2X79cAFwDXAjLIMcBPw1ahNczEDWJSZLwBPRMRq4AjgZxXnLkmSJOk1GD9+PEcfXXsy+9GPfpQrr7wSgJNPPhmAKVOm8IMf/ACAzZs3M2vWLFatWkVE8NJLLwHw7ne/m0suuYSuri5OPvlkJk6cyLJly1i+fDmHH344AL/97W/Zf//9K7uOyp40ZWY3cBnwFLViaTO14Xi/yMyXS7MuYFxZHgesLfu+XNrvWx/fwT7bRMSciOiMiM6enp7WX5AkSZKknbL9VN9b13fbbTcARo0axcsv10qDv/3bv+V973sfK1eu5NZbb932XqUPf/jDLFmyhD322IPp06dz5513kpnMmjWLFStWsGLFCh577DEuuOCCyq6jyuF5Y6g9JToQeDPwBmrD6yqRmfMzsyMzO9ra2qo6jSRJkqQmPfXUU/zsZ7UBYjfccAPHHHNMr203b97MuHG1ZyPXXXfdtvjjjz/OW9/6Vs4++2xmzJjBgw8+yLHHHstNN93E+vXrAdi4cSNPPvlkZddR5ZTjfwo8kZk9mfkS8APgaGDviNg6LLAd6C7L3cB4gLL9TdQmhNgW38E+kiRJkgapgw46iK997WscfPDBbNq0iU9+8pO9tv3c5z7Hueeey6GHHrrt6RPA4sWLecc73sHkyZNZuXIlM2fOZNKkSVx88cVMnTqVd73rXRx33HGsW7eusuuo8jdNTwFHRcTrgd8CxwKdwI+AD1KbQW8WcEtpv6Ss/6xsvzMzMyKWADdExJepPbGaCNxbYd6SJEmSWmD06NF8+9vffkVszZo125Y7Ojq46667gNpvl/7jP/5j27aLL74YgHnz5jFv3rzfO/app57Kqaee2vqkd6Cyoikzfx4RNwH3Ay8DDwDzgf8LLIqIi0vs2rLLtcC3ykQPG6nNmEdmPhwRi4FHynHmZuaWqvKWJEmSpHqVzp6XmecD528Xfpza7Hfbt30e+FAvx7kEuKTlCUqSJEmqxIQJE1i5cuVAp9ESVf6mSZIkSZKGPIsmSZIkSWrAokmSJEmSGrBokiRJkqQGKp0IQpKkwSoixgMLgbFAAvMz84qIuAD4C6CnND0vM28r+5wLzAa2AGdn5h0lPg24AhgFfCMzL+3Pa5GkwWLKXy9s6fGW/8PMptrdfvvtnHPOOWzZsoVPfOITO5yi/LWwaJIkjVQvA3+VmfdHxJ7A8ohYWrZdnpmX1TeOiEnUXodxCLX3Bv5LRPxR2fw14DigC7gvIpZk5iP9chWSNMJt2bKFuXPnsnTpUtrb2zn88MM58cQTmTRpUsvO4fA8SdKIlJnrMvP+svwc8CgwrsEuM4BFmflCZj4BrKb2Co0jgNWZ+Xhmvkjt5e0zqs1ekrTVvffey9ve9jbe+ta3suuuu3Laaadxyy23tPQcFk2SpBEvIiYAhwI/L6GzIuLBiFgQEWNKbBywtm63rhLrLS5J6gfd3d2MHz9+23p7ezvd3d0tPYdFkyRpRIuINwLfBz6Tmb8ErgH+EJgMrAO+1KLzzImIzojo7OnpefUdJEmDhkWTJGnEiohdqBVM38nMHwBk5jOZuSUzfwf8I7XhdwDdwPi63dtLrLf4K2Tm/MzsyMyOtra21l+MJI1Q48aNY+3a/3rg39XVxbhxrX3gb9EkSRqRIiKAa4FHM/PLdfED6pr9ObCyLC8BTouI3SLiQGAicC9wHzAxIg6MiF2pTRaxpD+uQZIEhx9+OKtWreKJJ57gxRdfZNGiRZx44oktPYez50mSRqqjgY8BD0XEihI7Dzg9IiZTm4Z8DfCXAJn5cEQsBh6hNvPe3MzcAhARZwF3UJtyfEFmPtyfFyJJg0WzU4S30ujRo/nqV7/K8ccfz5YtWzjzzDM55JBDWnuOlh5NkqQhIjN/CsQONt3WYJ9LgEt2EL+t0X6SpGpNnz6d6dOnV3Z8h+dJkiRJUgMWTZIkSZLUgEWTJEmSJDVQWdEUEQdFxIq6zy8j4jMRsU9ELI2IVeV7TGkfEXFlRKwuLxQ8rO5Ys0r7VRExq6qcJUmSJGl7lRVNmflYZk7OzMnAFOA3wM3APGBZZk4ElpV1gBOoTd86EZhD7eWCRMQ+wPnAkdTelXF+3dvZJUmSJKlS/TU871jgPzPzSWAGcH2JXw+cVJZnAAuz5h5g7/KujOOBpZm5MTM3AUuBaf2UtyRJkqQRrr+mHD8NuLEsj83MdWX5aWBsWR4HrK3bp6vEeotLkiRJGkSeuvCdLT3eW/7uoVdtc+aZZ/LDH/6Q/fffn5UrV75q+76o/ElTeTv6icD3tt+WmUnt5YGtOM+ciOiMiM6enp5WHFKSJEnSIPfxj3+c22+/vdJz9MfwvBOA+zPzmbL+TBl2R/leX+LdwPi6/dpLrLf4K2Tm/MzsyMyOtra2Fl+CJEmSpMHoPe95D/vss0+l5+iPoul0/mtoHsASYOsMeLOAW+riM8ssekcBm8swvjuAqRExpkwAMbXEJEmSJKlylf6mKSLeABwH/GVd+FJgcUTMBp4ETinx24DpwGpqM+2dAZCZGyPiIuC+0u7CzNxYZd6SJEmStFWlRVNm/hrYd7vYBmqz6W3fNoG5vRxnAbCgihwlSZIkqZH+mnJckiRJkoak/ppyXJIkSdIw18wU4a12+umnc9ddd/Hss8/S3t7OF77wBWbPnt3Sc1g0SZIkSRqybrzxxldv9Bo5PE+SJEmSGrBokiRJkqQGLJokSZIkqQGLJkmSJElqwKJJkiRJkhqwaJIkSZKkBpxyXJIkSVJLHH3V0S093t2fvvtV26xdu5aZM2fyzDPPEBHMmTOHc845p6V5WDRJkiRJGrJGjx7Nl770JQ477DCee+45pkyZwnHHHcekSZNadg6H50mSJEkasg444AAOO+wwAPbcc08OPvhguru7W3oOiyZJkiRJw8KaNWt44IEHOPLII1t6XIsmSZIkSUPer371Kz7wgQ/wla98hb322qulx7ZokiRJkjSkvfTSS3zgAx/gIx/5CCeffHLLj2/RJEmSJGnIykxmz57NwQcfzGc/+9lKzlHp7HkRsTfwDeAdQAJnAo8B3wUmAGuAUzJzU0QEcAUwHfgN8PHMvL8cZxbwN+WwF2fm9VXmLUmSJGnnNTNFeMvPeffdfOtb3+Kd73wnkydPBuCLX/wi06dPb9k5qp5y/Arg9sz8YETsCrweOA9YlpmXRsQ8YB7weeAEYGL5HAlcAxwZEfsA5wMd1Aqv5RGxJDM3VZy7JGkYi4jxwEJgLLX+ZX5mXlH6HW/uSdIQccwxx5CZlZ6jsuF5EfEm4D3AtQCZ+WJm/gKYAWztTK4HTirLM4CFWXMPsHdEHAAcDyzNzI2lUFoKTKsqb0nSiPEy8FeZOQk4CpgbEZOo3cxblpkTgWVlHV55c28OtZt71N3cOxI4Ajg/Isb054VIkqpV5W+aDgR6gG9GxAMR8Y2IeAMwNjPXlTZPU7vDBzAOWFu3f1eJ9RaXJKnPMnPd1idFmfkc8Ci1/sWbe5KkV6iyaBoNHAZck5mHAr/mv+7WAZC152gteZYWEXMiojMiOnt6elpxSEnSCBERE4BDgZ/jzT1J2ilVD41rhdeaY5VFUxfQlZk/L+s3USuinil35ijf68v2bmB83f7tJdZb/BUyc35mdmRmR1tbW0svRJI0fEXEG4HvA5/JzF/Wb/PmniQ1tvvuu7Nhw4ZBXThlJhs2bGD33Xfv8zEqmwgiM5+OiLURcVBmPgYcCzxSPrOAS8v3LWWXJcBZEbGI2rjwzZm5LiLuAL5YNz58KnBuVXlLkkaOiNiFWsH0ncz8QQk/ExEHlD6o2Zt7790uftf258rM+cB8gI6OjsH7/y4kaSe0t7fT1dXFYL8ZtPvuu9Pe3t7n/auePe/TwHfKzHmPA2dQe7q1OCJmA08Cp5S2t1GbkWg1tVmJzgDIzI0RcRFwX2l3YWZurDhvSdIwV2bDuxZ4NDO/XLdpCd7ck6Sm7LLLLhx44IEDnUblKi2aMnMFtanCt3fsDtomMLeX4ywAFrQ2O0nSCHc08DHgoYhYUWLnUSuWvLknSdqm6idNkiQNSpn5UyB62ezNPUnSNlVOBCFJkiRJQ55FkyRJkiQ1YNEkSZIkSQ1YNEmSJElSAxZNkiRJktSARZMkSZIkNWDRJEmSJEkNWDRJkiRJUgMWTZIkSZLUgEWTJEmSJDVg0SRJkiRJDVg0SZIkSVIDFk2SJEmS1IBFkyRJkiQ1YNEkSZIkSQ1YNEmSJElSA5UWTRGxJiIeiogVEdFZYvtExNKIWFW+x5R4RMSVEbE6Ih6MiMPqjjOrtF8VEbOqzFmSJEmS6vXHk6b3ZebkzOwo6/OAZZk5EVhW1gFOACaWzxzgGqgVWcD5wJHAEcD5WwstSZIkSaraQAzPmwFcX5avB06qiy/MmnuAvSPiAOB4YGlmbszMTcBSYFp/Jy1JkiRpZKq6aErgnyNieUTMKbGxmbmuLD8NjC3L44C1dft2lVhv8VeIiDkR0RkRnT09Pa28BkmSJEkj2OiKj39MZnZHxP7A0oj49/qNmZkRka04UWbOB+YDdHR0tOSYkiRJklTpk6bM7C7f64Gbqf0m6Zky7I7yvb407wbG1+3eXmK9xSVJkiSpcpUVTRHxhojYc+syMBVYCSwBts6ANwu4pSwvAWaWWfSOAjaXYXx3AFMjYkyZAGJqiUmSJElS5aocnjcWuDkitp7nhsy8PSLuAxZHxGzgSeCU0v42YDqwGvgNcAZAZm6MiIuA+0q7CzNzY4V5S5IkSdI2lRVNmfk48Mc7iG8Ajt1BPIG5vRxrAbCg1TlKkiRJ0qtpanheRCxrJiZJ0kCwn5IkVanhk6aI2B14PbBf+T1RlE17sYNpvyVJ6k/2U5Kk/vBqw/P+EvgM8GZgOf/VGf0S+GqFeUmS1Az7KUlS5RoWTZl5BXBFRHw6M6/qp5wkSWqK/ZQkqT80NRFEZl4VEX8CTKjfJzMXVpSXJElN60s/FRELgD8D1mfmO0rsAuAvgJ7S7LzMvK1sOxeYDWwBzs7MO0p8GnAFMAr4RmZe2tKLkyQNuKaKpoj4FvCHwApqnQVAAhZNkqQB18d+6jpqQ/i2b3N5Zl623fEnAacBh1AbCvgvEfFHZfPXgOOALuC+iFiSmY/0/WokSYNNs1OOdwCTyrTgkiQNNjvdT2XmjyNiQpPNZwCLMvMF4ImIWA0cUbatLq/ZICIWlbYWTZI0jDQ15TiwEvhvVSYiSdJr0Mp+6qyIeDAiFpQZ+aA2E9/aujZdJdZb/PdExJyI6IyIzp6enh01kSQNUs0+adoPeCQi7gVe2BrMzBMryUqSpJ3Tqn7qGuAiakP7LgK+BJzZigQzcz4wH6Cjo8ORG5I0hDRbNF1QZRKSJL1GF7TiIJn5zNbliPhH4IdltRsYX9e0vcRoEJckDRPNzp73r1UnIklSX7Wqn4qIAzJzXVn9c2rD/gCWADdExJepTQQxEbiX2nuhJkbEgdSKpdOAD7ciF0nS4NHs7HnPURuqALArsAvw68zcq6rEJElqVl/6qYi4EXgvsF9EdAHnA++NiMnlWGuovTyXzHw4IhZTm+DhZWBuZm4pxzkLuIPalOMLMvPhll+gJGlANfukac+tyxER1GYGOqqqpCRJ2hl96acy8/QdhK9t0P4S4JIdxG8Dbms6WUnSkNPs7HnbZM0/AcdXkI8kSa+J/ZQkqdWaHZ53ct3q66i9D+P5SjKSJGkn2U9JkqrU7Ox5/6tu+WVq47xntDwbSZL6xn5KklSZZn/TdEZfTxARo4BOoDsz/6zMMLQI2BdYDnwsM1+MiN2AhcAUYANwamauKcc4F5gNbAHOzsw7+pqPJGn4eS39lCRJr6ap3zRFRHtE3BwR68vn+xHR3uQ5zgEerVv/e+DyzHwbsIlaMUT53lTil5d2RMQkalO4HgJMA64uhZgkScBr7qckSWqo2YkgvkntHRVvLp9bS6yh0mH9T+AbZT2A9wM3lSbXAyeV5RllnbL92LoZkBZl5guZ+QSwGjiiybwlSSNDn/opSZKa0WzR1JaZ38zMl8vnOqCtif2+AnwO+F1Z3xf4RWa+XNa7gHFleRywFqBs31zab4vvYB9JkqDv/ZQkSa+q2aJpQ0R8NCJGlc9Hqf3uqFcR8WfA+sxc/pqzbEJEzImIzojo7Onp6Y9TSpIGj53upyRJalazRdOZwCnA08A64IPAx19ln6OBEyNiDbWJH94PXAHsHRFbJ6BoB7rLcjcwHqBsfxO1Dm9bfAf7bJOZ8zOzIzM72tq8uShJI0xf+ilJkprSbNF0ITArM9syc39qndMXGu2QmedmZntmTqA2kcOdmfkR4EfUOjOAWcAtZXlJWadsvzMzs8RPi4jdysx7E4F7m8xbkjQy7HQ/JUlSs5p9T9O7MnPT1pXM3BgRh/bxnJ8HFkXExcADwLUlfi3wrYhYDWykVmiRmQ9HxGLgEWrv3pibmVv6eG5J0vDUyn5KkqRXaLZoel1EjNnaIUXEPjuxL5l5F3BXWX6cHcx+l5nPAx/qZf9LgEuaPZ8kacR5Tf2UJEmNNNuhfAn4WUR8r6x/CIsYSdLgYT8lSapMU0VTZi6MiE5qkzkAnJyZj1SXliRJzbOfkiRVaWeG2D1C7XdFkiQNOvZTkqSqNDt7niRJkiSNSBZNkiRJktSARZMkSZIkNWDRJEmSJEkNWDRJkiRJUgMWTZIkSZLUgEWTJEmSJDVg0SRJkiRJDVg0SZIkSVIDFk2SJEmS1IBFkyRJkiQ1YNEkSZIkSQ1YNEmSJElSA5UVTRGxe0TcGxH/FhEPR8QXSvzAiPh5RKyOiO9GxK4lvltZX122T6g71rkl/lhEHF9VzpIkSZK0vSqfNL0AvD8z/xiYDEyLiKOAvwcuz8y3AZuA2aX9bGBTiV9e2hERk4DTgEOAacDVETGqwrwlSSNARCyIiPURsbIutk9ELI2IVeV7TIlHRFxZbuA9GBGH1e0zq7RfFRGzBuJaJEnVqqxoyppfldVdyieB9wM3lfj1wElleUZZp2w/NiKixBdl5guZ+QSwGjiiqrwlSSPGddRuxtWbByzLzInAsrIOcAIwsXzmANdArcgCzgeOpNY3nb+10JIkDR+V/qYpIkZFxApgPbAU+E/gF5n5cmnSBYwry+OAtQBl+2Zg3/r4DvaRJKlPMvPHwMbtwvU38La/sbew3BC8B9g7Ig4AjgeWZubGzNxEra/bvhCTJA1xlRZNmbklMycD7dTuwL29qnNFxJyI6IyIzp6enqpOI0ka3sZm5rqy/DQwtiz3dgOv6Rt79lOSNHT1y+x5mfkL4EfAu6ndnRtdNrUD3WW5GxgPULa/CdhQH9/BPvXnmJ+ZHZnZ0dbWVsl1SJJGjsxMasPKW3U8+ylJGqKqnD2vLSL2Lst7AMcBj1Irnj5Yms0CbinLS8o6ZfudpcNaApxWZtc7kNp48nuryluSNKI9U4bdUb7Xl3hvN/CaurEnSRraqnzSdADwo4h4ELiP2pjvHwKfBz4bEaup/Wbp2tL+WmDfEv8s5ce3mfkwsBh4BLgdmJuZWyrMW5I0ctXfwNv+xt7MMoveUcDmMozvDmBqRIwpE0BMLTFJ0jAy+tWb9E1mPggcuoP44+xg9rvMfB74UC/HugS4pNU5SpJGroi4EXgvsF9EdFGbBe9SYHFEzAaeBE4pzW8DplObwfU3wBkAmbkxIi6idnMQ4MLM3H5yCUnSEFdZ0SRJ0mCWmaf3sunYHbRNYG4vx1kALGhhapKkQaZfJoKQJEmSpKHKokmSJEmSGrBokiRJkqQGLJokSZIkqQGLJkmSJElqwKJJkiRJkhqwaJIkSZKkBiyaJEmSJKkBiyZJkiRJasCiSZIkSZIasGiSJEmSpAYsmiRJkiSpAYsmSZIkSWrAokmSJEmSGrBokiRJkqQGLJokSZIkqYHKiqaIGB8RP4qIRyLi4Yg4p8T3iYilEbGqfI8p8YiIKyNidUQ8GBGH1R1rVmm/KiJmVZWzJEmSJG2vyidNLwN/lZmTgKOAuRExCZgHLMvMicCysg5wAjCxfOYA10CtyALOB44EjgDO31poSZIkSVLVKiuaMnNdZt5flp8DHgXGATOA60uz64GTyvIMYGHW3APsHREHAMcDSzNzY2ZuApYC06rKW5IkSZLq9ctvmiJiAnAo8HNgbGauK5ueBsaW5XHA2rrdukqst/j255gTEZ0R0dnT09PS/CVJkiSNXJUXTRHxRuD7wGcy85f12zIzgWzFeTJzfmZ2ZGZHW1tbKw4pSZIkSdUWTRGxC7WC6TuZ+YMSfqYMu6N8ry/xbmB83e7tJdZbXJIkSZIqV+XseQFcCzyamV+u27QE2DoD3izglrr4zDKL3lHA5jKM7w5gakSMKRNATC0xSZIkSarc6AqPfTTwMeChiFhRYucBlwKLI2I28CRwStl2GzAdWA38BjgDIDM3RsRFwH2l3YWZubHCvCVJkiRpm8qKpsz8KRC9bD52B+0TmNvLsRYAC1qXnSRJkiQ1p19mz5MkSZKkocqiSZIkSZIasGiSJEmSpAYsmiRJ2k5ErImIhyJiRUR0ltg+EbE0IlaV7zElHhFxZUSsjogHI+Kwgc1ektRqFk2SJO3Y+zJzcmZ2lPV5wLLMnAgsK+sAJwATy2cOcE2/ZypJqpRFkyRJzZkBXF+WrwdOqosvzJp7gL23vsRdkjQ8WDRJkvT7EvjniFgeEXNKbGx56TrA08DYsjwOWFu3b1eJSZKGiSpfbitJ0lB1TGZ2R8T+wNKI+Pf6jZmZEZE7c8BSfM0BeMtb3tK6TCVJlfNJkyRJ28nM7vK9HrgZOAJ4Zuuwu/K9vjTvBsbX7d5eYtsfc35mdmRmR1tbW5XpS5JazKJJkqQ6EfGGiNhz6zIwFVgJLAFmlWazgFvK8hJgZplF7yhgc90wPknSMODwPEmSXmkscHNEQK2fvCEzb4+I+4DFETEbeBI4pbS/DZgOrAZ+A5zR/ylLkqpk0SRJUp3MfBz44x3ENwDH7iCewNx+SE2SNEAcnidJkiRJDVg0SZIkSVIDFk2SJEmS1IBFkyRJkiQ1UFnRFBELImJ9RKysi+0TEUsjYlX5HlPiERFXRsTqiHgwIg6r22dWab8qImbt6FySJEmSVJUqnzRdB0zbLjYPWJaZE4FlZR3gBGBi+cwBroFakQWcDxxJ7cWC528ttCRJkiSpP1RWNGXmj4GN24VnANeX5euBk+riC7PmHmDv8rb144GlmbkxMzcBS/n9QkySJEmSKtPfv2kaW/eW9KepvUAQYBywtq5dV4n1FpckSZKkfjFgE0GUlwFmq44XEXMiojMiOnt6elp1WEmSJEkjXH8XTc+UYXeU7/Ul3g2Mr2vXXmK9xX9PZs7PzI7M7Ghra2t54pIkSZJGpv4umpYAW2fAmwXcUhefWWbROwrYXIbx3QFMjYgxZQKIqSUmSZIkSf1idFUHjogbgfcC+0VEF7VZ8C4FFkfEbOBJ4JTS/DZgOrAa+A1wBkBmboyIi4D7SrsLM3P7ySUkSZIkqTKVFU2ZeXovm47dQdsE5vZynAXAghamJkmSJElNG7CJICRJkiRpKLBokiRJkqQGLJokSZIkqQGLJkmSJElqwKJJkiRJkhqwaJIkSZKkBiqbclySJEl6NVP+euFAp9Bny/9h5kCn0GdHX3X0QKfwmtz96bv79XwWTdIw99SF7xzoFDQEveXvHhroFCRp0BvSfeyYvQY6gyHF4XmSJEmS1IBFkyRJkiQ1YNEkSZIkSQ1YNEmSJElSAxZNkiRJktR7xPSCAAAEv0lEQVSARZMkSZIkNWDRJEmSJEkNWDRJkiRJUgNDpmiKiGkR8VhErI6IeQOdjyRJ9eynJGn4GhJFU0SMAr4GnABMAk6PiEkDm5UkSTX2U5I0vA2Jogk4AlidmY9n5ovAImDGAOckSdJW9lOSNIwNlaJpHLC2br2rxCRJGgzspyRpGBs90Am0SkTMAeaU1V9FxGMDmc8IsB/w7EAnMRjFZbMGOgU1z7/j3pwfrTjKH7TiIMOF/dTvq/gPxH/fvYizW/LvW/TL/8j5d9yLFv4dN/Vf41ApmrqB8XXr7SW2TWbOB+b3Z1IjWUR0ZmbHQOchvRb+HauF7KcGGf99azjw73jwGCrD8+4DJkbEgRGxK3AasGSAc5IkaSv7KUkaxobEk6bMfDkizgLuAEYBCzLz4QFOS5IkwH5Kkoa7IVE0AWTmbcBtA52HtnGIiYYD/47VMvZTg47/vjUc+Hc8SERmDnQOkiRJkjRoDZXfNEmSJEnSgLBo0k6LiGkR8VhErI6IeQOdj7SzImJBRKyPiJUDnYuk1rOf0lBnPzX4WDRpp0TEKOBrwAnAJOD0iJg0sFlJO+06YNpAJyGp9eynNExch/3UoGLRpJ11BLA6Mx/PzBeBRcCMAc5J2imZ+WNg40DnIakS9lMa8uynBh+LJu2sccDauvWuEpMkaTCwn5LUchZNkiRJktSARZN2Vjcwvm69vcQkSRoM7KcktZxFk3bWfcDEiDgwInYFTgOWDHBOkiRtZT8lqeUsmrRTMvNl4CzgDuBRYHFmPjywWUk7JyJuBH4GHBQRXRExe6BzktQa9lMaDuynBp/IzIHOQZIkSZIGLZ80SZIkSVIDFk2SJEmS1IBFkyRJkiQ1YNEkSZIkSQ1YNEmSJElSAxZNUoUiYu+I+FQ/nOe9EfEnVZ9HkjS82E9JzbFokqq1N9B0ZxQ1ffl3+V7AzkiStLPsp6Qm+J4mqUIRsQiYATwG/Ah4FzAG2AX4m8y8JSImUHsJ48+BKcB04E+BzwO/AP4NeCEzz4qINuDrwFvKKT4DdAP3AFuAHuDTmfmT/rg+SdLQZj8lNceiSapQ6Wh+mJnviIjRwOsz85cRsR+1DmQi8AfA48CfZOY9EfFm4P8BhwHPAXcC/1Y6oxuAqzPzpxHxFuCOzDw4Ii4AfpWZl/X3NUqShi77Kak5owc6AWkECeCLEfEe4HfAOGBs2fZkZt5Tlo8A/jUzNwJExPeAPyrb/hSYFBFbj7lXRLyxP5KXJA179lNSLyyapP7zEaANmJKZL0XEGmD3su3XTR7jdcBRmfl8fbCuc5Ikqa/sp6ReOBGEVK3ngD3L8puA9aUjeh+14Q47ch/wPyJiTBkq8YG6bf8MfHrrSkRM3sF5JElqlv2U1ASLJqlCmbkBuDsiVgKTgY6IeAiYCfx7L/t0A18E7gXuBtYAm8vms8sxHoyIR4D/XeK3An8eESsi4r9XdT2SpOHFfkpqjhNBSINQRLwxM39V7uDdDCzIzJsHOi9JksB+SiOPT5qkwemCiFgBrASeAP5pgPORJKme/ZRGFJ80SZIkSVIDPmmSJEmSpAYsmiRJkiSpAYsmSZIkSWrAokmSJEmSGrBokiRJkqQGLJokSZIkqYH/D5e2GSZgdpCqAAAAAElFTkSuQmCC\n"},"metadata":{}}]},{"metadata":{"_uuid":"d7922cec4a1654bf066b4dfa491f3e0ef3697822"},"cell_type":"markdown","source":"# Exploratory Plot 1\n**Plot the first six signals** (all 3 phases of the first 2 measurement IDs) \nMeasurement ID 1 - no faults \nMeasurement ID 2 -faults on all 3 phases \nEach row of plots has the original signal (left), high pass signal (middle), and high pass and DWT denoised signal (right). \nTwo rows are plotted for each signal: complete signal in the first row, and zoomed in signal in the second row. "},{"metadata":{"trusted":true,"_uuid":"a137e0bb2e22894ce620388954e24312224eb057","_kg_hide-input":false,"_kg_hide-output":false},"cell_type":"code","source":"train_length = 6\nfor i in range(train_length):\n signal_id = str(i)\n meta_row = metadata_train[metadata_train['signal_id'] == i]\n measurement = str(meta_row['id_measurement'].values[0])\n signal_id = str(meta_row['signal_id'].values[0])\n phase = str(meta_row['phase'].values[0])\n \n # Synchronize waveforms\n x_sync, dominant_frequency = sync_phase(subset_train[signal_id])\n \n # Apply high pass filter with low cutoff of 10kHz, this will remove the low frequency 50Hz sinusoidal motion in the signal\n x_hp = high_pass_filter(x_sync, low_cutoff=10000, sample_rate=sample_rate)\n \n # Apply denoising\n x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n \n slice_size = 10000\n font_size = 16\n \n fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(30, 10))\n \n ax[0, 0].plot(x_sync, alpha=0.5)\n ax[0, 0].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[0, 0].legend(['Original'], fontsize=font_size)\n \n # Show smaller slice of the signal to get a better idea of the effect the high pass frequency filter is having on the signal\n ax[1, 0].plot(x_sync[:slice_size], alpha=0.5)\n ax[1, 0].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[1, 0].legend([f\"Original n: {slice_size}\"], fontsize=font_size)\n \n ax[0, 1].plot(x_hp, 'r', alpha=0.5)\n ax[0, 1].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[0, 1].legend(['HP filter'], fontsize=font_size)\n ax[1, 1].plot(x_hp[:slice_size], 'r', alpha=0.5)\n ax[1, 1].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[1, 1].legend([f\"HP filter n: {slice_size}\"], fontsize=font_size)\n \n ax[0, 2].plot(x_dn, 'g', alpha=0.5)\n ax[0, 2].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[0, 2].legend(['HP filter and denoising'], fontsize=font_size)\n ax[1, 2].plot(x_dn[:slice_size], 'g', alpha=0.5)\n ax[1, 2].set_title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}\", fontsize=font_size)\n ax[1, 2].legend([f\"HP filter and denoising n: {slice_size}\"], fontsize=font_size)\n \n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"73bbb0d538e304a841b7c4cd457c180a24207f90"},"cell_type":"markdown","source":"# Exploratory Plot 2\n**Plot the first 21 signals (filtered and denoised) zoomed in on 2000 out of 800000 points to see peak detail and canceled peaks** \npeaks (and valleys) marked with a blue x \ncanceled false peaks marked with a red x \npeak widths and heights indicated by an orange line \nfaults on all three phases of measurement ID 1 only \n\n"},{"metadata":{"trusted":true,"_uuid":"4274b7f9f430097544518303568e195d03805bab","scrolled":false},"cell_type":"code","source":"train_length = subset_train.shape[1]\n\nfor i in range(train_length):\n signal_id = str(i)\n meta_row = metadata_train[metadata_train['signal_id'] == i]\n measurement = str(meta_row['id_measurement'].values[0])\n signal_id = str(meta_row['signal_id'].values[0])\n phase = str(meta_row['phase'].values[0])\n target = str(meta_row['target'].values[0])\n \n subset_train_row = subset_train[signal_id]\n \n # Apply high pass filter with low cutoff of 10kHz, this will remove the low frequency 50Hz sinusoidal motion in the signal\n x_hp = high_pass_filter(subset_train_row, low_cutoff=10000, sample_rate=sample_rate)\n \n # Apply denoising\n x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n \n start = 116000\n slice_size = 2000\n \n # Find peaks\n peaks, peak_properties = find_peaks(x_dn[start:start+slice_size], height=0.1, width=0, rel_height=0.2)\n valleys, valley_properties = find_peaks(-x_dn[start:start+slice_size], height=0.1, width=0, rel_height=0.2)\n all_peaks = np.sort(np.concatenate((peaks, valleys)))\n \n # Cancel false peaks\n true_pv, false_pv = cancel_false_peaks(x_dn[start:start+slice_size], all_peaks, min_height_fp=15)\n true_peaks = np.array(list(set(peaks) & set(true_pv)), dtype=np.int32) # intersection of true peaks/valleys and peaks\n true_valleys = np.array(list(set(valleys) & set(true_pv)), dtype=np.int32) # intersection of true peaks/valleys and valleys \n \n # Calculate heights and widths of true peaks\n peak_heights, _ , _ = signal.peak_prominences(x_dn[start:start+slice_size], true_peaks)\n peak_heights = np.vstack((x_dn[start+true_peaks], peak_heights)).min(axis=0)\n peak_widths, peak_width_heights, peak_left_ips , peak_right_ips = signal.peak_widths(x_dn[start:start+slice_size], true_peaks, rel_height=0.2)\n \n valley_heights, _ , _ = signal.peak_prominences(-x_dn[start:start+slice_size], true_valleys)\n valley_heights = np.vstack((-x_dn[start+true_valleys], valley_heights)).min(axis=0)\n valley_widths, valley_width_heights, valley_left_ips , valley_right_ips = signal.peak_widths(-x_dn[start:start+slice_size], true_valleys, rel_height=0.2)\n \n # Plot\n font_size = 16\n plt.figure(figsize=(30,8))\n plt.plot(x_dn[start:start+slice_size], 'g', alpha=0.5)\n plt.scatter(true_peaks, x_dn[start+true_peaks], marker=\"x\", color=\"b\", label=\"True Peaks\")\n plt.scatter(false_pv, x_dn[start+false_pv], marker=\"x\", color=\"red\", label=\"Cancelled Peaks\")\n plt.vlines(x=true_peaks, ymin=x_dn[start+true_peaks] - peak_heights, ymax = x_dn[start+true_peaks], color = \"C1\")\n plt.hlines(y=peak_width_heights, xmin=peak_left_ips, xmax=peak_right_ips, color = \"C1\")\n ## Valleys\n plt.scatter(true_valleys, x_dn[start+true_valleys], marker=\"x\", color=\"b\", label=\"True Valleys\")\n plt.vlines(x=true_valleys, ymin=x_dn[start+true_valleys] , ymax = x_dn[start+true_valleys] + valley_heights, color = \"C1\")\n plt.hlines(y=-valley_width_heights, xmin=valley_left_ips, xmax=valley_right_ips, color = \"C1\")\n \n plt.title(f\"measurement id: {measurement}, signal id: {signal_id}, phase: {phase}, target: {target}\", fontsize=font_size)\n plt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1854dfcd4b242a80044087f4493e8a53fbd7acce"},"cell_type":"markdown","source":"# Load Train Dataset\n**Note: dataset occupies 6.5GB of RAM, but ~15GB is needed while loading the dataset** \nIf you do not have this much RAM available, load and process the data in chunks as I have done with the larger test set below"},{"metadata":{"trusted":true,"_uuid":"e6e51c117b70d4a7a8c3dba730d1676a542adf60"},"cell_type":"code","source":"train = pq.read_pandas(data_dir + '/train.parquet').to_pandas()\ntrain.info()","execution_count":8,"outputs":[{"output_type":"stream","text":"\nRangeIndex: 800000 entries, 0 to 799999\nColumns: 8712 entries, 0 to 8711\ndtypes: int8(8712)\nmemory usage: 6.5 GB\n","name":"stdout"}]},{"metadata":{"_uuid":"e91404e186199b301121df0a138c8964847a2860"},"cell_type":"markdown","source":"# Train Set: Filter, Denoise, and Extract Features\n**Note: this can take over 30 minutes**\n"},{"metadata":{"trusted":true,"_uuid":"797f01f817d652a1f6a8520240331c0939b4bda5"},"cell_type":"code","source":"%%time\ntrain_length = train.shape[1] # 8712\ntrain_features = np.empty((train_length, 44))\nheight = (2, 100)\nrel_height = 0.2\n\nfor i in range(train_length):\n signal_id = str(i)\n \n # Synchronize phases\n x_sync, dominant_frequency = sync_phase(train[signal_id])\n \n # Apply high pass filter with low cutoff of 10kHz\n x_hp = high_pass_filter(x_sync, low_cutoff=10000, sample_rate=sample_rate)\n \n # Apply denoising\n x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n \n # Peak & valley features, divided into 4 parts of the waveform\n for j in range(4):\n subset_x_dn = x_dn[j*200000:j*200000+200000]\n p, _ = find_peaks(subset_x_dn, height=height, rel_height=rel_height)\n v, _ = find_peaks(-subset_x_dn, height=height, rel_height=rel_height)\n pv = np.sort(np.concatenate((p, v)))\n # If there are no peaks, or the dominant frequancy is not 50Hz, return zeros\n if pv.shape[0] == 0 or dominant_frequency != 50:\n train_features[i,j*11:j*11+11] = np.zeros(11) \n else:\n # Cancel false peaks\n pv_true, _ = cancel_false_peaks(subset_x_dn, pv, min_height_fp=15)\n # peak and valley features\n train_features[i,j*11] = pv_true.shape[0]\n train_features[i,j*11+1:j*11+11] = pv_features(subset_x_dn, pv_true, p, v, rel_height=rel_height)\n\n# Add feature/column names\n# \"p\" -> peaks \n# \"v\" -> valleys\n# \"1\",\"2\",\"3\",\"4\" -> four parts/sections of the waveform\n# \"n\" -> count\n# \"h\" -> height\n# \"w\" -> width\n# e.g. \"pv1_n\" -> peak and valley count in the first section, \"vh3_mean\" -> mean valley heigh in the third section\nfeature_columns = [\"pv1_n\", \"p1_n\", \"ph1_mean\", \"ph1_max\", \"pw1_mean\", \"pw1_max\", \"v1_n\", \"vh1_mean\", \"vh1_max\", \"vw1_mean\", \"vw1_max\", \n \"pv2_n\", \"p2_n\", \"ph2_mean\", \"ph2_max\", \"pw2_mean\", \"pw2_max\", \"v2_n\", \"vh2_mean\", \"vh2_max\", \"vw2_mean\", \"vw2_max\", \n \"pv3_n\", \"p3_n\", \"ph3_mean\", \"ph3_max\", \"pw3_mean\", \"pw3_max\", \"v3_n\", \"vh3_mean\", \"vh3_max\", \"vw3_mean\", \"vw3_max\", \n \"pv4_n\", \"p4_n\", \"ph4_mean\", \"ph4_max\", \"pw4_mean\", \"pw4_max\", \"v4_n\", \"vh4_mean\", \"vh4_max\", \"vw4_mean\", \"vw4_max\",\n ]\n# Combine peak/vally features with names\ntrain_features = pd.DataFrame(train_features, columns=feature_columns)\n# Add signal_id and id_measurement columns\ntrain_features.insert(loc=0, column=\"signal_id\", value=metadata_train.signal_id)\ntrain_features.insert(loc=1, column=\"id_measurement\", value=metadata_train.id_measurement)\n# Save to CSV\ntrain_features.to_csv('train_features.csv', index=False)","execution_count":10,"outputs":[{"output_type":"stream","text":"CPU times: user 32min 39s, sys: 2min 25s, total: 35min 5s\nWall time: 35min 5s\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"### Remove Train from Memory"},{"metadata":{"trusted":true,"_uuid":"12f1924ea4f9096422a7681d8281cb541bfd10f4"},"cell_type":"code","source":"del train","execution_count":11,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot Correlation Matrix"},{"metadata":{"trusted":true,"_uuid":"171dbb245dc3593c7c2dc693f0f18c8b8aecbca4"},"cell_type":"code","source":"features_corr = train_features[[\"p1_n\", \"ph1_mean\", \"ph1_max\", \"pw1_mean\", \"pw1_max\", \"v1_n\", \"vh1_mean\", \"vh1_max\", \"vw1_mean\", \"vw1_max\"]].copy()\n\nfeatures_corr.insert(loc=0, column=\"target\", value=metadata_train.target)\n\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\ncor = features_corr.corr()\nf, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\nsns.heatmap(cor, cmap=blues)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0d9c5ccbc2194e49c9a0dc585d914694841aa013"},"cell_type":"markdown","source":"## Plot Histograms of Features"},{"metadata":{"trusted":true,"_uuid":"92b835a9039a660dd487354c314c1a7e067b6a9f"},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(30, 10))\n \nax[0].hist(train_features[\"p1_n\"], bins='auto')\nax[0].set_title(\"Histogram of Peak Counts in 1st Section\", fontsize=16)\nax[1].hist(train_features[\"ph1_mean\"], bins='auto')\nax[1].set_title(\"Histogram of Mean Peak Heights in 1st Section\", fontsize=16)\nax[2].hist(train_features[\"ph1_max\"], bins='auto')\nax[2].set_title(\"Histogram of Max Peak Heights in 1st Section\", fontsize=16)","execution_count":16,"outputs":[{"output_type":"execute_result","execution_count":16,"data":{"text/plain":"Text(0.5, 1.0, 'Histogram of Max Peak Heights in 1st Section')"},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"
"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"# Test Set: Load, Filter, Denoise, and Extract Features \n**Note: this is done in 10 parts in order to not exceed available RAM** "},{"metadata":{"trusted":true,"_uuid":"243c312f0ed6e97c348d4e83f18aef10ec54a6cf"},"cell_type":"code","source":"%%time\n\nfirst_index = metadata_test.signal_id[0]\nlength_test = metadata_test.shape[0] #20337\nn_parts = 10\nstep_size = int(length_test/n_parts)\nlast_step = length_test % n_parts\n\n# Create a list of lists with start index and end index for each of the 10 parts and one for the last partial part\nstart_end = [[x, x+step_size] for x in range(first_index, length_test + first_index, step_size)]\nstart_end = start_end[:-1] + [[start_end[-1][0], start_end[-1][0] + last_step]]\n# print(start_end)\n\npeak_height = (2,100)\nrel_height = 0.2\ntest_features = np.empty((length_test, 44))\n\nfor start, end in start_end:\n \n subset_test = pq.read_pandas(data_dir + '/test.parquet', columns=[str(x) for x in range(start, end)]).to_pandas()\n \n for i in range(start, end):\n signal_id = str(i)\n \n # Synchronize phases\n x_sync, dominant_frequency = sync_phase(subset_test[signal_id])\n \n # Apply high pass filter with low cutoff of 10kHz\n x_hp = high_pass_filter(x_sync, low_cutoff=10000, sample_rate=sample_rate)\n \n # Apply denoising\n x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n \n # Peak & valley features, divided into 4 parts of the waveform\n for j in range(4):\n subset_x_dn = x_dn[j*200000:j*200000+200000]\n p, _ = find_peaks(subset_x_dn, height=peak_height, rel_height=rel_height)\n v, _ = find_peaks(-subset_x_dn, height=peak_height, rel_height=rel_height)\n pv = np.sort(np.concatenate((p, v)))\n # If there are no peaks, or the dominant frequancy is not 50Hz, return zeros\n if pv.shape[0] == 0 or dominant_frequency != 50:\n test_features[i-first_index, j*11:j*11+11] = np.zeros(11)\n else:\n # Cancel false peaks\n pv_true, _ = cancel_false_peaks(subset_x_dn, pv, min_height_fp=15)\n # Peak and valley features\n test_features[i-first_index, j*11] = pv_true.shape[0]\n test_features[i-first_index, j*11+1:j*11+11] = pv_features(subset_x_dn, pv_true, p, v, rel_height=rel_height)\n\n# Combine peak/vally features with names\ntest_features = pd.DataFrame(test_features, columns=feature_columns)\n# Add signal_id and id_measurement columns\ntest_features.insert(loc=0, column=\"signal_id\", value=metadata_test.signal_id)\ntest_features.insert(loc=1, column=\"id_measurement\", value=metadata_test.id_measurement)\n# Save features to CSV\ntest_features.to_csv('test_features.csv', index=False)\n","execution_count":18,"outputs":[{"output_type":"stream","text":"CPU times: user 1h 24min 36s, sys: 2min 36s, total: 1h 27min 12s\nWall time: 1h 27min 13s\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"# Fit Models\n**First Model: Random Forest** \n**Second Model: Random Forest with all three phases grouped together** \nThis is because about 80% of faults occur on all three phases in the train data, but this percentage was much lower for test set predictions from the first model \n**Stacked Ensemble Model: Use Logistic Regression to combine the other two models\n"},{"metadata":{},"cell_type":"markdown","source":"## Random Forest Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n# Define features to use (all features at this time)\nmodeled_features = [\"p1_n\", \"ph1_mean\", \"ph1_max\", \"pw1_mean\", \"pw1_max\", \"v1_n\", \"vh1_mean\", \"vh1_max\", \"vw1_mean\", \"vw1_max\",\n \"p2_n\", \"ph2_mean\", \"ph2_max\", \"pw2_mean\", \"pw2_max\", \"v2_n\", \"vh2_mean\", \"vh2_max\", \"vw2_mean\", \"vw2_max\",\n \"p3_n\", \"ph3_mean\", \"ph3_max\", \"pw3_mean\", \"pw3_max\", \"v3_n\", \"vh3_mean\", \"vh3_max\", \"vw3_mean\", \"vw3_max\",\n \"p4_n\", \"ph4_mean\", \"ph4_max\", \"pw4_mean\", \"pw4_max\", \"v4_n\", \"vh4_mean\", \"vh4_max\", \"vw4_mean\", \"vw4_max\"\n ]\nX_train = train_features[modeled_features]\n# Define target\ny_train = metadata_train.target\n\n# Random Forest Classifier\nrandom_state = 1\nclass_weight = dict({0:0.5, 1:2.0})\nRF_classifier = RandomForestClassifier(bootstrap=True, class_weight=class_weight, criterion='gini',\n max_depth=8, max_features='auto', max_leaf_nodes=None,\n min_impurity_decrease=0.0, min_impurity_split=None,\n min_samples_leaf=4, min_samples_split=10,\n min_weight_fraction_leaf=0.0, n_estimators=300, n_jobs=-1,\n oob_score=False, random_state=random_state, verbose=0, warm_start=False)\nRF_classifier.fit(X_train, y_train)\n\n# k-fold cross validation\nmcc_scorer = make_scorer(matthews_corrcoef)\nmcc = cross_val_score(estimator = RF_classifier,\n X = X_train,\n y = y_train,\n scoring=mcc_scorer,\n cv = 10)\nprint(\"MCC Mean = \", mcc.mean())\nprint(\"MCC SD = \", mcc.std())\n\n## Predict Test Set Target \nX_test = test_features[modeled_features]\ny_pred = RF_classifier.predict(X_test)","execution_count":27,"outputs":[{"output_type":"stream","text":"MCC Mean = 0.6014003698935049\nMCC SD = 0.08981433705031791\nCPU times: user 15.9 s, sys: 1.27 s, total: 17.1 s\nWall time: 31.9 s\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"## Grouped Random Forest: Three Phases Grouped Together"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Create grouped train features - take the mean of the three phases\ngroup_train_features = train_features.groupby([\"id_measurement\"]).mean()\ngroup_train_features.shape\ngroup_X_train = group_train_features[modeled_features]\n# Create grouped test features - take the mean of the three phases\ngroup_test_features = test_features.groupby([\"id_measurement\"]).mean()\ngroup_test_features.shape\n# Define target\ngroup_y_train = metadata_train.groupby([\"id_measurement\"]).median()[\"target\"] # take median target value\n\n# Grouped Random Forest Classifier\nrandom_state = 1\nclass_weight = dict({0:0.5, 1:2.0})\ngroup_RF_classifier = RandomForestClassifier(bootstrap=True, class_weight=class_weight, criterion='gini',\n max_depth=8, max_features='auto', max_leaf_nodes=None,\n min_impurity_decrease=0.0, min_impurity_split=None,\n min_samples_leaf=4, min_samples_split=10,\n min_weight_fraction_leaf=0.0, n_estimators=300, n_jobs=-1,\n oob_score=False, random_state=random_state, verbose=0, warm_start=False)\ngroup_RF_classifier.fit(group_X_train, group_y_train)\n\n# k-fold cross validation\ngroup_mcc = cross_val_score(estimator = group_RF_classifier,\n X = group_X_train,\n y = group_y_train,\n scoring=mcc_scorer,\n cv = 10)\nprint(\"MCC Mean = \", group_mcc.mean())\nprint(\"MCC SD = \", group_mcc.std())\n\n# Grouped Test Set Predictions\ngroup_X_test = group_test_features[modeled_features]\ngroup_y_pred = group_RF_classifier.predict(group_X_test)\ngroup_y_pred = group_y_pred.repeat(3)\n","execution_count":28,"outputs":[{"output_type":"stream","text":"MCC Mean = 0.6315854420698475\nMCC SD = 0.09014647661276673\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"# Stacked Ensemble (Meta) Model \n**Logistic Regression**"},{"metadata":{"trusted":true,"_uuid":"270093edf1f85dc4ee6f1407cd10a8732a323ab8"},"cell_type":"code","source":"# Use first two RF models to make predictions on the training sets (only made predictions on the test sets earlier)\nmeta_X_train1 = RF_classifier.predict(X_train)\nmeta_X_train2 = np.repeat(group_RF_classifier.predict(group_X_train), 3)\nmeta_X_train = np.column_stack([meta_X_train1, meta_X_train2])\n# Define target\ny_train = metadata_train.target\n\n# Logistic Regression (Classifier)\nlog_classifier = LogisticRegression(random_state = 0)\nlog_classifier.fit(meta_X_train, y_train)\n\n# k-fold cross validation\nmeta_mcc = cross_val_score(estimator = log_classifier,\n X = meta_X_train,\n y = y_train,\n scoring=mcc_scorer,\n cv = 10)\nprint(\"MCC Mean = \", meta_mcc.mean())\nprint(\"MCC SD = \", meta_mcc.std())\n\n# Predict Test Set Target\nmeta_X_test = np.column_stack([y_pred, group_y_pred])\nmeta_y_pred = log_classifier.predict(meta_X_test)\nmeta_y_pred.mean()\n\n# Prepare Kaggle submission\noutput = pd.read_csv(data_dir + '/sample_submission.csv')\noutput['target'] = meta_y_pred\noutput.to_csv('submission3.csv', index=False)\n","execution_count":30,"outputs":[{"output_type":"stream","text":"MCC Mean = 0.8083652690425829\nMCC SD = 0.06764916230251698\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"**My final submission had an MCC of 0.61178**"}],"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}