{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# VSB Power Line Fault Detection\nMedium voltage overhead power lines run for hundreds of miles to supply power to cities. These great distances make it expensive to manually inspect the lines for damage that doesn't immediately lead to a power outage, such as a tree branch hitting the line or a flaw in the insulator. These modes of damage lead to a phenomenon known as partial discharge — an electrical discharge which does not bridge the electrodes between an insulation system completely. Partial discharges slowly damage the power line, so left unrepaired they will eventually lead to a power outage or start a fire.\n\nYour challenge is to detect partial discharge patterns in signals acquired from these power lines with a new meter designed at the ENET Centre at VŠB. Effective classifiers using this data will make it possible to continuously monitor power lines for faults.\n\nENET Centre researches and develops renewable energy resources with the goal of reducing or eliminating harmful environmental impacts. Their efforts focus on developing technology solutions around transportation and processing of energy raw materials.\n\nBy developing a solution to detect partial discharge you’ll help reduce maintenance costs, and prevent power outages.\n\n> The full git repo for my approach to this project is at: https://github.com/jeffreyegan/VSB_Power_Line_Fault_Detection\n\n## Data Exploration and Preparation\n\n### Data Description\nFaults in electric transmission lines can lead to a destructive phenomenon called partial discharge. If left alone, partial discharges can damage equipment to the point that it stops functioning entirely. Your challenge is to detect partial discharges so that repairs can be made before any lasting harm occurs.\n\nEach signal contains 800,000 measurements of a power line's voltage, taken over 20 milliseconds. As the underlying electric grid operates at 50 Hz, this means each signal covers a single complete grid cycle. The grid itself operates on a 3-phase power scheme, and all three phases are measured simultaneously. \n\n### File Descriptions\n**metadata_[train/test].csv**\n*         id_measurement: the ID code for a trio of signals recorded at the same time.\n*         signal_id: the foreign key for the signal data. Each signal ID is unique across both train and test, so the first ID in train is '0' but the first ID in test is '8712'.\n*         phase: the phase ID code within the signal trio. The phases may or may not all be impacted by a fault on the line.\n*         target: 0 if the power line is undamaged, 1 if there is a fault.\n\n**[train/test].parquet** : The signal data. Each column contains one signal; 800,000 int8 measurements as exported with pyarrow.parquet version 0.11. Please note that this is different than our usual data orientation of one row per observation; the switch makes it possible loading a subset of the signals efficiently. If you haven't worked with Apache Parquet before, please refer to either the Python data loading starter kernel.\n \n**sample_submission.csv**: a valid sample submission.\n    \n### Data Background Information\nAs you may notice, the signal data comes from the real environment, not a lab, and they contain a lot of background noise. These signals are measured by our patented device with lower sampling rate (cost efficiency purpose) therefore I do not recommend to use any other publicly available dataset containing partial discharge patterns (PD patterns). - also we deployed the metering devices on more than 20 different locations. This implies that the spectrum of noise and quality of PD's are so different from each other, that the correct and robust classification is a still ongoing problem (the main motivation of this competition). The comparison and broadening of our view is also considered as beneficial and necessary in our research.\n\nPD pattern is therefore the main star tonight and there is a lot of literature about this phenomenon. I would recommend to read some papers, during my dissertation I tried a lot of different feature extraction models, but those based on fundamentals worked the best. In general the imbalanced dataset is very natural because PD pattern implies some degradation or damage of the observed system which is happening (fortunately) less often than the states when the system is operating correctly.\n\nIn our case, the measurements on the medium voltage overhead lines, PD pattern may look like this (pd_pattern.png). But because of a lot of various noise interference (overhead lines work as a huge antena grabbing all signals around), a lot of interpolated patterns may look similar (see samples.png).\n\nTo use any kind of wavelet transformation is very reasonable, butterworth filter was helpful for me to suppress the sine shape, DWT to obtain its close approximation - sometimes it is disrupted, and denoising with feature extractions are the alchemy of this competition.\n    \n### Load Data","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","jupyter":{"source_hidden":true}}},{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport pyarrow.parquet as pq\nfrom matplotlib import pyplot as plt\nimport seaborn as sns","metadata":{"_uuid":"ae82a48aceefec516c0ab896be7b293dca8c3ca2","execution":{"iopub.status.busy":"2022-05-17T13:46:49.535873Z","iopub.execute_input":"2022-05-17T13:46:49.536678Z","iopub.status.idle":"2022-05-17T13:46:49.544150Z","shell.execute_reply.started":"2022-05-17T13:46:49.536601Z","shell.execute_reply":"2022-05-17T13:46:49.542923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subset_train = pq.read_pandas(\"../input/vsb-power-line-fault-detection/train.parquet\", columns=[str(i) for i in range(9)]).to_pandas()","metadata":{"_uuid":"05c509c83859ba03d6f378a4007694a1bf4d15da","execution":{"iopub.status.busy":"2022-05-17T13:46:49.545784Z","iopub.execute_input":"2022-05-17T13:46:49.546538Z","iopub.status.idle":"2022-05-17T13:46:50.022085Z","shell.execute_reply.started":"2022-05-17T13:46:49.546471Z","shell.execute_reply":"2022-05-17T13:46:50.020722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subset_train.info()","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:46:50.023531Z","iopub.execute_input":"2022-05-17T13:46:50.023859Z","iopub.status.idle":"2022-05-17T13:46:50.265321Z","shell.execute_reply.started":"2022-05-17T13:46:50.023808Z","shell.execute_reply":"2022-05-17T13:46:50.264212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig=plt.figure(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\nplot_labels = ['Phase_0', 'Phase_1', 'Phase_2']\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\nplt.plot(list(range(len(subset_train))), subset_train[\"0\"], '-', label=plot_labels[0], color=blues[0])\nplt.plot(list(range(len(subset_train))), subset_train[\"1\"], '-', label=plot_labels[1], color=blues[1])\nplt.plot(list(range(len(subset_train))), subset_train[\"2\"], '-', label=plot_labels[2], color=blues[2])\nplt.ylim((-60, 60))\nplt.legend(loc='lower right')\nplt.title('Raw Signal Data without Partial Discharge Fault')\nplt.xlabel('Sample')\nplt.ylabel('Amplitude [bit]')\n\nfig=plt.figure(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\nplot_labels = ['Phase_0', 'Phase_1', 'Phase_2']\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\nplt.plot(list(range(len(subset_train))), subset_train[\"3\"], '-', label=plot_labels[0], color=blues[0])\nplt.plot(list(range(len(subset_train))), subset_train[\"4\"], '-', label=plot_labels[1], color=blues[1])\nplt.plot(list(range(len(subset_train))), subset_train[\"5\"], '-', label=plot_labels[2], color=blues[2])\nplt.ylim((-60, 60))\nplt.legend(loc='lower right')\nplt.title('Raw Signal Data with Partial Discharge Fault Present')\nplt.xlabel('Sample')\nplt.ylabel('Amplitude [bit]')","metadata":{"_uuid":"654e12a9239b6f24388a6b81dc6ad2423b3c455f","execution":{"iopub.status.busy":"2022-05-17T13:46:50.267035Z","iopub.execute_input":"2022-05-17T13:46:50.267319Z","iopub.status.idle":"2022-05-17T13:46:56.363199Z","shell.execute_reply.started":"2022-05-17T13:46:50.267269Z","shell.execute_reply":"2022-05-17T13:46:56.362324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pywt\nfrom statsmodels.robust import mad\nfrom scipy import signal\n \ndef waveletSmooth( x, wavelet=\"db4\", level=1, title=None ):\n    # calculate the wavelet coefficients\n    coeff = pywt.wavedec( x, wavelet, mode=\"per\" )\n    # calculate a threshold\n    sigma = mad( coeff[-level] )\n    # changing this threshold also changes the behavior,\n    # but I have not played with this very much\n    uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) )\n    coeff[1:] = ( pywt.threshold( i, value=uthresh, mode=\"soft\" ) for i in coeff[1:] )\n    # reconstruct the signal using the thresholded coefficients\n    y = pywt.waverec( coeff, wavelet, mode=\"per\" )\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    blues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\n    plt.plot( x, color=\"#66D7EB\", alpha=0.5, label=\"Original Signal\")\n    plt.plot( y, color=\"#51ACC5\", label=\"Transformed Signal\" )\n    plt.ylim((-60, 60))\n    plt.xlabel('Sample')\n    plt.ylabel('Amplitude')\n    plt.legend(loc='lower right')\n    if title:\n        ax.set_title(title)\n    ax.set_xlim((0,len(y)))\n    return y\ntitle0 = 'Discrete Wavelet Transform De-Noised Signal without Fault'\ntitle1 = 'Discrete Wavelet Transform De-Noised Signal with Fault'\nsignal_0 =  waveletSmooth(subset_train[\"0\"], wavelet=\"db4\", level=1, title=title0 )\nsignal_1 =  waveletSmooth(subset_train[\"3\"], wavelet=\"db4\", level=1, title=title1 )","metadata":{"_uuid":"56604a9daffd6790d2f74949c2e24eae4ce45eff","execution":{"iopub.status.busy":"2022-05-17T13:46:56.364615Z","iopub.execute_input":"2022-05-17T13:46:56.365063Z","iopub.status.idle":"2022-05-17T13:46:59.820261Z","shell.execute_reply.started":"2022-05-17T13:46:56.364999Z","shell.execute_reply":"2022-05-17T13:46:59.818860Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.optimize import leastsq\n\ndef fit_sinusoid(signal, title):\n    t = np.linspace(0, 2*np.pi, len(signal))  # data covers one period\n    guess_mean = np.mean(signal)\n    guess_std = 3*np.std(signal)/(2**0.5)/(2**0.5)\n    guess_phase = 0\n    guess_freq = 1\n    guess_amp = 20\n\n    # Define the function to optimize, in this case, we want to minimize the difference\n    # between the actual data and our \"guessed\" parameters\n    optimize_func = lambda x: x[0]*np.sin(x[1]*t+x[2]) + x[3] - signal\n    est_amp, est_freq, est_phase, est_mean = leastsq(optimize_func, [guess_amp, guess_freq, guess_phase, guess_mean])[0]\n\n    # recreate the fitted curve using the optimized parameters\n    signal_fit = est_amp*np.sin(est_freq*t+est_phase) + est_mean\n\n    # recreate the fitted curve using the optimized parameters\n    blues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    plt.plot(t, signal, '.', color=blues[0], label=\"Transformed Signal\")\n    plt.plot(t, signal_fit, label=\"Sinusoidal Fit\", color=blues[2])\n    plt.xlabel(r\"Sample (Scaled to a Samples over 2$\\pi$ Period)\")\n    plt.ylabel(\"Amplitude\")\n    plt.ylim((-60, 60))\n    plt.title(title)\n    plt.legend()\n    plt.show()\n    return signal_fit\n\ntitle0 = \"Sinusoidal Fit to the Discrete Wavelet Transform De-Noised Signal without Fault\"\ntitle1 = \"Sinusoidal Fit to the Discrete Wavelet Transform De-Noised Signal with Fault\"\nsignal_fit0 = fit_sinusoid(signal_0, title0)\nsignal_fit1 = fit_sinusoid(signal_1, title1)","metadata":{"_uuid":"e82d0330a5715ead9bed3e80286159f869c8202a","execution":{"iopub.status.busy":"2022-05-17T13:46:59.821712Z","iopub.execute_input":"2022-05-17T13:46:59.822292Z","iopub.status.idle":"2022-05-17T13:47:09.899863Z","shell.execute_reply.started":"2022-05-17T13:46:59.822202Z","shell.execute_reply":"2022-05-17T13:47:09.898661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def find_pd_probable(signal_fit, condition):\n    first_derivative = np.gradient(signal_fit)\n    return [i for i, elem in enumerate(first_derivative) if condition(elem)]\n\nhigh_prob_region0 = find_pd_probable(signal_fit0, lambda e: e>0)\nhigh_prob_region1 = find_pd_probable(signal_fit1, lambda e: e>0)","metadata":{"_uuid":"f241f002ec068acc0a9610bf39fb90f1da8638f4","execution":{"iopub.status.busy":"2022-05-17T13:47:09.901219Z","iopub.execute_input":"2022-05-17T13:47:09.901510Z","iopub.status.idle":"2022-05-17T13:47:11.278729Z","shell.execute_reply.started":"2022-05-17T13:47:09.901465Z","shell.execute_reply":"2022-05-17T13:47:11.277433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_pd_probable(signal, signal_fit, high_prob_region, title):\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    plt.plot(signal, color=blues[1], label=\"Transformed, De-Noised Signal\")\n    plt.plot(signal_fit, '--', color=blues[3], label=\"Sinusoidal Best Fit\")\n    plt.scatter(high_prob_region, signal_fit[high_prob_region], color=\"orange\", label=\"PD Probable Region\")\n    plt.legend()\n    plt.ylim((-60, 60))\n    plt.ylabel(\"Amplitude\")\n    plt.xlabel(\"Sample\")\n    plt.title(title)\n    return\n\ntitle0 = \"Sinusoidal Fit with PD-Probable Region Highlighted against the Transformed, De-Noised Signal without Fault\"\ntitle1 = \"Sinusoidal Fit with PD-Probable Region Highlighted against the Transformed, De-Noised Signal with Fault\"\nplot_pd_probable(signal_0, signal_fit0, high_prob_region0, title0)\nplot_pd_probable(signal_1, signal_fit1, high_prob_region1, title1)","metadata":{"_uuid":"ffabccad8aed977cbadfc25860d01d8347c00490","execution":{"iopub.status.busy":"2022-05-17T13:47:11.280500Z","iopub.execute_input":"2022-05-17T13:47:11.280890Z","iopub.status.idle":"2022-05-17T13:47:39.807712Z","shell.execute_reply.started":"2022-05-17T13:47:11.280814Z","shell.execute_reply":"2022-05-17T13:47:39.806374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def detrend_signal( signal, high_prob_idx, title=None ):\n    \n    x = np.diff(signal, n=1)\n    if max(high_prob_idx) == len(x):\n        high_prob_idx = high_prob_idx[0:-1]\n        \n    \n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    blues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\n    plt.plot( x, color=blues[1], label=\"De-Trended, Transformed Signal\")\n    plt.scatter(high_prob_idx, x[high_prob_idx], color=\"orange\", label=\"PD Probable Region of De-Trended, Transformed Signal\" )\n    plt.ylim((-60, 60))\n    plt.xlabel('Sample')\n    plt.ylabel('Amplitude')\n    plt.legend(loc='lower right')\n    if title:\n        ax.set_title(title)\n    ax.set_xlim((0,len(x)))\n    return x[high_prob_idx]\n\n\ntitle0 = 'Discrete Wavelet Transform De-Noised and De-Trended Signal without Fault'\ntitle1 = 'Discrete Wavelet Transform De-Noised and De-Trended Signal with Fault'\nsignal_0hp =  detrend_signal(signal_0, high_prob_region0, title0)\nsignal_1hp =  detrend_signal(signal_1, high_prob_region1, title1)","metadata":{"_uuid":"13ee6e7ad0b5c7a8b3875388de73e1d0b5222e4e","execution":{"iopub.status.busy":"2022-05-17T13:47:39.809369Z","iopub.execute_input":"2022-05-17T13:47:39.809757Z","iopub.status.idle":"2022-05-17T13:48:00.263933Z","shell.execute_reply.started":"2022-05-17T13:47:39.809691Z","shell.execute_reply":"2022-05-17T13:48:00.262768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_pd_prob_detrend(signal_hp, title):\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    plt.plot(signal_hp, color=blues[1])\n    plt.ylim((-60, 60))\n    plt.xlabel(\"Sample\")\n    plt.ylabel(\"Amplitude\")\n    plt.title(title)\n    return\n\ntitle0 = \"PD-Probable Region of De-Trended, Transformed Signal without Fault\"\ntitle1 = \"PD-Probable Region of De-Trended, Transformed Signal with Fault\"\nplot_pd_prob_detrend(signal_0hp, title0)\nplot_pd_prob_detrend(signal_1hp, title1)","metadata":{"_uuid":"493ae6928f5de784923fdb8cd1a9aa3377fb5038","execution":{"iopub.status.busy":"2022-05-17T13:48:00.265838Z","iopub.execute_input":"2022-05-17T13:48:00.266232Z","iopub.status.idle":"2022-05-17T13:48:01.907136Z","shell.execute_reply.started":"2022-05-17T13:48:00.266153Z","shell.execute_reply":"2022-05-17T13:48:01.905876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#import peakutils\ndef find_peaks(signal, title):\n    thresh = 0.6  # used for peakutils\n    thresh = 4.0  # used for fixed amplitude\n    min_d = 0\n\n    #peaks = peakutils.indexes(1.0*(signal), thres=thresh, min_dist=min_d)  # where peaks are\n    peaks = np.argwhere(signal > thresh)\n    #valleys = peakutils.indexes(-1.0*(signal), thres=thresh, min_dist=min_d)  # where peaks are\n    valleys = np.array(np.argwhere(signal < -thresh))\n    pois = np.sort(np.concatenate((peaks, valleys)))\n    \n    peak_indexes = []\n    for pk in pois:\n        #peak_indexes.append(pk)\n        peak_indexes.append(pk[0])\n\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    blues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\n    #plt.plot( x.diff(), color=\"#66D7EB\", alpha=0.5, label=\"De-Trended, Original Signal\")\n    plt.plot( signal, color=\"#51ACC5\", label=\"De-Trended, Transformed Signal\" )\n    plt.scatter(peak_indexes, signal[peak_indexes], marker=\"+\", color=\"red\", label=\"Peaks\" )\n    plt.ylim((-25, 20))\n    plt.ylim((-60, 60))\n    plt.xlim((102500, 102600))\n    plt.xlabel('Sample')\n    plt.ylabel('Amplitude')\n    plt.legend(loc='lower right')\n    if title:\n        ax.set_title(title)\n    ax.set_xlim((0,len(signal)))\n    print(len(peak_indexes))\n    return np.sort(peak_indexes)\n\ntitle0 = 'Detected Peaks in PD-Probable Region of the De-Trended, Transformed, and De-Noised Signal without Fault'\ntitle1 = 'Detected Peaks in PD-Probable Region of the De-Trended, Transformed, and De-Noised Signal with Fault'\npeak_indexes0 = find_peaks(signal_0hp, title=title0 )\npeak_indexes1 = find_peaks(signal_1hp, title=title1 )","metadata":{"_uuid":"70d413c889424da06dc18d73808208b850e84015","execution":{"iopub.status.busy":"2022-05-17T13:48:01.909270Z","iopub.execute_input":"2022-05-17T13:48:01.909705Z","iopub.status.idle":"2022-05-17T13:48:03.487965Z","shell.execute_reply.started":"2022-05-17T13:48:01.909629Z","shell.execute_reply":"2022-05-17T13:48:03.486782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cancel_false_peaks(signal, peak_indexes):\n    false_peak_indexes = []\n    max_sym_distance = 10  #\n    max_pulse_train = 20  # \n    max_height_ratio = 0.75  # \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]] > 0 and signal[peak_indexes[pk+1]] < 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:\n                if min(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]]))/max(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]])) > max_height_ratio:\n                    scrub = list(x for x in range(len(peak_indexes)) if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train)\n                    for x in scrub:\n                        false_peak_indexes.append(peak_indexes[x])\n\n            if (signal[peak_indexes[pk]] < 0 and signal[peak_indexes[pk+1]] > 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:\n                if min(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]]))/max(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]])) > max_height_ratio:\n                    scrub = list(x for x in range(len(peak_indexes)) if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train)\n                    for x in scrub:\n                        false_peak_indexes.append(peak_indexes[x])\n    return false_peak_indexes\n\nfalse_peak_indexes0 = cancel_false_peaks(signal_0hp, peak_indexes0)\nfalse_peak_indexes1 = cancel_false_peaks(signal_1hp, peak_indexes1)","metadata":{"_uuid":"7150f81df2f3902ea403acf932a230f487a56621","execution":{"iopub.status.busy":"2022-05-17T13:48:03.489320Z","iopub.execute_input":"2022-05-17T13:48:03.489915Z","iopub.status.idle":"2022-05-17T13:48:03.549288Z","shell.execute_reply.started":"2022-05-17T13:48:03.489796Z","shell.execute_reply":"2022-05-17T13:48:03.548431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cancel_high_amp_peaks(signal, peak_indexes, false_peak_indexes):\n    thresh = 30.0\n    #peaks = peakutils.indexes(1.0*(signal), thres=0.80, min_dist=0)\n    peaks = np.argwhere(signal > thresh)\n    #valleys = peakutils.indexes(-1.0*(signal), thres=0.80, min_dist=0) \n    valleys = np.argwhere(signal < -thresh)\n    hi_amp_pk_indexes = np.sort(np.concatenate((peaks, valleys)))\n    for pk_idx in hi_amp_pk_indexes:\n        if not pk_idx in false_peak_indexes:\n            false_peak_indexes.append(pk_idx[0])\n    return np.sort(false_peak_indexes)\n\nfalse_peak_indexes0 = cancel_high_amp_peaks(signal_0hp, peak_indexes0, false_peak_indexes0)\nfalse_peak_indexes1 = cancel_high_amp_peaks(signal_1hp, peak_indexes1, false_peak_indexes1)","metadata":{"_uuid":"e765da72990230a991d100da9abd5ac00372a147","execution":{"iopub.status.busy":"2022-05-17T13:48:03.550601Z","iopub.execute_input":"2022-05-17T13:48:03.551050Z","iopub.status.idle":"2022-05-17T13:48:03.579074Z","shell.execute_reply.started":"2022-05-17T13:48:03.550997Z","shell.execute_reply":"2022-05-17T13:48:03.577904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calcel Peaks Flagged as False, Find True Preaks\ndef cancel_flagged_peaks(peak_indexes, false_peak_indexes):\n    true_peak_indexes = set(peak_indexes) - set(false_peak_indexes)\n    return np.sort(list(true_peak_indexes))\n\ntrue_peak_indexes0 = cancel_flagged_peaks(peak_indexes0, false_peak_indexes0)\ntrue_peak_indexes1 = cancel_flagged_peaks(peak_indexes1, false_peak_indexes1)","metadata":{"_uuid":"a17608b9a4a8bf101203180b488b6273bc59395a","execution":{"iopub.status.busy":"2022-05-17T13:48:03.580959Z","iopub.execute_input":"2022-05-17T13:48:03.581278Z","iopub.status.idle":"2022-05-17T13:48:03.588228Z","shell.execute_reply.started":"2022-05-17T13:48:03.581219Z","shell.execute_reply":"2022-05-17T13:48:03.587240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_peaks(signal, true_peak_indexes, false_peak_indexes, title):\n    f, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\n    plt.plot( signal, color=\"#66D7EB\", label=\"De-Trended, Transformed Signal\" )  #5 alt color 1ACC5\n    plt.scatter(false_peak_indexes, signal[false_peak_indexes], marker=\"x\", color=\"red\", label=\"Cancelled Peaks\")\n    plt.scatter(true_peak_indexes, signal[true_peak_indexes], marker=\"+\", color=\"green\", label=\"True Peaks\")\n    plt.ylim((-25, 20))\n    plt.ylim((-60, 60))\n    plt.xlim((102500, 102600))\n    plt.xlabel('Sample')\n    plt.ylabel('Amplitude')\n    plt.legend(loc='lower right')\n    if title:\n        ax.set_title(title)\n    ax.set_xlim((0,len(signal)))\n    plt.show()\n    return\n\ntitle0 = 'Post Peak Processing on PD-Probable Region of De-Trended and De-Noised Signal without Fault'\ntitle1 = 'Post Peak Processing on PD-Probable Region of De-Trended and De-Noised Signal with Fault'\nplot_peaks(signal_0hp, true_peak_indexes0, false_peak_indexes0, title0)\nplot_peaks(signal_1hp, true_peak_indexes1, false_peak_indexes1, title1)","metadata":{"_uuid":"ec17a6a0d5c8487ab3039e06feb81e056bc2093d","execution":{"iopub.status.busy":"2022-05-17T13:48:03.589496Z","iopub.execute_input":"2022-05-17T13:48:03.590047Z","iopub.status.idle":"2022-05-17T13:48:05.238631Z","shell.execute_reply.started":"2022-05-17T13:48:03.590000Z","shell.execute_reply":"2022-05-17T13:48:05.237863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Full feature extraction code (without the per signal looping)","metadata":{"_uuid":"d395655c7c00d980b1b465d98e1841d87e4ef036"}},{"cell_type":"code","source":"import scipy\nfrom statsmodels.robust import mad\nimport collections\n\n#  Feature Extraction Functions\ndef calculate_entropy(signal):\n    counter_values = collections.Counter(signal).most_common()\n    probabilities = [elem[1]/len(signal) for elem in counter_values]\n    entropy=scipy.stats.entropy(probabilities)\n    return entropy\n\ndef calculate_statistics(signal):\n    n5 = np.nanpercentile(signal, 5)\n    n25 = np.nanpercentile(signal, 25)\n    n75 = np.nanpercentile(signal, 75)\n    n95 = np.nanpercentile(signal, 95)\n    median = np.nanpercentile(signal, 50)\n    mean = np.nanmean(signal)\n    std = np.nanstd(signal)\n    var = np.nanvar(signal)\n    rms = np.nanmean(np.sqrt(signal**2))\n    return [n5, n25, n75, n95, median, mean, std, var, rms]\n\ndef calculate_crossings(signal):\n    zero_crossing_indices = np.nonzero(np.diff(np.array(signal) > 0))[0]\n    no_zero_crossings = len(zero_crossing_indices)\n    mean_crossing_indices = np.nonzero(np.diff(np.array(signal) > np.nanmean(signal)))[0]\n    no_mean_crossings = len(mean_crossing_indices)\n    return [no_zero_crossings, no_mean_crossings]\n\ndef find_all_peaks(signal, threshold=0.7, min_distance=0):\n    #peaks = peakutils.indexes(1.0*(signal), thres=threshold, min_dist=min_distance)\n    peaks = np.argwhere(signal > 1.0*threshold)\n    #valleys = peakutils.indexes(-1.0*(signal), thres=threshold, min_dist=min_distance)\n    valleys = np.array(np.argwhere(signal < -1.0*threshold))\n    pois = np.sort(np.concatenate((peaks, valleys)))\n    peak_indexes = []\n    for pk in pois:\n        #peak_indexes.append(pk)\n        peak_indexes.append(pk[0])\n    return np.sort(peak_indexes), peaks, valleys\n\ndef calculate_peak_widths(peak_idxs):\n    tmp_w = 1\n    widths = []\n    for idx in range(1,len(peak_idxs)):\n        if peak_idxs[idx]-peak_idxs[idx-1] < 3:\n            tmp_w +=1\n        else:\n            widths.append(tmp_w)\n            tmp_w = 1\n    widths.append(tmp_w)\n    min_width = min(np.array(widths))\n    max_width = max(np.array(widths))\n    mean_width = np.nanmean(np.array(widths))\n    num_true_peaks = len(widths)\n    return min_width, max_width, mean_width, num_true_peaks\n\ndef cancel_false_peaks(signal, peak_indexes):\n    false_peak_indexes = []\n    max_sym_distance = 10  #\n    max_pulse_train = 500  # \n    max_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]] > 0 and signal[peak_indexes[pk+1]] < 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:  # opposite polarity and within symmetric check distance\n                if min(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]]))/max(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]])) > max_height_ratio:  # ratio of opposing polarity check\n                    scrub = list(x for x in range(len(peak_indexes)) if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train)  # build pulse train\n                    for x in scrub:\n                        false_peak_indexes.append(peak_indexes[x])\n            if (signal[peak_indexes[pk]] < 0 and signal[peak_indexes[pk+1]] > 0) and (peak_indexes[pk+1] - peak_indexes[pk]) < max_sym_distance:\n                if min(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]]))/max(abs(signal[peak_indexes[pk]]),abs(signal[peak_indexes[pk+1]])) > max_height_ratio:\n                    scrub = list(x for x in range(len(peak_indexes)) if peak_indexes[pk] <= peak_indexes[x] <= peak_indexes[pk]+max_pulse_train)\n                    for x in scrub:\n                        false_peak_indexes.append(peak_indexes[x])\n    return false_peak_indexes\n\ndef cancel_high_amp_peaks(signal, peak_indexes, false_peak_indexes):\n    threshold = 50  # amplitude threshld for determining high amplitude peaks for cancellation\n    #peaks = peakutils.indexes(1.0*(signal), thres=0.80, min_dist=0)\n    peaks = np.argwhere(signal > 1.0*threshold)\n    #valleys = peakutils.indexes(-1.0*(signal), thres=0.80, min_dist=0)\n    valleys = np.argwhere(signal < -1.0*threshold)\n    hi_amp_pk_indexes = np.sort(np.concatenate((peaks, valleys)))\n    for pk_idx in hi_amp_pk_indexes:\n        if not pk_idx[0] in false_peak_indexes:\n            false_peak_indexes.append(pk_idx[0])\n    return false_peak_indexes\n\ndef cancel_flagged_peaks(peak_indexes, false_peak_indexes):\n    true_peak_indexes = list(set(peak_indexes) - set(false_peak_indexes))\n    true_peak_indexes.sort()\n    return true_peak_indexes\n\ndef low_high_peaks(signal, true_peak_indexes, hi_idx, lo_idx):\n    if np.size(np.array(hi_idx))> 0:\n        lhr = 1.0*np.size(np.array(lo_idx))/np.size(np.array(hi_idx))\n    else:\n        lhr = 0.0\n    hi_true = 0\n    lo_true = 0\n    for x in true_peak_indexes:\n        if signal[x] > 0.0:\n            hi_true += 1\n        else:\n            lo_true += 1\n    if hi_true > 0:\n        lhrt = 1.0*lo_true/hi_true\n    else:\n        lhrt = 0.0\n    return [hi_idx, lo_idx, lhr, hi_true, lo_true, lhrt]\n\ndef calculate_peaks(signal, true_peak_indexes):  # Peak Characteristics on True Peaks\n    peak_values = signal[true_peak_indexes]\n    num_detect_peak = len(true_peak_indexes)\n    if num_detect_peak > 0:\n        min_height = min(peak_values)\n        max_height = max(peak_values)\n        mean_height = np.nanmean(peak_values)\n        min_width, max_width, mean_width, num_true_peaks = calculate_peak_widths(true_peak_indexes)\n        return [min_height, max_height, mean_height, min_width, max_width, mean_width, num_detect_peak, num_true_peaks]\n    else:\n        return [0, 0, 0, 0, 0, 0, 0, 0]\n\ndef get_features(signal, signal_id, threshold, min_distance): # Extract features from the signal and build an array of them\n    peak_indexes, hi_idx, lo_idx = find_all_peaks(signal, threshold, min_distance)\n    #print(\"Now processing signal_id: \"+str(signal_id)+\" with peak detection threshold at \"+str(threshold)+\" yielding \"+str(len(peak_indexes))+\" peaks at \"+datetime.now().strftime('%Y-%m-%d %H:%M:%S'))\n    false_peak_indexes = cancel_false_peaks(signal, peak_indexes)\n    false_peak_indexes = cancel_high_amp_peaks(signal, peak_indexes, false_peak_indexes)\n    #false_peak_indexes = []  # Don't cancel peaks\n    true_peak_indexes = cancel_flagged_peaks(peak_indexes, false_peak_indexes)\n\n    entropy = calculate_entropy(signal)\n    statistics = calculate_statistics(signal)\n    crossings = calculate_crossings(signal)\n    peaks = calculate_peaks(signal, true_peak_indexes)\n    low_high_stats = low_high_peaks(signal, true_peak_indexes, hi_idx, lo_idx)\n    return [entropy] + statistics + crossings + peaks + low_high_stats\n\nthreshold = 0.45  # Threshold for db_bit fixed threshold, not peak_detector package\nmin_distance = 0  # Minimum distance in samples between detected peaks\n\nsignal_id = 0  # in this example use signal_id = 0\ndata_type = \"train\"  # signal_id 0 belongs to the training set, not the test set in this example\nsignal_features = get_features(signal_0hp, signal_id, threshold, min_distance)","metadata":{"_uuid":"d654849dd6d85510be2926cfb58c89b2d172dc2e","execution":{"iopub.status.busy":"2022-05-17T13:48:05.239909Z","iopub.execute_input":"2022-05-17T13:48:05.240287Z","iopub.status.idle":"2022-05-17T13:48:05.913427Z","shell.execute_reply.started":"2022-05-17T13:48:05.240244Z","shell.execute_reply":"2022-05-17T13:48:05.912538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initialize Feature Dataframe\nfeature_matrix_columns = [\"signal_id\", \"entropy\", \"n5\", \"n25\", \"n75\", \"n95\", \"median\", \"mean\", \"std\", \"var\", \"rms\", \"no_zero_crossings\", \"no_mean_crossings\", \"min_height\", \"max_height\", \"mean_height\", \"min_width\", \"max_width\", \"mean_width\", \"num_detect_peak\", \"num_true_peaks\", \"hi_count\", \"lo_count\", \"low_high_ratio\", \"hi_true\", \"lo_true\", \"low_high_ratio_true\", \"fault\"]\nfeature_matrix = pd.DataFrame([], columns=feature_matrix_columns)\n\n# Store Extracted Features om the Feature MAtrix\n#if data_type.lower() == \"train\":  # Stage Feature Array for Addition to Feature Matrix\n#    df_features = pd.DataFrame([[signal_id] + signal_features + [df_meta.target[df_meta.signal_id == signal_id].values[0]]], columns=feature_matrix_columns)\n#else:  # for test data\n#    df_features = pd.DataFrame([[signal_id] + signal_features + [np.NaN]], columns=feature_matrix_columns)\n#feature_matrix = feature_matrix.append(df_features, ignore_index=True)  # Append Feature Matrix Data Frame\n\n# Store feature matrix to CSV (commented out for this notebook, but left in for example.)\n#feature_matrix.to_csv(\"extracted_features/\"+data_type+\"_features.csv\", sep=\",\")","metadata":{"_uuid":"22b9e084aeb1261b351390fe390c7701b0759e6a","execution":{"iopub.status.busy":"2022-05-17T13:48:05.914644Z","iopub.execute_input":"2022-05-17T13:48:05.915078Z","iopub.status.idle":"2022-05-17T13:48:05.929176Z","shell.execute_reply.started":"2022-05-17T13:48:05.915030Z","shell.execute_reply":"2022-05-17T13:48:05.927752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Explore Extracted Feature Data\nOutside of this notebook, features have been extracted for all signals in the training and test data sets as described above and stored into CSV files. \n\nThis section examines extracted features in the training set, data useful for input into a machine learning model.","metadata":{"_uuid":"f97cd4554266b377e5b5fc92dc152b5fe67decd4"}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\ndf_all = pd.read_csv('../input/threshold/train_featuresHiLo_thresh_4.5_db4.csv') #, sep='|', compression = 'gzip', dtype='unicode', error_bad_lines=False)\nfeatures = [\"entropy\", \"n5\", \"n25\", \"n75\", \"n95\", \"median\", \"mean\", \"std\", \"var\", \"rms\", \"no_zero_crossings\", \"no_mean_crossings\", \"min_height\", \"max_height\", \"mean_height\", \"min_width\", \"max_width\", \"mean_width\", \"num_detect_peak\", \"num_true_peaks\", \"low_high_ratio\", \"hi_true\", \"lo_true\", \"low_high_ratio_true\"]\ntarget = [\"fault\"]\n\ndf = df_all[features+target]\ndf.describe()\ndf[\"fault\"].value_counts()","metadata":{"_uuid":"bd392455ddc1901ef79494f8f6d04158e1c649e6","execution":{"iopub.status.busy":"2022-05-17T13:48:05.930756Z","iopub.execute_input":"2022-05-17T13:48:05.931109Z","iopub.status.idle":"2022-05-17T13:48:06.295281Z","shell.execute_reply.started":"2022-05-17T13:48:05.931045Z","shell.execute_reply":"2022-05-17T13:48:06.294140Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from matplotlib import pyplot as plt\nimport seaborn as sns\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\ncor = df.corr()\nf, ax = plt.subplots(figsize=(14, 8), dpi= 120, facecolor='w', edgecolor='k')\nsns.heatmap(cor, cmap=blues)","metadata":{"_uuid":"3e7eaacf4366d03fc7bda4667eadd389fe0af1ad","execution":{"iopub.status.busy":"2022-05-17T13:48:06.296789Z","iopub.execute_input":"2022-05-17T13:48:06.297117Z","iopub.status.idle":"2022-05-17T13:48:07.660711Z","shell.execute_reply.started":"2022-05-17T13:48:06.297067Z","shell.execute_reply":"2022-05-17T13:48:07.659919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Detail\nsns.pairplot(df[[\"std\", \"max_width\", \"hi_true\", \"fault\"]], hue=\"fault\", palette=blues, height=3, diag_kind=\"kde\", diag_kws=dict(shade=True, bw=.05, vertical=False) )\n\n# Full (takes a while to execute the full pairplot, uncomment out if you want to dig in)\n# sns.pairplot(df, hue=\"fault\", palette=blues, height=3, diag_kind=\"kde\", diag_kws=dict(shade=True, bw=.05, vertical=False) )","metadata":{"_uuid":"a228e85db6f770ead7ddd35380d868e448c92304","execution":{"iopub.status.busy":"2022-05-17T13:48:07.661824Z","iopub.execute_input":"2022-05-17T13:48:07.662189Z","iopub.status.idle":"2022-05-17T13:48:14.431321Z","shell.execute_reply.started":"2022-05-17T13:48:07.662149Z","shell.execute_reply":"2022-05-17T13:48:14.430091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#pd.plotting.scatter_matrix(df, alpha=0.2, figsize=(10, 10), diagonal='kde')","metadata":{"_uuid":"ccdc8b2bbdbe6bcafe2b9d988ce25ad783e7d99c","execution":{"iopub.status.busy":"2022-05-17T13:48:14.433537Z","iopub.execute_input":"2022-05-17T13:48:14.433975Z","iopub.status.idle":"2022-05-17T13:48:14.439967Z","shell.execute_reply.started":"2022-05-17T13:48:14.433908Z","shell.execute_reply":"2022-05-17T13:48:14.438412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Machine Learning Model\nOnce all signals havebeen run through signal processing and features extracted and stored in a new CSV file. Load that CSV file into a data frame.","metadata":{"_uuid":"a9fa9b8db819a362a105f10bdcfec3c0767153c1"}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport warnings\n\n## Plotting libraries\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n## Sklearn Libraries\nfrom sklearn.utils import shuffle\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.preprocessing import label_binarize\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.metrics import make_scorer\nfrom sklearn.metrics import matthews_corrcoef\nfrom sklearn.metrics import f1_score, confusion_matrix, roc_curve, auc, classification_report, recall_score, precision_recall_curve\n\n\nclass Create_ensemble(object):\n    def __init__(self, n_splits, base_models):\n        self.n_splits = n_splits\n        self.base_models = base_models\n\n    def predict(self, X, y, T):\n        X = np.array(X)\n        y = np.array(y)\n        T = np.array(T)\n        no_class = len(np.unique(y))\n\n        folds = list(StratifiedKFold(n_splits=self.n_splits, shuffle=True, \n                                     random_state = random_state).split(X, y))\n\n        train_proba = np.zeros((X.shape[0], no_class))\n        test_proba = np.zeros((T.shape[0], no_class))\n        \n        train_pred = np.zeros((X.shape[0], len(self.base_models)))\n        test_pred = np.zeros((T.shape[0], len(self.base_models)* self.n_splits))\n        f1_scores = np.zeros((len(self.base_models), self.n_splits))\n        recall_scores = np.zeros((len(self.base_models), self.n_splits))\n        \n        test_col = 0\n        for i, clf in enumerate(self.base_models):\n            \n            for j, (train_idx, valid_idx) in enumerate(folds):\n                \n                X_train = X[train_idx]\n                Y_train = y[train_idx]\n                X_valid = X[valid_idx]\n                Y_valid = y[valid_idx]\n                \n                clf.fit(X_train, Y_train)\n                \n                valid_pred = clf.predict(X_valid)\n                recall  = recall_score(Y_valid, valid_pred, average='macro')\n                f1 = f1_score(Y_valid, valid_pred, average='macro')\n                mcc = matthews_corrcoef(Y_valid, valid_pred)\n                \n                recall_scores[i][j] = recall\n                f1_scores[i][j] = f1\n                \n                train_pred[valid_idx, i] = valid_pred\n                test_pred[:, test_col] = clf.predict(T)\n                test_col += 1\n                \n                ## Probabilities\n                valid_proba = clf.predict_proba(X_valid)\n                train_proba[valid_idx, :] = valid_proba\n                test_proba  += clf.predict_proba(T)\n                \n                print( \"Model- {} and CV- {} recall: {}, f1_score: {}, mcc: {}\".format(i, j, recall, f1, mcc))\n                \n            test_proba /= self.n_splits\n            \n        return train_proba, test_proba, train_pred, test_pred, clf\n\ndef re_predict(data, threshods):\n    argmax = np.argmax(data)\n    if argmax == 1:\n        if data[argmax] >= threshods[argmax] : \n            return 1\n        else:\n            return 0\n    else:  # argmax == 0 \n        if data[argmax] >= threshods[argmax] : \n            return 0\n        else:\n            return 1","metadata":{"_uuid":"dae18bfb73d710293b3cde3078d15f9889dfb3c4","execution":{"iopub.status.busy":"2022-05-17T13:48:14.443611Z","iopub.execute_input":"2022-05-17T13:48:14.443978Z","iopub.status.idle":"2022-05-17T13:48:14.465383Z","shell.execute_reply.started":"2022-05-17T13:48:14.443909Z","shell.execute_reply":"2022-05-17T13:48:14.464033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Data\nfeatures = [\"entropy\", \"n5\", \"n25\", \"n75\", \"n95\", \"median\", \"mean\", \"std\", \"var\", \"rms\", \"no_zero_crossings\", \"no_mean_crossings\", \"min_height\", \"max_height\", \"mean_height\", \"min_width\", \"max_width\", \"mean_width\", \"num_detect_peak\", \"num_true_peaks\", \"low_high_ratio\", \"hi_true\", \"lo_true\", \"low_high_ratio_true\"]\ntarget = [\"fault\"]\n\ndata_file = \"../input/threshold/train_featuresHiLo_thresh_4.5_db4.csv\"\ndf_train = pd.read_csv(data_file)\ntrain = df_train[features + target]\n\ndata_file = \"../input/thresholdtest/test_featuresHiLo_thresh_4.5_db4.csv\"\ndf_test = pd.read_csv(data_file)\ntest = df_test[features]","metadata":{"_uuid":"955b27515dd15319a7bf2f849f4753ee3de192bb","execution":{"iopub.status.busy":"2022-05-17T13:48:14.467049Z","iopub.execute_input":"2022-05-17T13:48:14.467362Z","iopub.status.idle":"2022-05-17T13:48:15.373298Z","shell.execute_reply.started":"2022-05-17T13:48:14.467313Z","shell.execute_reply":"2022-05-17T13:48:15.371864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[\"fault\"].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:15.375048Z","iopub.execute_input":"2022-05-17T13:48:15.375778Z","iopub.status.idle":"2022-05-17T13:48:15.386481Z","shell.execute_reply.started":"2022-05-17T13:48:15.375452Z","shell.execute_reply":"2022-05-17T13:48:15.385433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from sklearn.utils import class_weight\n\n# class_weights = class_weight.compute_class_weight('balanced',np.unique(train['fault'].values),train['fault'].values)\n# print(class_weights)","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:15.388408Z","iopub.execute_input":"2022-05-17T13:48:15.389125Z","iopub.status.idle":"2022-05-17T13:48:15.397502Z","shell.execute_reply.started":"2022-05-17T13:48:15.388731Z","shell.execute_reply":"2022-05-17T13:48:15.395991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"random_state = 1\nclass_weight = dict({0:0.5, 1:2.0})\n\nrdf = 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,\n        random_state=random_state, verbose=0, warm_start=False)\n\nbase_models = [rdf]\nn_splits = 5\nlgb_stack = Create_ensemble(n_splits = n_splits, base_models = base_models)        \n\nxtrain = train.drop(['fault'], axis=1)\nytrain = train['fault'].values\n\ntrain_proba, test_proba, train_pred, test_pred, clf = lgb_stack.predict(xtrain, ytrain, test)","metadata":{"_uuid":"140d7e260f3da464db65ee5911f5bdedec1881a6","execution":{"iopub.status.busy":"2022-05-17T13:48:15.400049Z","iopub.execute_input":"2022-05-17T13:48:15.400824Z","iopub.status.idle":"2022-05-17T13:48:28.227377Z","shell.execute_reply.started":"2022-05-17T13:48:15.400753Z","shell.execute_reply":"2022-05-17T13:48:28.226470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('\\nPerformance Metrics after Weighted Random Forest Cross Validation')\nprint('1. The F-1 score of the model {}\\n'.format(f1_score(ytrain, train_pred, average='macro')))\nprint('2. The recall score of the model {}\\n'.format(recall_score(ytrain, train_pred, average='macro')))\nprint('3. The Matthews Correlation Coefficient: {}\\n'.format(matthews_corrcoef(ytrain, train_pred)))\nprint('4. Classification report \\n {} \\n'.format(classification_report(ytrain, train_pred)))\nprint('5. Confusion matrix \\n {} \\n'.format(confusion_matrix(ytrain, train_pred)))","metadata":{"_uuid":"8ef8120993dd99eceb62717d1aaa7b96d43f9e05","execution":{"iopub.status.busy":"2022-05-17T13:48:28.229026Z","iopub.execute_input":"2022-05-17T13:48:28.229294Z","iopub.status.idle":"2022-05-17T13:48:28.319048Z","shell.execute_reply.started":"2022-05-17T13:48:28.229249Z","shell.execute_reply":"2022-05-17T13:48:28.317841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# histogram of important features\nimp = clf.feature_importances_\nimp, features = zip(*sorted(zip(imp, features)))\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\nf, ax = plt.subplots(figsize=(14, 6), dpi= 120, facecolor='w', edgecolor='k')\nplt.barh(range(len(features)), imp, color=blues[1], align=\"center\")\nplt.yticks(range(len(features)), features)\nplt.xlabel(\"Importance of Features\")\nplt.ylabel(\"Features\")\nplt.title(\"Importance of Each Feature in Classifier Model\")\nplt.show()","metadata":{"_uuid":"5f5655f3728e7f9c3f8ea0efaaa3b0c24a9d7819","execution":{"iopub.status.busy":"2022-05-17T13:48:28.321040Z","iopub.execute_input":"2022-05-17T13:48:28.321618Z","iopub.status.idle":"2022-05-17T13:48:29.058520Z","shell.execute_reply.started":"2022-05-17T13:48:28.321388Z","shell.execute_reply":"2022-05-17T13:48:29.057172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# histogram of predicted probabilities\nblues = [\"#66D7EB\", \"#51ACC5\", \"#3E849E\", \"#2C5F78\", \"#1C3D52\", \"#0E1E2B\"]\nf, ax = plt.subplots(figsize=(14, 6), dpi= 120, facecolor='w', edgecolor='k')\nnclasses = 2\ntitles = [\"Probabilities for No Partial Discharge Fault Present\", \"Probabilities for Partial Discharge Fault Present\"]\nfor i in range(nclasses):\n    plt.subplot(1, nclasses, i+1)\n    plt.hist(train_proba[:, i], bins=50, histtype='bar', rwidth=0.95, color=blues[1])\n    plt.xlim(0,1)\n    plt.title(titles[i])\n    plt.xlabel('Probability')\n    plt.ylabel('Frequency')\nplt.tight_layout()\nplt.show()","metadata":{"_uuid":"82754c171277c0b04be8c414a1eefc55ee88d815","execution":{"iopub.status.busy":"2022-05-17T13:48:29.060396Z","iopub.execute_input":"2022-05-17T13:48:29.060808Z","iopub.status.idle":"2022-05-17T13:48:30.754360Z","shell.execute_reply.started":"2022-05-17T13:48:29.060736Z","shell.execute_reply":"2022-05-17T13:48:30.753434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = label_binarize(ytrain, classes=[0, 1])\n_, _, th1 = roc_curve(y[:, 0], train_proba[:, 0])\n_, _, th2 = roc_curve(y[:, 0], train_proba[:, 1])\nprint('\\nMedian Detection Thresholds for Fault Detection')  # use for setting reprediction thresholds\nprint(np.median(th1))\nprint(np.median(th2))","metadata":{"_uuid":"492c56cb7bd6ff1efc0e0d9d495f57f2f53a42e4","execution":{"iopub.status.busy":"2022-05-17T13:48:30.755940Z","iopub.execute_input":"2022-05-17T13:48:30.756252Z","iopub.status.idle":"2022-05-17T13:48:30.770798Z","shell.execute_reply.started":"2022-05-17T13:48:30.756196Z","shell.execute_reply":"2022-05-17T13:48:30.769853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from plotnine import *\nimport plotnine\n\n# Create the ROC curve\nfpr, tpr, thresholds = roc_curve(y[:, 0], train_proba[:, 1])\n\n# Plot the ROC curve\ndf_fpr_tpr = pd.DataFrame({'FPR':fpr, 'TPR':tpr, 'Threshold':thresholds})\ndf_fpr_tpr.head()\n\n# Create the data viz\nplotnine.options.figure_size = (8, 4.8)\n(\n    ggplot(data = df_fpr_tpr)+\n    geom_point(aes(x = 'FPR',\n                   y = 'TPR'),\n               size = 0.4)+\n    geom_line(aes(x = 'FPR',\n                  y = 'TPR'))+\n    labs(title = 'ROC Curve')+\n    xlab('False Positive Rate')+\n    ylab('True Positive Rate')+\n    theme_minimal()\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:30.772341Z","iopub.execute_input":"2022-05-17T13:48:30.772757Z","iopub.status.idle":"2022-05-17T13:48:31.462031Z","shell.execute_reply.started":"2022-05-17T13:48:30.772672Z","shell.execute_reply":"2022-05-17T13:48:31.461156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from plotnine import *\nimport plotnine\n\n# Create the ROC curve\nfpr, tpr, thresholds = roc_curve(y[:, 0], train_proba[:, 0])\n\n# Plot the ROC curve\ndf_fpr_tpr = pd.DataFrame({'FPR':fpr, 'TPR':tpr, 'Threshold':thresholds})\ndf_fpr_tpr.head()\n\n# Create the data viz\nplotnine.options.figure_size = (8, 4.8)\n(\n    ggplot(data = df_fpr_tpr)+\n    geom_point(aes(x = 'FPR',\n                   y = 'TPR'),\n               size = 0.4)+\n    geom_line(aes(x = 'FPR',\n                  y = 'TPR'))+\n    labs(title = 'ROC Curve')+\n    xlab('False Positive Rate')+\n    ylab('True Positive Rate')+\n    theme_minimal()\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:31.463793Z","iopub.execute_input":"2022-05-17T13:48:31.464060Z","iopub.status.idle":"2022-05-17T13:48:32.153163Z","shell.execute_reply.started":"2022-05-17T13:48:31.464016Z","shell.execute_reply":"2022-05-17T13:48:32.152142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the G-mean\ngmean = np.sqrt(tpr * (1 - fpr))\n\n# Find the optimal threshold\nindex = np.argmax(gmean)\nthresholdOpt = round(thresholds[index], ndigits = 4)\ngmeanOpt = round(gmean[index], ndigits = 4)\nfprOpt = round(fpr[index], ndigits = 4)\ntprOpt = round(tpr[index], ndigits = 4)\nprint('Best Threshold: {} with G-Mean: {}'.format(thresholdOpt, gmeanOpt))\nprint('FPR: {}, TPR: {}'.format(fprOpt, tprOpt))\n\n# Create data viz\nplotnine.options.figure_size = (8, 4.8)\n(\n    ggplot(data = df_fpr_tpr)+\n    geom_point(aes(x = 'FPR',\n                   y = 'TPR'),\n               size = 0.4)+\n    # Best threshold\n    geom_point(aes(x = fprOpt,\n                   y = tprOpt),\n               color = '#981220',\n               size = 4)+\n    geom_line(aes(x = 'FPR',\n                  y = 'TPR'))+\n    geom_text(aes(x = fprOpt,\n                  y = tprOpt),\n              label = 'Optimal threshold \\n for class: {}'.format(thresholdOpt),\n              nudge_x = 0.14,\n              nudge_y = -0.10,\n              size = 10,\n              fontstyle = 'italic')+\n    labs(title = 'ROC Curve')+\n    xlab('False Positive Rate (FPR)')+\n    ylab('True Positive Rate (TPR)')+\n    theme_minimal()\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:32.154817Z","iopub.execute_input":"2022-05-17T13:48:32.155083Z","iopub.status.idle":"2022-05-17T13:48:37.239144Z","shell.execute_reply.started":"2022-05-17T13:48:32.155039Z","shell.execute_reply":"2022-05-17T13:48:37.238190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:37.240928Z","iopub.execute_input":"2022-05-17T13:48:37.241591Z","iopub.status.idle":"2022-05-17T13:48:37.246441Z","shell.execute_reply.started":"2022-05-17T13:48:37.241507Z","shell.execute_reply":"2022-05-17T13:48:37.245716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\n# Array for finding the optimal threshold\nthresholds = np.arange(0.0, 1.0, 0.0001)\nfscore = np.zeros(shape=(len(thresholds)))\nprint('Length of sequence: {}'.format(len(thresholds)))\n\n# Fit the model\nfor index, elem in enumerate(thresholds):\n    # Corrected probabilities\n    y_pred_prob = (train_proba[:, 1] > elem).astype('int')\n    # Calculate the f-score\n    fscore[index] = f1_score(y[:, 0], y_pred_prob)\n\n# Find the optimal threshold\nindex = np.argmax(fscore)\nthresholdOpt = round(thresholds[index], ndigits = 4)\nfscoreOpt = round(fscore[index], ndigits = 4)\nprint('Best Threshold: {} with F-Score: {}'.format(thresholdOpt, fscoreOpt))\n\n# Plot the threshold tuning\ndf_threshold_tuning = pd.DataFrame({'F-score':fscore,\n                                    'Threshold':thresholds})\ndf_threshold_tuning.head()\n\nplotnine.options.figure_size = (8, 4.8)\n(\n    ggplot(data = df_threshold_tuning)+\n    geom_point(aes(x = 'Threshold',\n                   y = 'F-score'),\n               size = 0.4)+\n    # Best threshold\n    geom_point(aes(x = thresholdOpt,\n                   y = fscoreOpt),\n               color = '#981220',\n               size = 4)+\n    geom_line(aes(x = 'Threshold',\n                   y = 'F-score'))+\n    # Annotate the text\n    geom_text(aes(x = thresholdOpt,\n                  y = fscoreOpt),\n              label = 'Optimal threshold \\n for class: {}'.format(thresholdOpt),\n              nudge_x = 0,\n              nudge_y = -0.10,\n              size = 10,\n              fontstyle = 'italic')+\n    labs(title = 'Threshold Tuning Curve')+\n    xlab('Threshold')+\n    ylab('F-score')+\n    theme_minimal()\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:48:37.248027Z","iopub.execute_input":"2022-05-17T13:48:37.248297Z","iopub.status.idle":"2022-05-17T13:49:42.228251Z","shell.execute_reply.started":"2022-05-17T13:48:37.248250Z","shell.execute_reply":"2022-05-17T13:49:42.227218Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"threshold = [np.mean(th1), np.mean(th2)]\nnew_pred = []\nfor i in range(train_pred.shape[0]):\n    new_pred.append(re_predict(train_proba[i, :], threshold))\nprint('\\nPerformance Metrics after Over Prediction')\nprint('1. The F-1 score of the model {}\\n'.format(f1_score(ytrain, new_pred, average='macro')))\nprint('2. The recall score of the model {}\\n'.format(recall_score(ytrain, new_pred, average='macro')))\nprint('3. The Matthews Correlation Coefficient: {}\\n'.format(matthews_corrcoef(ytrain, new_pred)))\nprint('4. Classification report \\n {} \\n'.format(classification_report(ytrain, new_pred)))\nprint('5. Confusion matrix \\n {} \\n'.format(confusion_matrix(ytrain, new_pred)))","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:49:42.230222Z","iopub.execute_input":"2022-05-17T13:49:42.230920Z","iopub.status.idle":"2022-05-17T13:49:42.359999Z","shell.execute_reply.started":"2022-05-17T13:49:42.230845Z","shell.execute_reply":"2022-05-17T13:49:42.358791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Repredict with new thresholds","metadata":{"_uuid":"87d765e4ef847d66fb8b0a06b33f8b2db27e723f"}},{"cell_type":"code","source":"threshold = [np.median(th1), np.median(th2)]\nnew_pred = []\nfor i in range(train_pred.shape[0]):\n    new_pred.append(re_predict(train_proba[i, :], threshold))\nprint('\\nPerformance Metrics after Over Prediction')\nprint('1. The F-1 score of the model {}\\n'.format(f1_score(ytrain, new_pred, average='macro')))\nprint('2. The recall score of the model {}\\n'.format(recall_score(ytrain, new_pred, average='macro')))\nprint('3. The Matthews Correlation Coefficient: {}\\n'.format(matthews_corrcoef(ytrain, new_pred)))\nprint('4. Classification report \\n {} \\n'.format(classification_report(ytrain, new_pred)))\nprint('5. Confusion matrix \\n {} \\n'.format(confusion_matrix(ytrain, new_pred)))","metadata":{"execution":{"iopub.status.busy":"2022-05-17T13:49:42.362201Z","iopub.execute_input":"2022-05-17T13:49:42.362621Z","iopub.status.idle":"2022-05-17T13:49:42.493780Z","shell.execute_reply.started":"2022-05-17T13:49:42.362526Z","shell.execute_reply":"2022-05-17T13:49:42.492721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create submission file","metadata":{"_uuid":"78a3ae3d3beadbc89b67ddee01c5f2620e4a95d5"}},{"cell_type":"markdown","source":"From, the forums it appears that most people with models scoring greater than 0.6 are using deep learning methods and neural networks instead of relying on feature-based approaches that require signal processing.\n\n## Improvements\nThere are notably at least two promising paths to improve performance of partial discharge fault detection models, this model included.\n\n### Transition Data Collection and Processing from Single Cycle to a N Cycle Accumulator\nMoving from collecting and processing one 50Hz grid cycle to sampling N consecutive samples would benefit partial discharge faults in a number of ways. First and foremost, the sample size is larger. But more importantly it offers the opportunity to improve detection of persistent partial discharge patterns while mitigating sporadic effects like corona discharge. At a minimum, measurement data should be tagged with a collector ID so that the same length of covered conductor can be assessed over time. This data should come for free since presumably, if a fault is detected, the power company is going to want to know where to dispatch repair crews.\n\n### Leverage Grouped Faults in the Prediction Model\nWhile not universally true, when a measurement contains a fault in one of the signal phases, a fault is present on at least one of the other phases as well. In the provided, labeled training data 90.1% of the time, this observation holds true. It may be possible to implement a M of N voting system for better fault detection but care must be taken to not introduce additional error and bias.  \n","metadata":{"_uuid":"7f8de859f25390654c41ab53c9aa8e9303f74df3"}}]}