{"cells":[{"metadata":{"trusted":true,"_uuid":"6dfda152237bf6a34931d3f68cfce3e166552f26"},"cell_type":"markdown","source":"**The purpose of this notebook is to demonstrate the use of autoencoding to extract relevant data from a signal. **\n\nTraining usually consists of taking data that has been manually classified, and using it to train an artificial neural network (ann). When each example contains large amounts of data, we sometimes extract features before training (transformation), or use multi-layer ann to automatically extract features. \n\nThe former has issues because we may discard important information that takes place at a tiny detailed scale of the data, in favor of dimensionality reduction.\n\nThe latter can be cumbersome and time consuming for two reasons; first due to the large size of input layers - because we may need information from the small scale, and second due to the large depth of the ann - because large scale info may be important to correctly classify the data.\n\nInstead, this notebook proposes to train a artificial neural network to simply compress the data, without regard to the manual classification - thereby learning the patterns that are common to all of the training signals. We then subtract this \"common data\" from the individual training signal, and only use the remaining portion (the residual) to perform classification.\n\n### **Theory:**\n#### You are probably familiar with existing methods to extract basic time-varying signal information through the use of transformation. **Fourier** transformation is often used for periodic signals, where the convolution is with various period sine and cosine functions. **Wavelet** transformation is often used for time limited signals, where convolution is with a wavelet shape that has been selected by hand. Note that with wavelet decomposition, the scale of the wavelet shape is changed, similar to the way the period of the Fourier functions are changed. **Convolutional** networks seek to learn a set of optimal wavelets. Unlike traditional wavelet transformation, convolutional networks do not change the scale of the wavelet shape, and instead rely on layers of the network to learn patterns of various scale. In all of these cases, a compressed \"lossy\" representation of the original signal is created. Dense networks can also be used to learn patterns within time varying signals, and also create a compressed version of the original signal.\n\n#### A **Dense** neural network layer is used to compress the original signal of nearly 1MB down to a few floating point numbers (less than 1kB). Using this compressed information, a lossy version of the original signal can be reconstructed. The difference between the original signal, and the reconstructed signal is called the **Residual**. The process can be repeated in an attempt to compress the residual, similar to the method of successive wavelet decomposition.\n\n#### Along the way, we can also see interesting patterns emerge. By using visualization, we can compare the original signal to successive de-compressed versions, and see what is common across all signals, and also what is different between them.\n\n### **Revision Info:**\nPAN Dec. 21, 2018\n#### Forked from Panchajanya Banerjee (Pancham) - First Steps EDA https://www.kaggle.com/delayedkarma/first-steps-eda\n#### V09 fixed typo in feature extraction code\n#### V08 - tuned dense nn\n#### V07 switched from random forest as final classifier to Dense NN\n#### V05 and V06 fixing memory overflow errors in Kaggle, added comments to make it more \"educational\" (hopefully)\n#### V04 - Improved residual compression MSE by using multi-layer, added better feature extraction of 2nd residual using example code from VSB Power LSTM attention https://www.kaggle.com/braquino/vsb-power-lstm-attention by Bruno Aquino\n#### V02 and V03 - Added 2nd residual features\n#### V01 - Intended for release to the public as an educational tool Jan 2019\n\n### **Training and Testing Data Info:**\nThe data being used - comes from the VSB Power contest on kaggle.com (January, 2019)\n#### Each id_measurement [train 0-2903, test 2904-9682]  consists of three signal_id's [train 0-8711, test 8712-29048] where each signal_id is associated with one of three phase's [0-2] also, the training set also provides a target [0-1] for each signal_id indicating a fault or not \n#### There are 800,000 int8 samples for each signal_id [train 0-8711, test 8712-29048]\n\n### **Possible next steps:**\n#### Sloppy \"cut and paste\" coding should be replaced by function calls\n#### The final residual data can be further analyzed for feature extraction\n#### If enough folks are interested, I can make the documentation more rigorous with citations and formulas"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# LOADING UP PYTHON COMPONENTS\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\nimport pyarrow.parquet as pq\nimport os\nprint(os.listdir(\"../input\"))\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport tensorflow as tf\nfrom tensorflow.python.keras.models import Sequential, Model, load_model\nfrom tensorflow.python.keras.layers import Dense","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# VIEWING THE FIRST SIX TRAINING SIGNALS\n# The first row of three shows the three phases of a power line that has been hand-classified as not having a fault\n# The second row of three shows the three phases of a power line that has a fault\nmeta_train = pd.read_csv('../input/metadata_train.csv')\n# %%time \n# Read in the first two signals (three phases each) for display.  \n# Each column contains one signal\nsubset_train = pq.read_pandas('../input/train.parquet', columns=[str(i) for i in range(6)]).to_pandas()\n# Comparing a good signal (forst row, all three phases)\n# with a bad signal (second row, all three phases)\nfig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2,3, figsize=(12,12))\nax1.plot(subset_train['0']) ;\nax2.plot(subset_train['1']) ;\nax3.plot(subset_train['2']) ;\nax4.plot(subset_train['3']) ;\nax5.plot(subset_train['4']) ;\nax6.plot(subset_train['5']) ;\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ceb6cd605fe4db0b91c92a6567aec39166a2497b"},"cell_type":"code","source":"# HAND-CLASSIFIED INFORMATION ABOUT THE ABOVE SIGNALS\nprint(\"Classifications of above\")\nmeta_train[0:6]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d085923ea0804d45c9750f0154122a00c6736bad"},"cell_type":"code","source":"%%time\n### Load Raw Training Data\nbegin_col = 0\n# number of test examples to use\n# num_to_use = 8712\n# use a smaller subset because of Kaggle RAM limitations \nnum_to_use = 2001\n# num_to_use = 30\n\nfilename = '../input/train.parquet'\n\nX = pq.read_pandas(filename, columns=[str(j + begin_col) for j in range(num_to_use)]).to_pandas().values.transpose()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ca0a541bb8e8102928843cc3e422d462ba241882"},"cell_type":"code","source":"# the data we will use to train the autoencoder\nX.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f1ff70f5f2fef9958edc81e76f90148c06574873"},"cell_type":"code","source":"# The autoencoder network\ncompression_1 = 5 \n\nmodel = Sequential()\nmodel.add(Dense(compression_1, activation='relu', input_shape=(800000,), name='compress'))\n# this output layer has to have 800,000 neurons, and needs to be linearly activated\nmodel.add(Dense(800000, activation='linear'))\nmodel.compile(loss='mean_squared_error', optimizer='adam')\n# Diplay the model summary\nprint(\"model summary\")\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e8b661ab9852c4829a5521339ff77b023dc5ba49","scrolled":false},"cell_type":"code","source":"# Train the autoencoder to minimize mean square error\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\n\nearlystopper = EarlyStopping(patience=2, verbose=1) \ncheckpointer = ModelCheckpoint('VSBautoassoc', verbose=1, save_best_only=True)\nresults = model.fit(X, X, validation_split=0.1, batch_size=50, epochs=300, \n                    callbacks=[earlystopper, checkpointer])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a020f1b89aae91b2d94189bc55963fb8c520bec1"},"cell_type":"code","source":"model = load_model('VSBautoassoc')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"81f59d68809a929d07cf4b1416a7a4292321f9a1"},"cell_type":"code","source":"# The front half of the autoencoder is the \"coder\" part of this CODEC pair\n# We can use the coder to convert the very large signal data (800000 integers) into a much smaller compressed version\n# Load just the compression model\nc_model = Sequential()\nc_model.add(Dense(compression_1, activation='relu', input_shape=(800000,), name='compress'))\n\nc_model.load_weights('VSBautoassoc', by_name=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e1a9455c2a3f9ada9b37cc403080c016ebf1e772"},"cell_type":"code","source":"# Since the autoencoder both codes and decodes, we can use it to create the lossy approximation of the original signal\n# Looking at those first six signals again, overlaying their compressed version on top of the original signal\nX_decompressed = model.predict(X[0:6])\nX_decompressed.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7bed3b22851e6e767c015e29ce2ca78b45d05ce1"},"cell_type":"code","source":"# Here we are comparing the original signal to its lossy reconstruction\n# Notice how only the elements most common to ALL the signals in the training set are represented\nfig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2,3, figsize=(12,12))\nax1.plot(subset_train['0']) ;\nax2.plot(subset_train['1']) ;\nax3.plot(subset_train['2']) ;\nax4.plot(subset_train['3']) ;\nax5.plot(subset_train['4']) ;\nax6.plot(subset_train['5']) ;\nax1.plot(X_decompressed[0]) ;\nax2.plot(X_decompressed[1]) ;\nax3.plot(X_decompressed[2]) ;\nax4.plot(X_decompressed[3]) ;\nax5.plot(X_decompressed[4]) ;\nax6.plot(X_decompressed[5]) ;","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ab06801693446b1e48e558065bdfbe8aed597988"},"cell_type":"code","source":"# the residual is what remains after subtracting the recreated signal from the original.\n# In theorythis should contain the information that makes the training signals DIFFERENT from each other\n# we will be repeating this iteratively, but for now, let's look at the residual\nX_residual = X[0:6] - X_decompressed[0:6]\nX_residual.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eeed0af65b065b597f6a54e1a9f237522403e21c"},"cell_type":"code","source":"# Here we see the residual superimposed on the original signal\n# As expected, the residual removes the unimportant information, and leaves the possibly important noise spikes\n# Nevertheless, a strange seventh harmonic (ripple with seven peaks) shows up on these first two sets of 3-phase samples (?)\nfig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2,3, figsize=(12,12))\nax1.plot(subset_train['0']) ;\nax2.plot(subset_train['1']) ;\nax3.plot(subset_train['2']) ;\nax4.plot(subset_train['3']) ;\nax5.plot(subset_train['4']) ;\nax6.plot(subset_train['5']) ;\nax1.plot(X_residual[0]) ;\nax2.plot(X_residual[1]) ;\nax3.plot(X_residual[2]) ;\nax4.plot(X_residual[3]) ;\nax5.plot(X_residual[4]) ;\nax6.plot(X_residual[5]) ;","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9bbb473350dbe096fc43f4e73993dadb06c903cc"},"cell_type":"code","source":"# we will now create a new autoencoder, to compress those residual signals\n# We are repeating this a second time in case there is still useless information that is very common between all signals\n# There is little risk of \"throwing away the baby with the bath water\" because the \"codes\" or compressed version of the info\n# will still be kept and used for classification training.\n#\n# so we begin by gathering all the residuals\n### CONVERT X into FIRST Residuals (to save memory)\nX.shape\nbegin_col = 0\nbatch = 600 \nnum_batches = int((num_to_use) / batch)\nremainder = int((num_to_use) % batch)\nX_decompressed = np.zeros([batch, 800000])\nif (num_batches > 0) :\n    for ix in range (num_batches) :\n        X_decompressed[int(0):int(batch)] = model.predict(X[int(ix*batch):int((ix+1)*batch)])\n        X[int(ix*batch):int((ix+1)*batch)] = X[int(ix*batch):int((ix+1)*batch)] - X_decompressed[int(0):int(batch)]\nif (remainder > 0) :\n    ix = num_batches\n    X_decompressed[int(0):int(remainder)] = model.predict(X[int(ix*batch):int((ix*batch)+remainder)])\n    X[int(ix*batch):int((ix*batch)+remainder)] = X[int(ix*batch):int((ix*batch)+remainder)] - X_decompressed[int(0):int(remainder)]\n\nX_decompressed = 5 # quick and dirty data=cleanup\nX.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"09b3e6a5819c531ceec100222fe6626c07d0c388"},"cell_type":"code","source":"# The residuals are a bit more complex than the original signals, and so we'll use a multilayer network for that compression\n### As with the original signal, we compress the residuals using an autoencoder\ncompression_2 = 20\n\nr_model = Sequential()\nr_model.add(Dense(compression_2 * 2, activation='relu', input_shape=(800000,), name='layer1'))\nr_model.add(Dense(compression_2 * 4, activation='relu', name='layer2'))\nr_model.add(Dense(compression_2, activation='relu', name='compressed'))\nr_model.add(Dense(compression_2 * 4, activation='relu', name='layer4'))\nr_model.add(Dense(compression_2 * 2, activation='relu', name='layer5'))\nr_model.add(Dense(800000, activation='linear', name='output'))\n\nr_model.compile(loss='mean_squared_error', optimizer='adam')\n# Diplay the model summary\nprint(\"model summary\")\nr_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b72133dd3d87415e34d2a008adf9aac3c2a91aa5"},"cell_type":"code","source":"### Train the residual compressor\nearlystopper = EarlyStopping(patience=2, verbose=1) \ncheckpointer = ModelCheckpoint('VSB_r_autoassoc', verbose=1, save_best_only=True)\nresults = r_model.fit(X, X, validation_split=0.1, batch_size=50, epochs=300, \n                    callbacks=[earlystopper, checkpointer])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d6f7f3029262635615c7ed6922b63947946cba2b"},"cell_type":"code","source":"r_model = load_model('VSB_r_autoassoc')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"790a5368d99a4e38a29ab4dac7ea88d66549c8f8"},"cell_type":"code","source":"# Load just the compression portion of the autoencoder\nc_r_model = Sequential()\nc_r_model.add(Dense(compression_2 * 2, activation='relu', input_shape=(800000,), name='layer1'))\nc_r_model.add(Dense(compression_2 * 4, activation='relu', name='layer2'))\nc_r_model.add(Dense(compression_2, activation='relu', name='compressed'))\n\nc_r_model.load_weights('VSB_r_autoassoc', by_name=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"943224d6e296f3d9c93528bb66a5fa6f24b69ab7"},"cell_type":"code","source":"# Use the complete autoencoder to create lossy compressed versions of the residuals\nX_decompressed = r_model.predict(X[0:6])\nX_decompressed.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1a26d6a7cc40e56eafda7ae1927ba6197dfaf13c"},"cell_type":"code","source":"# Display these decompressed residuals, showing difference from the actual residuals\n# Here we should see any remaining common elements represented by the residual approximation\nfig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2,3, figsize=(12,12))\nax1.plot(X[0]) ;\nax2.plot(X[1]) ;\nax3.plot(X[2]) ;\nax4.plot(X[3]) ;\nax5.plot(X[4]) ;\nax6.plot(X[5]) ;\nax1.plot(X_decompressed[0]) ;\nax2.plot(X_decompressed[1]) ;\nax3.plot(X_decompressed[2]) ;\nax4.plot(X_decompressed[3]) ;\nax5.plot(X_decompressed[4]) ;\nax6.plot(X_decompressed[5]) ;","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f926a0b9ea0b2a8108fd8d47f0ab41b57a50682c"},"cell_type":"code","source":"# What about the residual of the residual? This is also called the 2nd decomposition, in wavelet transformation lingo\n# This is (hopefully) the distilled information that makes every signal different from all the others\n# becuase it was not captured by the \"approximations\" learned by the previous two autoencoders\nX_residual = X[0:6] - X_decompressed[0:6]\nX_residual.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e4ed15ed2ab07bb0f60ba76eef75691926488658"},"cell_type":"code","source":"# This 2nd residual is now displayed, superimposed on the original signal\n# Visually / qualitatively, it should look like the \"important\" information, if our process is working correctly\nfig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2,3, figsize=(12,12))\nax1.plot(subset_train['0']) ;\nax2.plot(subset_train['1']) ;\nax3.plot(subset_train['2']) ;\nax4.plot(subset_train['3']) ;\nax5.plot(subset_train['4']) ;\nax6.plot(subset_train['5']) ;\nax1.plot(X_residual[0]) ;\nax2.plot(X_residual[1]) ;\nax3.plot(X_residual[2]) ;\nax4.plot(X_residual[3]) ;\nax5.plot(X_residual[4]) ;\nax6.plot(X_residual[5]) ;","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2b0ab5fbf2dcaf69f07b140dbebb215cfdee4ea0"},"cell_type":"code","source":"# Some global parameters\n# number of non-overlapping time slices within the signal that will be used to extract feature information\nwindows = 80\n# number of features per time slice\nwin_feat = 7\n# total number of features extracted from the signal\nnumber_of_features = int(compression_1 + compression_2 + (windows * win_feat))\n\n# Now batch convert original signals into training feature vectors\n# the feature vectors consist of three parts:\n# * compressed signal (a handful of floating point numbers)\n# * compressed first residual (some more floating point numbers)\n# * N*F features of the N windows of 800,000 / N samples of the 2nd residual \n#     put another way, this last part is F floating point numbers containing statistical info about the 2nd residual\n#     after we chop up that 2nd residual into N different non-overlapping time windows\ndef extract_features(begin_col, loc_num_to_use, filename) :  \n    batch = 300 \n    num_batches = int((loc_num_to_use) / batch)\n    remainder = int((loc_num_to_use) % batch)\n    win_size = int(800000 / windows)\n    ### Create a pandas data frame\n    loc_X = np.zeros((int(loc_num_to_use), number_of_features))\n    if num_batches > 0:\n        for ix in range (num_batches) :\n            # load a batch of signals\n            x1 = pq.read_pandas(filename, columns=[str(ix * batch + j + begin_col) for j in range(batch)]).to_pandas().values.transpose()\n            # compress them into 20 data values\n            c_x1 = c_model.predict(x1)\n            # recreate them from the compressed data\n            x1_r = model.predict(x1)\n            # residual is the original signal minus the recreated one\n            res_x1 = x1 - x1_r\n            # compress the residual\n            c_res_x1 = c_r_model.predict(res_x1)\n            # recreate the residuals from the compressed data\n            res_x1_r = r_model.predict(res_x1)\n            # second residual is the residual minus the recreated residual\n            res_2_x1 = res_x1 - res_x1_r\n            for j in range (0,batch) :\n                i = ix * batch + j\n                loc_X[i,0:compression_1] = c_x1[j]\n                loc_X[i,compression_1:compression_1 + compression_2] = c_res_x1[j]\n                for win in range (windows) :\n                    # start and end of window in signal data\n                    win_start = win * win_size\n                    win_end = win_start + win_size\n                    # start of windows features in feature array (loc_X)\n                    win_fs = win * win_feat\n                    ### V04 ----------\n                    ### Using example code from VSB Power LSTM attention https://www.kaggle.com/braquino/vsb-power-lstm-attention by Bruno Aquino\n                    # Mean\n                    loc_X[i,compression_1 + compression_2 + win_fs + 0] = mean = res_2_x1[j,win_start:win_end].mean()\n                    # Standard Deviation = sqrt(variance)\n                    loc_X[i,compression_1 + compression_2 + win_fs + 1] = std  = res_2_x1[j,win_start:win_end].std()\n                    # top of standard deviation range\n                    loc_X[i,compression_1 + compression_2 + win_fs + 2] = mean + std\n                    # bottom of standard deviation range\n                    loc_X[i,compression_1 + compression_2 + win_fs + 3] = mean - std\n                    # calculate a handful of percentiles\n                    pct_calc = np.percentile(res_2_x1[j,win_start:win_end], [0, 1, 25, 50, 75, 99, 100])\n                    # max range of percentiles\n                    loc_X[i,compression_1 + compression_2 + win_fs + 4] = pct_calc[-1] - pct_calc[0]\n                    # coefficient of variation (standard deviation divided by mean)\n                    loc_X[i,compression_1 + compression_2 + win_fs + 5] = std / mean\n                    # A measure of asymmetry (75th percentile subtracted from mean)\n                    loc_X[i,compression_1 + compression_2 + win_fs + 6] = mean - pct_calc[4]\n                    # the seven percentile values calcuated earlier\n                    ### end of example code for suggested features to extract\n    ix = num_batches\n    # load a batch of signals\n    x1 = pq.read_pandas(filename, columns=[str(ix * batch + j + begin_col) for j in range(batch)]).to_pandas().values.transpose()\n    # compress them into 20 data values\n    c_x1 = c_model.predict(x1)\n    # recreate them from the compressed data\n    x1_r = model.predict(x1)\n    # residual is the original signal minus the recreated one\n    res_x1 = x1 - x1_r\n    # compress the residual\n    c_res_x1 = c_r_model.predict(res_x1)\n    # recreate the residuals from the compressed data\n    res_x1_r = r_model.predict(res_x1)\n    # second residual is the resodual minus the recreated residual\n    res_2_x1 = res_x1 - res_x1_r\n                   \n    for j in range (0,remainder) :\n        i = ix * batch + j\n        loc_X[i,0:compression_1] = c_x1[j]\n        loc_X[i,compression_1:compression_1 + compression_2] = c_res_x1[j]\n        for win in range (windows) :\n            win_start = win * win_size\n            win_end = win_start + win_size\n            # start of windows features in feature array (loc_X)\n            win_fs = win * win_feat\n            ### V04 ----------\n            ### Using example code from VSB Power LSTM attention https://www.kaggle.com/braquino/vsb-power-lstm-attention by Bruno Aquino\n            # Mean\n            loc_X[i,compression_1 + compression_2 + win_fs + 0] = mean = res_2_x1[j,win_start:win_end].mean()\n            # Standard Deviation = sqrt(variance)\n            loc_X[i,compression_1 + compression_2 + win_fs + 1] = std  = res_2_x1[j,win_start:win_end].std()\n            # top of standard deviation range\n            loc_X[i,compression_1 + compression_2 + win_fs + 2] = mean + std\n            # bottom of standard deviation range\n            loc_X[i,compression_1 + compression_2 + win_fs + 3] = mean - std\n            # calculate a handful of percentiles\n            pct_calc = np.percentile(res_2_x1[j,win_start:win_end], [0, 1, 25, 50, 75, 99, 100])\n            # max range of percentiles\n            loc_X[i,compression_1 + compression_2 + win_fs + 4] = pct_calc[-1] - pct_calc[0]\n            # coefficient of variation (standard deviation divided by mean)\n            loc_X[i,compression_1 + compression_2 + win_fs + 5] = std / mean\n            # A measure of asymmetry (75th percentile subtracted from mean)\n            loc_X[i,compression_1 + compression_2 + win_fs + 6] = mean - pct_calc[4]\n            # the seven percentile values calcuated earlier\n            ### end of example code for suggested features to extract\n                \n    return loc_X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"99f4350a3ee9d3d9015eea6a94373b13ff1c354d"},"cell_type":"code","source":"%%time\n# Now we are recreating what we did above for the first six samples, but for the entire training set\n# encoding the signals, encoding the residuals, and gathering statistical info about each 2nd residual\nnum_to_use = 8712\n# num_to_use = 30\n\nfrom sklearn import preprocessing\n\n# extract features\nX_unscaled = extract_features(0, num_to_use, '../input/train.parquet')\n\nscaler = preprocessing.StandardScaler().fit(X_unscaled)\nX = scaler.transform(X_unscaled)\n\n# correct classifications from training set\ny = np.zeros((num_to_use))\nfor i in range(0, int(num_to_use)):\n    y[i] = meta_train.target[i]\n    \nprint(y.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f2e8f6c58c01ec316502582357c0a19e5835e68c"},"cell_type":"code","source":"# This is not really necessary, but it gives me a good hint if the problem is being solved correctly\n# How often do power line faults occur on all three phases simultaneously, versus on fewer than all 3?\ntriples = 0\ndoubles = 0\nsingles = 0\nfor i in range(0,int(num_to_use),3) : \n    if (meta_train.target[i] and meta_train.target[i+1] and meta_train.target[i+2] ):\n#        print('triple',meta_train.signal_id[i], meta_train.phase[i] )\n        triples = triples + 1\n    elif (meta_train.target[i] + meta_train.target[i+1] + meta_train.target[i+2] == 2):\n#        print('double',meta_train.signal_id[i], meta_train.phase[i])\n        doubles = doubles + 1\n    elif (meta_train.target[i] + meta_train.target[i+1] + meta_train.target[i+2] == 1):\n#        print('single',meta_train.signal_id[i], meta_train.phase[i])\n        singles = singles + 1\n\nprint('triples', triples, 'doubles', doubles, 'singles', singles)\nprint('sanity check: ', 'total faults', meta_train.target[0:int(num_to_use)].sum(), ' sum of above ', 3 * triples + 2 * doubles + singles)\n# plt.plot(meta_train.target)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c971d263193c500574a7d6044271eca762628a63"},"cell_type":"code","source":"\"\"\"%%time\n\n\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import train_test_split\n\nfirst_model = RandomForestRegressor(n_estimators=30, min_samples_leaf=30, \n                                    random_state=1).fit(X, y)\n# Environment Set-Up for feedback system.\nfrom learntools.core import binder\nbinder.bind(globals())\nfrom learntools.ml_insights.ex2 import *\nprint(\"Training Complete\")\npredicted_y = first_model.predict(X)\n\"\"\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"faf820c5e6bfb2aeddb454ef8491da33719fe4e0"},"cell_type":"code","source":"# Now we train a multi-layer ann to learn the correct classifications\nclass_model = Sequential()\nclass_model.add(Dense(number_of_features, activation='relu', input_shape=(number_of_features,)))\nclass_model.add(Dense(number_of_features * 2, activation='relu'))\nclass_model.add(Dense(number_of_features * 4, activation='relu'))\nclass_model.add(Dense(number_of_features * 2, activation='relu'))\nclass_model.add(Dense(number_of_features, activation='relu'))\nclass_model.add(Dense(1, activation='sigmoid', name='output'))\n\nclass_model.compile(loss='binary_crossentropy', optimizer='adam')\n# Diplay the model summary\nprint(\"model summary\")\nclass_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82a6a52776433bff1f23a1d94ada9f30a7cdd998"},"cell_type":"code","source":"### Train the residual compressor\nearlystopper = EarlyStopping(patience=5, verbose=1) \ncheckpointer = ModelCheckpoint('VSB_classifier', verbose=1, save_best_only=True)\nresults = class_model.fit(X, y, validation_split=0.2, batch_size=50, epochs=300, \n                    callbacks=[earlystopper, checkpointer])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0f92c3ec6072d77d86dbd3798c0e6f9fab3b0979"},"cell_type":"code","source":"class_model = load_model('VSB_classifier')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"40fe41a0652eeca370dda60dd836a8b95e0bd087"},"cell_type":"code","source":"predicted_y = class_model.predict(X)\npredicted_y.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ff1f617577f4e12d6ca7f759979c3e75ae955f88"},"cell_type":"code","source":"### See if the best_threshold is indeed the best by using the Matthews Correlation Coefficient\n### The Matthews correlation coefficient\n### MCC=( TP * TN - FP * FN ) / sqrt  ( (TP+FP) * (TP+FN) * (TN+FP) * (TN+FN) )\n\nimport math \n\nbest_thresh = 0\nnum_threshes = 1000           # number of different thresholds to try\nmin_thresh = (1 / num_threshes) / 2 # the minimum threshold to try\nbest_MCC = 0\nTPs = np.zeros(num_threshes)\nTNs = np.zeros(num_threshes)\nFPs = np.zeros(num_threshes)\nFNs = np.zeros(num_threshes)\nMCCs = np.zeros(num_threshes)\nthreshes = np.zeros(num_threshes)\n# print (\"threshold, Data count, TP, TN, FP, FN, TP + TN + FP + FN, MCC\")\nfor this_thresh in range (num_threshes):\n    thresh = round(min_thresh + (this_thresh / num_threshes), 3)\n    TP = 0\n    TN = 0\n    FP = 0\n    FN = 0\n    for i in range(0, int(num_to_use)):\n        if (y[i]) :\n            if (predicted_y[i] > thresh) :\n                TP = TP + 1\n            else :\n                FN = FN + 1\n        else:\n            if (predicted_y[i] > thresh) :\n                FP = FP + 1\n            else :\n                TN = TN + 1\n        if (math.sqrt ( (TP + FP) * (TP + FN) * (TN + FP) * (TN + FP))):\n            MCC = ((TP * TN) - (FP - FN)) / math.sqrt ( (TP + FP) * (TP + FN) * (TN + FP) * (TN + FP) )\n        else :\n            MCC = 0\n    TPs[this_thresh] = TP\n    TNs[this_thresh] = TN\n    FPs[this_thresh] = FP\n    FNs[this_thresh] = FN\n    MCCs[this_thresh] = MCC\n    threshes[this_thresh] = thresh\n#    print(thresh, num_to_use * 3, TP, TN, FP, FN, (TP+TN+FP+FN), MCC)\n    if (MCC >= best_MCC):\n        best_MCC = MCC\n        best_thresh = thresh\n\nfig, ((ax1, ax2), (ax3, ax4), (ax5,ax6)) = plt.subplots(3,2, figsize=(20, 20))\nax1.plot(threshes, TPs) \nax1.set_title(\"True Positives\");\nax2.plot(threshes, TNs) ;\nax2.set_title(\"True Negatives\");\nax3.plot(threshes, FPs) ;\nax3.set_title(\"False Positives\");\nax4.plot(threshes, FNs) ;\nax4.set_title(\"False Negatives\");\nax5.plot(threshes, MCCs) ;\nax5.set_title(\"The Matthews correlation coefficient\");\nax6.plot(threshes) ;\nax6.set_title(\"Threshod Values\");\n        \nprint(' ')\nprint('best threshold, ', best_thresh, ' yielding best MCC, ', best_MCC)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"65055b4d4e0aabb311ca529d446a4b61e9addb60"},"cell_type":"code","source":"y = predicted_y\n# See how many doubles and triples we got on the test data\n# How often did we predict power line faults occur on all three phases simultaneously, versus on fewer than all 3?\ntriples = 0\ndoubles = 0\nsingles = 0\nfor i in range(0,int(num_to_use),3) : \n    y[i] = int(y[i] > best_thresh) \n    y[i+1] = int(y[i+1] > best_thresh)\n    y[i+2] = int(y[i+2] > best_thresh)\n    num_phases_faulty = y[i] + y[i+1] + y[i+2]\n    if (num_phases_faulty == 3):\n#        print('triple',meta_train.signal_id[i], meta_train.phase[i] )\n        triples = triples + 1\n    elif (num_phases_faulty == 2):\n#        print('double',meta_train.signal_id[i], meta_train.phase[i])\n        doubles = doubles + 1\n    elif (num_phases_faulty == 1):\n#        print('single',meta_train.signal_id[i], meta_train.phase[i])\n        singles = singles + 1\n\nprint('triples', triples, 'doubles', doubles, 'singles', singles)\nprint('sanity check: ', 'total faults', y.sum(), ' sum of above ', 3 * triples + 2 * doubles + singles)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2fcb8001ed787354a7efa185da6be2da3ab94b6e"},"cell_type":"code","source":"# quick and dirty RAM cleanup\nX_unscaled = 5\nX = 5\ny = 5\npredicted_y = 5\nTPs = 5\nTNs = 5\nFPs = 5\nFNs = 5\nMCCs = 5\nthreshes = 5\nresults = 5\nX_residual = 5\nX_decompressed = 5","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f827365b2b96d174a5dd033406aeac276afde986"},"cell_type":"code","source":"%%time \n# Now load the test data\nmeta_test = pd.read_csv('../input/metadata_test.csv')\n\n# Calculate Test Value Parameters\n\n# number of test examples to use \nnum_to_use = 20337\n# use a smaller subset for testing if you don't want to wait (must be a multiple of 3)\n# num_to_use = 30\n\nX = 5\nX_unscaled = 5\n\nX_unscaled = extract_features(8712, num_to_use, '../input/test.parquet')\n\nX = scaler.transform(X_unscaled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8412016cd69798748e74be9f364fa7705b920cc9"},"cell_type":"code","source":"y = class_model.predict(X)\n\n# How often did we predict power line faults occur on all three phases simultaneously, versus on fewer than all 3?\ntriples = 0\ndoubles = 0\nsingles = 0\nfor i in range(0,int(num_to_use),3) : \n    y[i] = int(y[i] > best_thresh) \n    y[i+1] = int(y[i+1] > best_thresh)\n    y[i+2] = int(y[i+2] > best_thresh)\n    num_phases_faulty = y[i] + y[i+1] + y[i+2]\n    if (num_phases_faulty == 3):\n#        print('triple',meta_train.signal_id[i], meta_train.phase[i] )\n        triples = triples + 1\n    elif (num_phases_faulty == 2):\n#        print('double',meta_train.signal_id[i], meta_train.phase[i])\n        doubles = doubles + 1\n    elif (num_phases_faulty == 1):\n#        print('single',meta_train.signal_id[i], meta_train.phase[i])\n        singles = singles + 1\n\nprint('triples', triples, 'doubles', doubles, 'singles', singles)\nprint('sanity check: ', 'total faults', y.sum(), ' sum of above ', 3 * triples + 2 * doubles + singles)\n# plt.plot(meta_train.target)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"df30f3c42f11c39ab8ad134b238395d3b5ddbd7c"},"cell_type":"code","source":"output = pd.DataFrame({\"signal_id\":meta_test.signal_id[0:int(num_to_use * 3)]})\noutput[\"target\"] = pd.Series(y[:,0]) \noutput['signal_id'] = output['signal_id'].astype(np.int64)\noutput['target'] = output['target'].astype(np.int64)\noutput.to_csv(\"submission.csv\", index=False)\noutput","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}