{"cells":[{"metadata":{"_uuid":"e154a47bf09b8770980486e87786317a1b3038e1"},"cell_type":"markdown","source":"

VSB Power Line Fault Detection

\n
\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. Each 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

\nThe Challenge is to detect partial discharge patterns in signals acquired from these power lines with a new meter designed at the ENET Centre at VSB. Effective classifiers using this data will make it possible to continuously monitor power lines for faults.\n\n
Inspired from Bruno Aquino's Kernel:\nhttps://www.kaggle.com/braquino/5-fold-lstm-attention-fully-commented-0-694
\n\n
Feature Engineering : Raphael Vallat's entropy repository AND Capsule Layer : keras-contrib
\n\n
\nContent\n-  Import Library\n-  Define Matthews Correlation\n-  Create Attention Layer\n-  Feature Engineering\n-  Import Data\n-  Define Transformations\n-  Prepare Training Data\n-  Building the RNN Model Architecture\n-  Applying Stratified K Fold Technique\n-  Identify the Best Threshold\n-  Create Submission File\n
"},{"metadata":{"_uuid":"9882a9623b6a1b939d465a51ea381e56d9dc4f9c"},"cell_type":"markdown","source":"
Import Library
"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport pyarrow.parquet as pq # Used to read the data\nimport os \nimport numpy as np\nfrom keras.layers import * # Keras is the most friendly Neural Network library, this Kernel use a lot of layers classes\nfrom keras.models import Model\nfrom tqdm import tqdm # Processing time measurement\nfrom sklearn.model_selection import train_test_split \nfrom keras import backend as K # The backend give us access to tensorflow operations and allow us to create the Attention class\nfrom keras import optimizers # Allow us to access the Adam class to modify some parameters\nfrom sklearn.model_selection import GridSearchCV, StratifiedKFold,RepeatedStratifiedKFold # Used to use Kfold to train our model\nfrom keras.callbacks import * # This object helps the model to train in a smarter way, avoiding overfitting\nfrom sklearn.preprocessing import MinMaxScaler\nfrom numba import jit\nfrom math import log, floor\nfrom sklearn.neighbors import KDTree\nfrom scipy.signal import periodogram, welch\nimport tensorflow as tf","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6e6379386e44afc69bee8895a52da22199e888fb","trusted":true},"cell_type":"code","source":"# select how many folds will be created\nN_SPLITS = 5\n# it is just a constant with the measurements data size\nsample_size = 800000","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"38d0c7c3b97be5cd3e5b72135f5a7de6933de953"},"cell_type":"markdown","source":"
Define Matthews Correlation
"},{"metadata":{"_uuid":"c3340ee96becb5ca8f075d9c44b7df383ddba5ee","trusted":true},"cell_type":"code","source":"# It is the official metric used in this competition\n# below is the declaration of a function used inside the keras model, calculation with K (keras backend / thensorflow)\ndef matthews_correlation(y_true, y_pred):\n '''Calculates the Matthews correlation coefficient measure for quality\n of binary classification problems.\n '''\n y_pred = tf.convert_to_tensor(y_pred, np.float32)\n y_true = tf.convert_to_tensor(y_true, np.float32)\n \n y_pred_pos = K.round(K.clip(y_pred, 0, 1))\n y_pred_neg = 1 - y_pred_pos\n\n y_pos = K.round(K.clip(y_true, 0, 1))\n y_neg = 1 - y_pos\n\n tp = K.sum(y_pos * y_pred_pos)\n tn = K.sum(y_neg * y_pred_neg)\n\n fp = K.sum(y_neg * y_pred_pos)\n fn = K.sum(y_pos * y_pred_neg)\n\n numerator = (tp * tn - fp * fn)\n denominator = K.sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))\n\n return numerator / (denominator + K.epsilon())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"240ff24f0e5c2d976131af097f4becb883a11f95"},"cell_type":"markdown","source":"
Create Attention Layer
"},{"metadata":{"_uuid":"eda7ea366117d1ce8e5fce69e5bba333821d8b48","trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/suicaokhoailang/lstm-attention-baseline-0-652-lb\n\nclass Attention(Layer):\n def __init__(self, step_dim,\n W_regularizer=None, b_regularizer=None,\n W_constraint=None, b_constraint=None,\n bias=True, **kwargs):\n self.supports_masking = True\n self.init = initializers.get('glorot_uniform')\n\n self.W_regularizer = regularizers.get(W_regularizer)\n self.b_regularizer = regularizers.get(b_regularizer)\n\n self.W_constraint = constraints.get(W_constraint)\n self.b_constraint = constraints.get(b_constraint)\n\n self.bias = bias\n self.step_dim = step_dim\n self.features_dim = 0\n super(Attention, self).__init__(**kwargs)\n\n def build(self, input_shape):\n assert len(input_shape) == 3\n\n self.W = self.add_weight((input_shape[-1],),\n initializer=self.init,\n name='{}_W'.format(self.name),\n regularizer=self.W_regularizer,\n constraint=self.W_constraint)\n self.features_dim = input_shape[-1]\n\n if self.bias:\n self.b = self.add_weight((input_shape[1],),\n initializer='zero',\n name='{}_b'.format(self.name),\n regularizer=self.b_regularizer,\n constraint=self.b_constraint)\n else:\n self.b = None\n\n self.built = True\n\n def compute_mask(self, input, input_mask=None):\n return None\n\n def call(self, x, mask=None):\n features_dim = self.features_dim\n step_dim = self.step_dim\n\n eij = K.reshape(K.dot(K.reshape(x, (-1, features_dim)),\n K.reshape(self.W, (features_dim, 1))), (-1, step_dim))\n\n if self.bias:\n eij += self.b\n\n eij = K.tanh(eij)\n\n a = K.exp(eij)\n\n if mask is not None:\n a *= K.cast(mask, K.floatx())\n\n a /= K.cast(K.sum(a, axis=1, keepdims=True) + K.epsilon(), K.floatx())\n\n a = K.expand_dims(a)\n weighted_input = x * a\n return K.sum(weighted_input, axis=1)\n\n def compute_output_shape(self, input_shape):\n return input_shape[0], self.features_dim","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8217f085a654d709b625394d9c6aec586d5fb6e9"},"cell_type":"markdown","source":"
Feature Engineering
"},{"metadata":{"trusted":true,"_uuid":"518408ef5066aa933142fdf9cc221ead5aa8b51b"},"cell_type":"code","source":"def squash(x, axis=-1):\n \"\"\"\n Squash activation function (generally used in Capsule layers).\n \"\"\"\n s_squared_norm = K.sum(K.square(x), axis, keepdims=True) + K.epsilon()\n scale = K.sqrt(s_squared_norm) / (0.5 + s_squared_norm)\n return scale * x\n\n\nclass Capsule(Layer):\n \"\"\"Capsule Layer implementation in Keras\n This implementation is based on Dynamic Routing of Capsules,\n Geoffrey Hinton et. al.\n The Capsule Layer is a Neural Network Layer which helps\n modeling relationships in image and sequential data better\n than just CNNs or RNNs. It achieves this by understanding\n the spatial relationships between objects (in images)\n or words (in text) by encoding additional information\n about the image or text, such as angle of rotation,\n thickness and brightness, relative proportions etc.\n This layer can be used instead of pooling layers to\n lower dimensions and still capture important information\n about the relationships and structures within the data.\n A normal pooling layer would lose a lot of\n this information.\n This layer can be used on the output of any layer\n which has a 3-D output (including batch_size). For example,\n in image classification, it can be used on the output of a\n Conv2D layer for Computer Vision applications. Also,\n it can be used on the output of a GRU or LSTM Layer\n (Bidirectional or Unidirectional) for NLP applications.\n The default activation function is 'linear'. But, this layer\n is generally used with the 'squash' activation function\n (recommended). To use the squash activation function, do :\n from keras_contrib.activations import squash\n capsule = Capsule(num_capsule=10,\n dim_capsule=10,\n routings=3,\n share_weights=True,\n activation=squash)\n # Example usage :\n 1). COMPUTER VISION\n input_image = Input(shape=(None, None, 3))\n conv_2d = Conv2D(64,\n (3, 3),\n activation='relu')(input_image)\n capsule = Capsule(num_capsule=10,\n dim_capsule=16,\n routings=3,\n activation='relu',\n share_weights=True)(conv_2d)\n 2). NLP\n maxlen = 72\n max_features = 120000\n input_text = Input(shape=(maxlen,))\n embedding = Embedding(max_features,\n embed_size,\n weights=[embedding_matrix],\n trainable=False)(input_text)\n bi_gru = Bidirectional(GRU(64,\n return_seqeunces=True))(embedding)\n capsule = Capsule(num_capsule=5,\n dim_capsule=5,\n routings=4,\n activation='sigmoid',\n share_weights=True)(bi_gru)\n # Arguments\n num_capsule : Number of Capsules (int)\n dim_capsules : Dimensions of the vector output of each Capsule (int)\n routings : Number of dynamic routings in the Capsule Layer (int)\n share_weights : Whether to share weights between Capsules or not\n (boolean)\n activation : Activation function for the Capsules\n regularizer : Regularizer for the weights of the Capsules\n initializer : Initializer for the weights of the Caspules\n constraint : Constraint for the weights of the Capsules\n # Input shape\n 3D tensor with shape:\n (batch_size, input_num_capsule, input_dim_capsule)\n [any 3-D Tensor with the first dimension as batch_size]\n # Output shape\n 3D tensor with shape:\n (batch_size, num_capsule, dim_capsule)\n # References\n - [Dynamic-Routing-Between-Capsules]\n (https://arxiv.org/pdf/1710.09829.pdf)\n - [Keras-Examples-CIFAR10-CNN-Capsule]\"\"\"\n\n def __init__(self,\n num_capsule,\n dim_capsule,\n routings=3,\n share_weights=True,\n initializer='glorot_uniform',\n activation=squash,\n regularizer=None,\n constraint=None,\n **kwargs):\n super(Capsule, self).__init__(**kwargs)\n self.num_capsule = num_capsule\n self.dim_capsule = dim_capsule\n self.routings = routings\n self.share_weights = share_weights\n\n self.activation = activations.get(activation)\n self.regularizer = regularizers.get(regularizer)\n self.initializer = initializers.get(initializer)\n self.constraint = constraints.get(constraint)\n\n def build(self, input_shape):\n input_dim_capsule = input_shape[-1]\n if self.share_weights:\n self.W = self.add_weight(name='capsule_kernel',\n shape=(1,\n input_dim_capsule,\n self.num_capsule *\n self.dim_capsule),\n initializer=self.initializer,\n regularizer=self.regularizer,\n constraint=self.constraint,\n trainable=True)\n else:\n input_num_capsule = input_shape[-2]\n self.W = self.add_weight(name='capsule_kernel',\n shape=(input_num_capsule,\n input_dim_capsule,\n self.num_capsule *\n self.dim_capsule),\n initializer=self.initializer,\n regularizer=self.regularizer,\n constraint=self.constraint,\n trainable=True)\n\n self.build = True\n\n def call(self, inputs):\n if self.share_weights:\n u_hat_vectors = K.conv1d(inputs, self.W)\n else:\n u_hat_vectors = K.local_conv1d(inputs, self.W, [1], [1])\n\n # u_hat_vectors : The spatially transformed input vectors (with local_conv_1d)\n\n batch_size = K.shape(inputs)[0]\n input_num_capsule = K.shape(inputs)[1]\n u_hat_vectors = K.reshape(u_hat_vectors, (batch_size,\n input_num_capsule,\n self.num_capsule,\n self.dim_capsule))\n\n u_hat_vectors = K.permute_dimensions(u_hat_vectors, (0, 2, 1, 3))\n routing_weights = K.zeros_like(u_hat_vectors[:, :, :, 0])\n\n for i in range(self.routings):\n capsule_weights = K.softmax(routing_weights, 1)\n outputs = K.batch_dot(capsule_weights, u_hat_vectors, [2, 2])\n if K.ndim(outputs) == 4:\n outputs = K.sum(outputs, axis=1)\n if i < self.routings - 1:\n outputs = K.l2_normalize(outputs, -1)\n routing_weights = K.batch_dot(outputs, u_hat_vectors, [2, 3])\n if K.ndim(routing_weights) == 4:\n routing_weights = K.sum(routing_weights, axis=1)\n\n return self.activation(outputs)\n\n def compute_output_shape(self, input_shape):\n return (None, self.num_capsule, self.dim_capsule)\n\n def get_config(self):\n config = {'num_capsule': self.num_capsule,\n 'dim_capsule': self.dim_capsule,\n 'routings': self.routings,\n 'share_weights': self.share_weights,\n 'activation': activations.serialize(self.activation),\n 'regularizer': regularizers.serialize(self.regularizer),\n 'initializer': initializers.serialize(self.initializer),\n 'constraint': constraints.serialize(self.constraint)}\n\n base_config = super(Capsule, self).get_config()\n return dict(list(base_config.items()) + list(config.items()))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0c94824ab8a2f0e7ee38811a6effecd7e2b6477d"},"cell_type":"markdown","source":"
Import Data
"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# just load train data\ndf_train = pd.read_csv('../input/metadata_train.csv')\n# set index, it makes the data access much faster\ndf_train = df_train.set_index(['id_measurement', 'phase'])\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"26df6c7fbfecd537404866faec13d1238ae3ebc6","trusted":true},"cell_type":"code","source":"# in other notebook I have extracted the min and max values from the train data, the measurements\nmax_num = 127\nmin_num = -128","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"87d74c8547a580e74dd945600b361a96ce8aa0d8"},"cell_type":"markdown","source":"
Define Transformations
"},{"metadata":{"_uuid":"7b0717b14bcfcba1f48d33c8161ae51c778687af","trusted":true},"cell_type":"code","source":"# This function standardize the data from (-128 to 127) to (-1 to 1)\n# Theoretically it helps in the NN Model training, but I didn't tested without it\ndef min_max_transformation(ts, min_data, max_data, range_needed=(-1,1)):\n if min_data < 0:\n ts_std = (ts + abs(min_data)) / (max_data + abs(min_data))\n else:\n ts_std = (ts - min_data) / (max_data - min_data)\n if range_needed[0] < 0: \n return ts_std * (range_needed[1] + abs(range_needed[0])) + range_needed[0]\n else:\n return ts_std * (range_needed[1] - range_needed[0]) + range_needed[0]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c6137bbbe75c3a1509a5f98e08805dbbd492aa37","trusted":true},"cell_type":"code","source":"def transform_ts(ts, n_dim=160, min_max=(-1,1)):\n # convert data into -1 to 1\n ts_std = min_max_transformation(ts, min_data=min_num, max_data=max_num)\n # bucket or chunk size, 5000 in this case (800000 / 160)\n bucket_size = int(sample_size / n_dim)\n # new_ts will be the container of the new data\n new_ts = []\n # this for iteract any chunk/bucket until reach the whole sample_size (800000)\n for i in range(0, sample_size, bucket_size):\n # cut each bucket to ts_range\n ts_range = ts_std[i:i + bucket_size]\n # calculate each feature\n mean = ts_range.mean()\n std = ts_range.std() # standard deviation\n std_top = mean + std # I have to test it more, but is is like a band\n std_bot = mean - std\n # I think that the percentiles are very important, it is like a distribuiton analysis from eath chunk\n percentil_calc = np.percentile(ts_range, [0, 1, 25, 50, 75, 99, 100]) \n max_range = percentil_calc[-1] - percentil_calc[0] # this is the amplitude of the chunk\n relative_percentile = percentil_calc - mean # maybe it could heap to understand the asymmetry\n new_ts.append(np.concatenate([np.asarray([mean, std, std_top, std_bot, max_range]),percentil_calc, relative_percentile]))\n return np.asarray(new_ts)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"48f1bb0248e8ec5a8e20e48061f71692491d2c5e"},"cell_type":"markdown","source":"
Prepare Training Data
"},{"metadata":{"_uuid":"7460e718a605803f1d9e4fbec61750a0deb02a47","trusted":true},"cell_type":"code","source":"# this function take a piece of data and convert using transform_ts(), but it does to each of the 3 phases\n# if we would try to do in one time, could exceed the RAM Memmory\ndef prepare_data(start, end):\n\n praq_train = pq.read_pandas('../input/train.parquet', columns=[str(i) for i in range(start, end)]).to_pandas()\n X = []\n y = []\n \n for id_measurement in tqdm(df_train.index.levels[0].unique()[int(start/3):int(end/3)]):\n X_signal = []\n # for each phase of the signal\n for phase in [0,1,2]:\n # extract from df_train both signal_id and target to compose the new data sets\n signal_id, target = df_train.loc[id_measurement].loc[phase]\n # but just append the target one time, to not triplicate it\n if phase == 0:\n y.append(target)\n # extract and transform data into sets of features\n X_signal.append(transform_ts(praq_train[str(signal_id)]))\n # concatenate all the 3 phases in one matrix\n X_signal = np.concatenate(X_signal, axis=1)\n # add the data to X\n X.append(X_signal)\n X = np.asarray(X)\n y = np.asarray(y)\n return X, y\n\n#Reference: https://www.kaggle.com/jiweiliu/lgb-2-leaves-augment?scriptVersionId=11603471\ndef augment(x,y,t=2):\n xs,xn = [],[]\n for i in range(t):\n mask = y>0\n x1 = x[mask].copy()\n ids = np.arange(x1.shape[0])\n for c in range(x1.shape[1]):\n np.random.shuffle(ids)\n x1[:,c] = x1[ids][:,c]\n xs.append(x1)\n\n for i in range(t//2):\n mask = y==0\n x1 = x[mask].copy()\n ids = np.arange(x1.shape[0])\n for c in range(x1.shape[1]):\n np.random.shuffle(ids)\n x1[:,c] = x1[ids][:,c]\n xn.append(x1)\n\n xs = np.vstack(xs)\n xn = np.vstack(xn)\n ys = np.ones(xs.shape[0])\n yn = np.zeros(xn.shape[0])\n x = np.vstack([x,xs,xn])\n y = np.concatenate([y,ys,yn])\n return x,y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"52dc826ab9ee1dd56c9fb29bd5c1b2d26b5928bf","trusted":true},"cell_type":"code","source":"# this code is very simple, divide the total size of the df_train into two sets and process it\nX = []\ny = []\ndef load_all():\n total_size = len(df_train)\n for ini, end in [(0, int(total_size/2)), (int(total_size/2), total_size)]:\n X_temp, y_temp = prepare_data(ini, end)\n X_aug, y_aug = augment(X_temp,y_temp,t=2)\n X.append(X_aug)\n y.append(y_aug)\nload_all()\nX = np.concatenate(X)\ny = np.concatenate(y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9134417664c4713f75e7b427891fa0041ac70c98"},"cell_type":"code","source":"def petrosian_fd(x):\n \"\"\"Petrosian fractal dimension.\n Parameters\n ----------\n x : list or np.array\n One dimensional time series\n Returns\n -------\n pfd : float\n Petrosian fractal dimension\n Notes\n -----\n The Petrosian algorithm can be used to provide a fast computation of\n the FD of a signal by translating the series into a binary sequence.\n The Petrosian fractal dimension of a time series :math:`x` is defined by:\n .. math:: \\dfrac{log_{10}(N)}{log_{10}(N) +\n log_{10}(\\dfrac{N}{N+0.4N_{\\Delta}})}\n where :math:`N` is the length of the time series, and\n :math:`N_{\\Delta}` is the number of sign changes in the binary sequence.\n Original code from the pyrem package by Quentin Geissmann.\n References\n ----------\n .. [1] A. Petrosian, Kolmogorov complexity of finite sequences and\n recognition of different preictal EEG patterns, in , Proceedings of the\n Eighth IEEE Symposium on Computer-Based Medical Systems, 1995,\n pp. 212-217.\n .. [2] Goh, Cindy, et al. \"Comparison of fractal dimension algorithms for\n the computation of EEG biomarkers for dementia.\" 2nd International\n Conference on Computational Intelligence in Medicine and Healthcare\n (CIMED2005). 2005.\n Examples\n --------\n Petrosian fractal dimension.\n >>> import numpy as np\n >>> from entropy import petrosian_fd\n >>> np.random.seed(123)\n >>> x = np.random.rand(100)\n >>> print(petrosian_fd(x))\n 1.0505\n \"\"\"\n n = len(x)\n # Number of sign changes in the first derivative of the signal\n diff = np.ediff1d(x)\n N_delta = (diff[1:-1] * diff[0:-2] < 0).sum()\n return np.log10(n) / (np.log10(n) + np.log10(n / (n + 0.4 * N_delta)))\n\n\ndef katz_fd(x):\n \"\"\"Katz Fractal Dimension.\n Parameters\n ----------\n x : list or np.array\n One dimensional time series\n Returns\n -------\n kfd : float\n Katz fractal dimension\n Notes\n -----\n The Katz Fractal dimension is defined by:\n .. math:: FD_{Katz} = \\dfrac{log_{10}(n)}{log_{10}(d/L)+log_{10}(n)}\n where :math:`L` is the total length of the time series and :math:`d`\n is the Euclidean distance between the first point in the\n series and the point that provides the furthest distance\n with respect to the first point.\n Original code from the mne-features package by Jean-Baptiste Schiratti\n and Alexandre Gramfort.\n References\n ----------\n .. [1] Esteller, R. et al. (2001). A comparison of waveform fractal\n dimension algorithms. IEEE Transactions on Circuits and Systems I:\n Fundamental Theory and Applications, 48(2), 177-183.\n .. [2] Goh, Cindy, et al. \"Comparison of fractal dimension algorithms for\n the computation of EEG biomarkers for dementia.\" 2nd International\n Conference on Computational Intelligence in Medicine and Healthcare\n (CIMED2005). 2005.\n Examples\n --------\n Katz fractal dimension.\n >>> import numpy as np\n >>> from entropy import katz_fd\n >>> np.random.seed(123)\n >>> x = np.random.rand(100)\n >>> print(katz_fd(x))\n 5.1214\n \"\"\"\n x = np.array(x)\n dists = np.abs(np.ediff1d(x))\n ll = dists.sum()\n ln = np.log10(np.divide(ll, dists.mean()))\n aux_d = x - x[0]\n d = np.max(np.abs(aux_d[1:]))\n return np.divide(ln, np.add(ln, np.log10(np.divide(d, ll))))\n\n@jit('UniTuple(float64, 2)(float64[:], float64[:])', nopython=True)\ndef _linear_regression(x, y):\n \"\"\"Fast linear regression using Numba.\n Parameters\n ----------\n x, y : ndarray, shape (n_times,)\n Variables\n Returns\n -------\n slope : float\n Slope of 1D least-square regression.\n intercept : float\n Intercept\n \"\"\"\n n_times = x.size\n sx2 = 0\n sx = 0\n sy = 0\n sxy = 0\n for j in range(n_times):\n sx2 += x[j] ** 2\n sx += x[j]\n sxy += x[j] * y[j]\n sy += y[j]\n den = n_times * sx2 - (sx ** 2)\n num = n_times * sxy - sx * sy\n slope = num / den\n intercept = np.mean(y) - slope * np.mean(x)\n return slope, intercept\n\n\n@jit('i8[:](f8, f8, f8)', nopython=True)\ndef _log_n(min_n, max_n, factor):\n \"\"\"\n Creates a list of integer values by successively multiplying a minimum\n value min_n by a factor > 1 until a maximum value max_n is reached.\n Used for detrended fluctuation analysis (DFA).\n Function taken from the nolds python package\n (https://github.com/CSchoel/nolds) by Christopher Scholzel.\n Parameters\n ----------\n min_n (float):\n minimum value (must be < max_n)\n max_n (float):\n maximum value (must be > min_n)\n factor (float):\n factor used to increase min_n (must be > 1)\n Returns\n -------\n list of integers:\n min_n, min_n * factor, min_n * factor^2, ... min_n * factor^i < max_n\n without duplicates\n \"\"\"\n max_i = int(floor(log(1.0 * max_n / min_n) / log(factor)))\n ns = [min_n]\n for i in range(max_i + 1):\n n = int(floor(min_n * (factor ** i)))\n if n > ns[-1]:\n ns.append(n)\n return np.array(ns, dtype=np.int64)\n\n@jit('float64(float64[:], int32)')\ndef _higuchi_fd(x, kmax):\n \"\"\"Utility function for `higuchi_fd`.\n \"\"\"\n n_times = x.size\n lk = np.empty(kmax)\n x_reg = np.empty(kmax)\n y_reg = np.empty(kmax)\n for k in range(1, kmax + 1):\n lm = np.empty((k,))\n for m in range(k):\n ll = 0\n n_max = floor((n_times - m - 1) / k)\n n_max = int(n_max)\n for j in range(1, n_max):\n ll += abs(x[m + j * k] - x[m + (j - 1) * k])\n ll /= k\n ll *= (n_times - 1) / (k * n_max)\n lm[m] = ll\n # Mean of lm\n m_lm = 0\n for m in range(k):\n m_lm += lm[m]\n m_lm /= k\n lk[k - 1] = m_lm\n x_reg[k - 1] = log(1. / k)\n y_reg[k - 1] = log(m_lm)\n higuchi, _ = _linear_regression(x_reg, y_reg)\n return higuchi\n\n\ndef higuchi_fd(x, kmax=10):\n \"\"\"Higuchi Fractal Dimension.\n Parameters\n ----------\n x : list or np.array\n One dimensional time series\n kmax : int\n Maximum delay/offset (in number of samples).\n Returns\n -------\n hfd : float\n Higuchi Fractal Dimension\n Notes\n -----\n Original code from the mne-features package by Jean-Baptiste Schiratti\n and Alexandre Gramfort.\n The `higuchi_fd` function uses Numba to speed up the computation.\n References\n ----------\n .. [1] Higuchi, Tomoyuki. \"Approach to an irregular time series on the\n basis of the fractal theory.\" Physica D: Nonlinear Phenomena 31.2\n (1988): 277-283.\n Examples\n --------\n Higuchi Fractal Dimension\n >>> import numpy as np\n >>> from entropy import higuchi_fd\n >>> np.random.seed(123)\n >>> x = np.random.rand(100)\n >>> print(higuchi_fd(x))\n 2.051179\n \"\"\"\n x = np.asarray(x, dtype=np.float64)\n kmax = int(kmax)\n return _higuchi_fd(x, kmax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"282a1118464ab49b86f3392a41d42e60a6853924"},"cell_type":"code","source":"def _embed(x, order=3, delay=1):\n \"\"\"Time-delay embedding.\n Parameters\n ----------\n x : 1d-array, shape (n_times)\n Time series\n order : int\n Embedding dimension (order)\n delay : int\n Delay.\n Returns\n -------\n embedded : ndarray, shape (n_times - (order - 1) * delay, order)\n Embedded time-series.\n \"\"\"\n N = len(x)\n if order * delay > N:\n raise ValueError(\"Error: order * delay should be lower than x.size\")\n if delay < 1:\n raise ValueError(\"Delay has to be at least 1.\")\n if order < 2:\n raise ValueError(\"Order has to be at least 2.\")\n Y = np.zeros((order, N - (order - 1) * delay))\n for i in range(order):\n Y[i] = x[i * delay:i * delay + Y.shape[1]]\n return Y.T\n\nall = ['perm_entropy', 'spectral_entropy', 'svd_entropy', 'app_entropy',\n 'sample_entropy']\n\n\ndef perm_entropy(x, order=3, delay=1, normalize=False):\n \"\"\"Permutation Entropy.\n Parameters\n ----------\n x : list or np.array\n One-dimensional time series of shape (n_times)\n order : int\n Order of permutation entropy\n delay : int\n Time delay\n normalize : bool\n If True, divide by log2(order!) to normalize the entropy between 0\n and 1. Otherwise, return the permutation entropy in bit.\n Returns\n -------\n pe : float\n Permutation Entropy\n Notes\n -----\n The permutation entropy is a complexity measure for time-series first\n introduced by Bandt and Pompe in 2002 [1]_.\n The permutation entropy of a signal :math:`x` is defined as:\n .. math:: H = -\\sum p(\\pi)log_2(\\pi)\n where the sum runs over all :math:`n!` permutations :math:`\\pi` of order\n :math:`n`. This is the information contained in comparing :math:`n`\n consecutive values of the time series. It is clear that\n :math:`0 ≤ H (n) ≤ log_2(n!)` where the lower bound is attained for an\n increasing or decreasing sequence of values, and the upper bound for a\n completely random system where all :math:`n!` possible permutations appear\n with the same probability.\n The embedded matrix :math:`Y` is created by:\n .. math:: y(i)=[x_i,x_{i+delay}, ...,x_{i+(order-1) * delay}]\n .. math:: Y=[y(1),y(2),...,y(N-(order-1))*delay)]^T\n References\n ----------\n .. [1] Bandt, Christoph, and Bernd Pompe. \"Permutation entropy: a\n natural complexity measure for time series.\" Physical review letters\n 88.17 (2002): 174102.\n Examples\n --------\n 1. Permutation entropy with order 2\n >>> from entropy import perm_entropy\n >>> x = [4, 7, 9, 10, 6, 11, 3]\n >>> # Return a value in bit between 0 and log2(factorial(order))\n >>> print(perm_entropy(x, order=2))\n 0.918\n 2. Normalized permutation entropy with order 3\n >>> from entropy import perm_entropy\n >>> x = [4, 7, 9, 10, 6, 11, 3]\n >>> # Return a value comprised between 0 and 1.\n >>> print(perm_entropy(x, order=3, normalize=True))\n 0.589\n \"\"\"\n x = np.array(x)\n ran_order = range(order)\n hashmult = np.power(order, ran_order)\n # Embed x and sort the order of permutations\n sorted_idx = _embed(x, order=order, delay=delay).argsort(kind='quicksort')\n # Associate unique integer to each permutations\n hashval = (np.multiply(sorted_idx, hashmult)).sum(1)\n # Return the counts\n _, c = np.unique(hashval, return_counts=True)\n # Use np.true_divide for Python 2 compatibility\n p = np.true_divide(c, c.sum())\n pe = -np.multiply(p, np.log2(p)).sum()\n if normalize:\n pe /= np.log2(factorial(order))\n return pe\n\n\ndef spectral_entropy(x, sf, method='fft', nperseg=None, normalize=False):\n \"\"\"Spectral Entropy.\n Parameters\n ----------\n x : list or np.array\n One-dimensional time series of shape (n_times)\n sf : float\n Sampling frequency\n method : str\n Spectral estimation method ::\n 'fft' : Fourier Transform (via scipy.signal.periodogram)\n 'welch' : Welch periodogram (via scipy.signal.welch)\n nperseg : str or int\n Length of each FFT segment for Welch method.\n If None, uses scipy default of 256 samples.\n normalize : bool\n If True, divide by log2(psd.size) to normalize the spectral entropy\n between 0 and 1. Otherwise, return the spectral entropy in bit.\n Returns\n -------\n se : float\n Spectral Entropy\n Notes\n -----\n Spectral Entropy is defined to be the Shannon Entropy of the Power\n Spectral Density (PSD) of the data:\n .. math:: H(x, sf) = -\\sum_{f=0}^{f_s/2} PSD(f) log_2[PSD(f)]\n Where :math:`PSD` is the normalised PSD, and :math:`f_s` is the sampling\n frequency.\n References\n ----------\n .. [1] Inouye, T. et al. (1991). Quantification of EEG irregularity by\n use of the entropy of the power spectrum. Electroencephalography\n and clinical neurophysiology, 79(3), 204-210.\n Examples\n --------\n 1. Spectral entropy of a pure sine using FFT\n >>> from entropy import spectral_entropy\n >>> import numpy as np\n >>> sf, f, dur = 100, 1, 4\n >>> N = sf * duration # Total number of discrete samples\n >>> t = np.arange(N) / sf # Time vector\n >>> x = np.sin(2 * np.pi * f * t)\n >>> print(np.round(spectral_entropy(x, sf, method='fft'), 2)\n 0.0\n 2. Spectral entropy of a random signal using Welch's method\n >>> from entropy import spectral_entropy\n >>> import numpy as np\n >>> np.random.seed(42)\n >>> x = np.random.rand(3000)\n >>> print(spectral_entropy(x, sf=100, method='welch'))\n 9.939\n 3. Normalized spectral entropy\n >>> print(spectral_entropy(x, sf=100, method='welch', normalize=True))\n 0.995\n \"\"\"\n x = np.array(x)\n # Compute and normalize power spectrum\n if method == 'fft':\n _, psd = periodogram(x, sf)\n elif method == 'welch':\n _, psd = welch(x, sf, nperseg=nperseg)\n psd_norm = np.divide(psd, psd.sum())\n se = -np.multiply(psd_norm, np.log2(psd_norm)).sum()\n if normalize:\n se /= np.log2(psd_norm.size)\n return se\n\n\ndef svd_entropy(x, order=3, delay=1, normalize=False):\n \"\"\"Singular Value Decomposition entropy.\n Parameters\n ----------\n x : list or np.array\n One-dimensional time series of shape (n_times)\n order : int\n Order of permutation entropy\n delay : int\n Time delay\n normalize : bool\n If True, divide by log2(order!) to normalize the entropy between 0\n and 1. Otherwise, return the permutation entropy in bit.\n Returns\n -------\n svd_e : float\n SVD Entropy\n Notes\n -----\n SVD entropy is an indicator of the number of eigenvectors that are needed\n for an adequate explanation of the data set. In other words, it measures\n the dimensionality of the data.\n The SVD entropy of a signal :math:`x` is defined as:\n .. math::\n H = -\\sum_{i=1}^{M} \\overline{\\sigma}_i log_2(\\overline{\\sigma}_i)\n where :math:`M` is the number of singular values of the embedded matrix\n :math:`Y` and :math:`\\sigma_1, \\sigma_2, ..., \\sigma_M` are the\n normalized singular values of :math:`Y`.\n The embedded matrix :math:`Y` is created by:\n .. math:: y(i)=[x_i,x_{i+delay}, ...,x_{i+(order-1) * delay}]\n .. math:: Y=[y(1),y(2),...,y(N-(order-1))*delay)]^T\n Examples\n --------\n 1. SVD entropy with order 2\n >>> from entropy import svd_entropy\n >>> x = [4, 7, 9, 10, 6, 11, 3]\n >>> # Return a value in bit between 0 and log2(factorial(order))\n >>> print(svd_entropy(x, order=2))\n 0.762\n 2. Normalized SVD entropy with order 3\n >>> from entropy import svd_entropy\n >>> x = [4, 7, 9, 10, 6, 11, 3]\n >>> # Return a value comprised between 0 and 1.\n >>> print(svd_entropy(x, order=3, normalize=True))\n 0.687\n \"\"\"\n x = np.array(x)\n mat = _embed(x, order=order, delay=delay)\n W = np.linalg.svd(mat, compute_uv=False)\n # Normalize the singular values\n W /= sum(W)\n svd_e = -np.multiply(W, np.log2(W)).sum()\n if normalize:\n svd_e /= np.log2(order)\n return svd_e\n\n\ndef _app_samp_entropy(x, order, metric='chebyshev', approximate=True):\n \"\"\"Utility function for `app_entropy`` and `sample_entropy`.\n \"\"\"\n _all_metrics = KDTree.valid_metrics\n if metric not in _all_metrics:\n raise ValueError('The given metric (%s) is not valid. The valid '\n 'metric names are: %s' % (metric, _all_metrics))\n phi = np.zeros(2)\n r = 0.2 * np.std(x, axis=-1, ddof=1)\n\n # compute phi(order, r)\n _emb_data1 = _embed(x, order, 1)\n if approximate:\n emb_data1 = _emb_data1\n else:\n emb_data1 = _emb_data1[:-1]\n count1 = KDTree(emb_data1, metric=metric).query_radius(emb_data1, r,\n count_only=True\n ).astype(np.float64)\n # compute phi(order + 1, r)\n emb_data2 = _embed(x, order + 1, 1)\n count2 = KDTree(emb_data2, metric=metric).query_radius(emb_data2, r,\n count_only=True\n ).astype(np.float64)\n if approximate:\n phi[0] = np.mean(np.log(count1 / emb_data1.shape[0]))\n phi[1] = np.mean(np.log(count2 / emb_data2.shape[0]))\n else:\n phi[0] = np.mean((count1 - 1) / (emb_data1.shape[0] - 1))\n phi[1] = np.mean((count2 - 1) / (emb_data2.shape[0] - 1))\n return phi\n\n\n@jit('f8(f8[:], i4, f8)', nopython=True)\ndef _numba_sampen(x, mm=2, r=0.2):\n \"\"\"\n Fast evaluation of the sample entropy using Numba.\n \"\"\"\n n = x.size\n n1 = n - 1\n mm += 1\n mm_dbld = 2 * mm\n\n # Define threshold\n r *= x.std()\n\n # initialize the lists\n run = [0] * n\n run1 = run[:]\n r1 = [0] * (n * mm_dbld)\n a = [0] * mm\n b = a[:]\n p = a[:]\n\n for i in range(n1):\n nj = n1 - i\n\n for jj in range(nj):\n j = jj + i + 1\n if abs(x[j] - x[i]) < r:\n run[jj] = run1[jj] + 1\n m1 = mm if mm < run[jj] else run[jj]\n for m in range(m1):\n a[m] += 1\n if j < n1:\n b[m] += 1\n else:\n run[jj] = 0\n for j in range(mm_dbld):\n run1[j] = run[j]\n r1[i + n * j] = run[j]\n if nj > mm_dbld - 1:\n for j in range(mm_dbld, nj):\n run1[j] = run[j]\n\n m = mm - 1\n\n while m > 0:\n b[m] = b[m - 1]\n m -= 1\n\n b[0] = n * n1 / 2\n a = np.array([float(aa) for aa in a])\n b = np.array([float(bb) for bb in b])\n p = np.true_divide(a, b)\n return -log(p[-1])\n\n\ndef app_entropy(x, order=2, metric='chebyshev'):\n \"\"\"Approximate Entropy.\n Parameters\n ----------\n x : list or np.array\n One-dimensional time series of shape (n_times)\n order : int (default: 2)\n Embedding dimension.\n metric : str (default: chebyshev)\n Name of the metric function used with\n :class:`~sklearn.neighbors.KDTree`. The list of available\n metric functions is given by: ``KDTree.valid_metrics``.\n Returns\n -------\n ae : float\n Approximate Entropy.\n Notes\n -----\n Original code from the mne-features package.\n Approximate entropy is a technique used to quantify the amount of\n regularity and the unpredictability of fluctuations over time-series data.\n Smaller values indicates that the data is more regular and predictable.\n The value of :math:`r` is set to :math:`0.2 * std(x)`.\n Code adapted from the mne-features package by Jean-Baptiste Schiratti\n and Alexandre Gramfort.\n References\n ----------\n .. [1] Richman, J. S. et al. (2000). Physiological time-series analysis\n using approximate entropy and sample entropy. American Journal of\n Physiology-Heart and Circulatory Physiology, 278(6), H2039-H2049.\n 1. Approximate entropy with order 2.\n >>> from entropy import app_entropy\n >>> import numpy as np\n >>> np.random.seed(1234567)\n >>> x = np.random.rand(3000)\n >>> print(app_entropy(x, order=2))\n 2.075\n \"\"\"\n phi = _app_samp_entropy(x, order=order, metric=metric, approximate=True)\n return np.subtract(phi[0], phi[1])\n\n\ndef sample_entropy(x, order=2, metric='chebyshev'):\n \"\"\"Sample Entropy.\n Parameters\n ----------\n x : list or np.array\n One-dimensional time series of shape (n_times)\n order : int (default: 2)\n Embedding dimension.\n metric : str (default: chebyshev)\n Name of the metric function used with KDTree. The list of available\n metric functions is given by: `KDTree.valid_metrics`.\n Returns\n -------\n se : float\n Sample Entropy.\n Notes\n -----\n Sample entropy is a modification of approximate entropy, used for assessing\n the complexity of physiological time-series signals. It has two advantages\n over approximate entropy: data length independence and a relatively\n trouble-free implementation. Large values indicate high complexity whereas\n smaller values characterize more self-similar and regular signals.\n Sample entropy of a signal :math:`x` is defined as:\n .. math:: H(x, m, r) = -log\\dfrac{C(m + 1, r)}{C(m, r)}\n where :math:`m` is the embedding dimension (= order), :math:`r` is\n the radius of the neighbourhood (default = :math:`0.2 * std(x)`),\n :math:`C(m + 1, r)` is the number of embedded vectors of length\n :math:`m + 1` having a Chebyshev distance inferior to :math:`r` and\n :math:`C(m, r)` is the number of embedded vectors of length\n :math:`m` having a Chebyshev distance inferior to :math:`r`.\n Note that if metric == 'chebyshev' and x.size < 5000 points, then the\n sample entropy is computed using a fast custom Numba script. For other\n metric types or longer time-series, the sample entropy is computed using\n a code from the mne-features package by Jean-Baptiste Schiratti\n and Alexandre Gramfort (requires sklearn).\n References\n ----------\n .. [1] Richman, J. S. et al. (2000). Physiological time-series analysis\n using approximate entropy and sample entropy. American Journal of\n Physiology-Heart and Circulatory Physiology, 278(6), H2039-H2049.\n Examples\n --------\n 1. Sample entropy with order 2.\n >>> from entropy import sample_entropy\n >>> import numpy as np\n >>> np.random.seed(1234567)\n >>> x = np.random.rand(3000)\n >>> print(sample_entropy(x, order=2))\n 2.192\n 2. Sample entropy with order 3 using the Euclidean distance.\n >>> from entropy import sample_entropy\n >>> import numpy as np\n >>> np.random.seed(1234567)\n >>> x = np.random.rand(3000)\n >>> print(sample_entropy(x, order=3, metric='euclidean'))\n 2.725\n \"\"\"\n x = np.asarray(x, dtype=np.float64)\n if metric == 'chebyshev' and x.size < 5000:\n return _numba_sampen(x, mm=order, r=0.2)\n else:\n phi = _app_samp_entropy(x, order=order, metric=metric,\n approximate=False)\n return -np.log(np.divide(phi[1], phi[0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"125197f288f8ec0bbeea4bd9702f3aeb1ef5dc8c"},"cell_type":"code","source":"def entropy_and_fractal_dim(x):\n return [perm_entropy(x), svd_entropy(x), app_entropy(x), sample_entropy(x), petrosian_fd(x), katz_fd(x), higuchi_fd(x)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8db969dcc7b41a65ddbe3fe9c001ee12574c04b1"},"cell_type":"code","source":"signals = X.reshape((len(X), X.shape[1]*X.shape[2]))\nfeatures = []\n\nfor signal in signals:\n features.append(entropy_and_fractal_dim(signal))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cd9c1e6a1233c57ccbc14a3619ccf1bf76f835df"},"cell_type":"code","source":"features = np.array(features).reshape((len(features), 7))\nscaler = MinMaxScaler(feature_range=(0, 1))\nscaler.fit(features)\nfeatures = scaler.transform(features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"68eeeb6f88422dfb36d55dbab87afdfb2b689829"},"cell_type":"code","source":"print(X.shape, y.shape)\n# save data into file, a numpy specific format\nnp.save(\"X.npy\",X)\nnp.save(\"y.npy\",y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"524f3ddfbf9aaf53483a0c49501e00015e88c969"},"cell_type":"markdown","source":"
Building the RNN Model Architecture
"},{"metadata":{"trusted":true,"_uuid":"7464aa9d8b1bd941235e56b1b2d3ad02fe17878a"},"cell_type":"code","source":"# This is NN LSTM Model creation\ndef model_lstm(input_shape, feat_shape):\n inp = Input(shape=(input_shape[1], input_shape[2],))\n feat = Input(shape=(feat_shape[1],))\n\n bi_lstm_1 = Bidirectional(CuDNNLSTM(128, return_sequences=True), merge_mode='concat')(inp)\n bi_lstm_2 = Bidirectional(CuDNNGRU(64, return_sequences=True), merge_mode='concat')(bi_lstm_1)\n \n attention = Attention(input_shape[1])(bi_lstm_2)\n \n x = concatenate([attention, feat], axis=1)\n x = Dense(64, activation=\"relu\")(x)\n x = Dense(1, activation=\"sigmoid\")(x)\n \n model = Model(inputs=[inp, feat], outputs=x)\n model.compile(loss='binary_crossentropy', optimizer='adam', metrics=[matthews_correlation])\n \n return model","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0b0c205f394c57d65b3d5704d06603508dd87c7e"},"cell_type":"markdown","source":"
Applying Stratified K Fold Technique
"},{"metadata":{"trusted":true,"_uuid":"d056fe026b7578ab106d3aba5aef7b032b7d8cf5"},"cell_type":"code","source":"# Here is where the training happens\n# First, create a set of indexes of the 5 folds\nsplits = list(StratifiedKFold(n_splits=N_SPLITS, shuffle=True, random_state=2019).split(X, y))\npreds_val = []\ny_val = []\n# Then, iteract with each fold\n# If you dont know, enumerate(['a', 'b', 'c']) returns [(0, 'a'), (1, 'b'), (2, 'c')]\nfor idx, (train_idx, val_idx) in enumerate(splits):\n K.clear_session() # I dont know what it do, but I imagine that it \"clear session\" :)\n print(\"Beginning fold {}\".format(idx+1))\n # use the indexes to extract the folds in the train and validation data\n train_X, train_feat, train_y, val_X, val_feat, val_y = X[train_idx], features[train_idx], y[train_idx], X[val_idx], features[val_idx], y[val_idx]\n # instantiate the model for this fold\n model = model_lstm(train_X.shape, features.shape)\n # This checkpoint helps to avoid overfitting. It just save the weights of the model if it delivered an\n # validation matthews_correlation greater than the last one.\n ckpt = ModelCheckpoint('weights_{}.h5'.format(idx), save_best_only=True, save_weights_only=True, verbose=1, monitor='val_matthews_correlation', mode='max')\n # Train, train, train\n model.fit([train_X, train_feat], train_y, batch_size=128, epochs=50, validation_data=([val_X, val_feat], val_y), callbacks=[ckpt])\n # loads the best weights saved by the checkpoint\n model.load_weights('weights_{}.h5'.format(idx))\n # Add the predictions of the validation to the list preds_val\n preds_val.append(model.predict([val_X, val_feat], batch_size=512))\n # and the val true y\n y_val.append(val_y)\n\n# concatenates all and prints the shape \npreds_val = np.concatenate(preds_val)[...,0]\ny_val = np.concatenate(y_val)\npreds_val.shape, y_val.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5d8e87bb034aed8f8b315b8ce8e8da393aa6df80"},"cell_type":"code","source":"import tensorflow as tf\n\ndef threshold_search(y_true, y_proba):\n best_threshold = 0\n best_score = 0\n for threshold in tqdm([i * 0.01 for i in range(100)]):\n score = K.eval(matthews_correlation(y_true.astype(np.float64), (y_proba > threshold).astype(np.float64)))\n if score > best_score:\n best_threshold = threshold\n best_score = score\n search_result = {'threshold': best_threshold, 'matthews_correlation': best_score}\n return search_result","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"92cb38e80553c4c4bc43bb1e39c6bd993b65b275"},"cell_type":"markdown","source":"
Identify the Best Threshold
"},{"metadata":{"_uuid":"c252624065ee8551ae88a2c10b343e4fee3ecc80","trusted":true},"cell_type":"code","source":"best_threshold = threshold_search(y_val, preds_val)['threshold']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8fc5b3dceb6f9a2b3882af4342db975c4f487062","trusted":true},"cell_type":"code","source":"best_threshold","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ae9bd3fa9d8c0781c0708846bb7f2a9f9e6cbd3c","trusted":true},"cell_type":"code","source":"meta_test = pd.read_csv('../input/metadata_test.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3eb186d032f79c99ffba05dd1a7fabb77e13cec5","trusted":true},"cell_type":"code","source":"meta_test = meta_test.set_index(['signal_id'])\nmeta_test.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6f8e94387f625bff0a9a6289e1ee038908bc5856","trusted":true},"cell_type":"code","source":"first_sig = meta_test.index[0]\nn_parts = 10\nmax_line = len(meta_test)\npart_size = int(max_line / n_parts)\nlast_part = max_line % n_parts\nprint(first_sig, n_parts, max_line, part_size, last_part, n_parts * part_size + last_part)\n# Here we create a list of lists with start index and end index for each of the 10 parts and one for the last partial part\nstart_end = [[x, x+part_size] for x in range(first_sig, max_line + first_sig, part_size)]\nstart_end = start_end[:-1] + [[start_end[-1][0], start_end[-1][0] + last_part]]\nprint(start_end)\nX_test = []\n# now, very like we did above with the train data, we convert the test data part by part\n# transforming the 3 phases 800000 measurement in matrix (160,57)\nfor start, end in start_end:\n subset_test = pq.read_pandas('../input/test.parquet', columns=[str(i) for i in range(start, end)]).to_pandas()\n for i in tqdm(subset_test.columns):\n id_measurement, phase = meta_test.loc[int(i)]\n subset_test_col = subset_test[i]\n subset_trans = transform_ts(subset_test_col)\n X_test.append([i, id_measurement, phase, subset_trans])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af9aa6b2b8f8a2beda1a02ff998e3072fcad8d06","trusted":true},"cell_type":"code","source":"X_test_input = np.asarray([np.concatenate([X_test[i][3],X_test[i+1][3], X_test[i+2][3]], axis=1) for i in range(0,len(X_test), 3)])\nnp.save(\"X_test.npy\",X_test_input)\nX_test_input.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eee50e9c6e6e10791563c7ca278c4539e216a162"},"cell_type":"code","source":"signals = X_test_input.reshape((len(X_test_input), X_test_input.shape[1]*X_test_input.shape[2]))\nfeatures_test = []\n\nfor signal in signals:\n features_test.append(entropy_and_fractal_dim(signal))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"979fd9a6a88f514921e694f3ab31de347c6d5462"},"cell_type":"code","source":"features_test = np.array(features_test).reshape((len(features_test), 7))\nfeatures_test = scaler.transform(features_test)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"222d3d99e51f4ef6629a6c47864440721c9b4be2"},"cell_type":"markdown","source":"
Create Submission File
"},{"metadata":{"_uuid":"cfd265d3e07c4cc1679d2c4d55fe7de631c813e7","trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv')\nprint(len(submission))\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2f7342296138f6bfd3e9cedd029e1035de3b98fc","trusted":true},"cell_type":"code","source":"preds_test = []\nfor i in range(N_SPLITS):\n model.load_weights('weights_{}.h5'.format(i))\n pred = model.predict([X_test_input, features_test], batch_size=300, verbose=1)\n pred_3 = []\n for pred_scalar in pred:\n for i in range(3):\n pred_3.append(pred_scalar)\n preds_test.append(pred_3)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4b212616b85a0c704ad8dd11d2ea88c818a4d2a2","trusted":true},"cell_type":"code","source":"preds_test = (np.squeeze(np.mean(preds_test, axis=0)) > best_threshold).astype(np.int)\npreds_test.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"347a76deca88ff3834c2b372391ada6b30a6ae08","trusted":true},"cell_type":"code","source":"submission['target'] = preds_test\nsubmission.to_csv('submission.csv', index=False)\nsubmission.head()","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}