{"cells":[{"metadata":{"_uuid":"e154a47bf09b8770980486e87786317a1b3038e1"},"cell_type":"markdown","source":"### This kernel is forked from [this wonderful kernel](https://www.kaggle.com/braquino/5-fold-lstm-attention-fully-commented-0-694) by Bruno Aquino. I have just stacked the recurrent layers together. Additionally, I have done some feature engineering based on time series entropy and fractal dimension. This achieves an LB score of 0.699.\n\nFeature Engineering : [Raphael Vallat's entropy repository](https://github.com/raphaelvallat/entropy/blob/master/entropy/utils.py)\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","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\nimport gc\nimport pywt\nimport scipy\nfrom keras.layers import * # Keras is the most friendly Neural Network library, this Kernel use a lot of layers classes\nfrom keras.models import *\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 # 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 keras import activations\nfrom keras import regularizers\nfrom keras import initializers\nfrom keras import constraints\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\nfrom keras.engine import Layer\nfrom keras.engine import InputSpec\nfrom keras.objectives import categorical_crossentropy\nfrom keras.objectives import sparse_categorical_crossentropy\nfrom scipy import stats\nfrom statsmodels.robust import mad\nfrom scipy import signal\nfrom scipy.signal import butter\ngc.enable()\n# import gplearn\n# from gplearn.genetic import SymbolicTransformer\n# from gplearn.functions import make_function","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f2d087dcc22505502c97227497bde9f4a5a92fb2"},"cell_type":"code","source":"import warnings\n\n# Suppress pandas future warnings, I am using different library versions locally\n# that do not raise warnings.\nwarnings.simplefilter(action='ignore', category=FutureWarning)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6e6379386e44afc69bee8895a52da22199e888fb"},"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":{"trusted":true,"_uuid":"c3340ee96becb5ca8f075d9c44b7df383ddba5ee"},"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_true = K.cast(y_true, float)\n    y_pred = K.cast(y_pred, float)\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":{"trusted":true,"_uuid":"eda7ea366117d1ce8e5fce69e5bba333821d8b48"},"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":{"_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":{"trusted":true,"_uuid":"3315b491ff127ac7a351c48a01ddd6ef8f200791"},"cell_type":"code","source":"len(df_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"26df6c7fbfecd537404866faec13d1238ae3ebc6"},"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":{"trusted":true,"_uuid":"7b0717b14bcfcba1f48d33c8161ae51c778687af"},"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_transf(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":{"trusted":true,"_uuid":"c6137bbbe75c3a1509a5f98e08805dbbd492aa37"},"cell_type":"code","source":"# This is one of the most important peace of code of this Kernel\n# Any power line contain 3 phases of 800000 measurements, or 2.4 millions data \n# It would be praticaly impossible to build a NN with an input of that size\n# The ideia here is to reduce it each phase to a matrix of <n_dim> bins by n features\n# Each bean is a set of 5000 measurements (800000 / 160), so the features are extracted from this 5000 chunk data.\ndef transform_ts(ts, n_dim=160, min_max=(-1,1)):\n    #Apply HPF\n    #x_hp = high_pass_filter(ts, low_cutoff=10000, sample_rate=sample_rate)\n    # Apply denoising\n    #x_dn = denoise_signal(ts, wavelet='haar', level=1)\n    #qqq = pd.Series(data = x_dn)\n    #ts = qqq\n    #del x_dn,qqq\n    #gc.collect()\n    #print('Cycle done')\n    \n    # convert data into -1 to 1\n    ts_std = min_max_transf(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\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        # now, we just add all the features to new_ts and convert it to np.array\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":{"trusted":true,"_uuid":"cf0839d7a351c66b42b6ed01c63cc4c8ffc95f6c"},"cell_type":"code","source":"# 800,000 data points taken over 20 ms\n# Grid operates at 50hz, 0.02 * 50 = 1, so 800k samples in 20 milliseconds will capture one complete cycle\nn_samples = 800000\n\n# Sample duration is 20 miliseconds\nsample_duration = 0.02\n\n# Sample rate is the number of samples in one second\n# Sample rate will be 40mhz\nsample_rate = n_samples * (1 / sample_duration)\n\ndef maddest(d, axis=None):\n    \"\"\"\n    Mean Absolute Deviation\n    \"\"\"\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef high_pass_filter(x, low_cutoff=1000, sample_rate=sample_rate):\n    \"\"\"\n    From @randxie https://github.com/randxie/Kaggle-VSB-Baseline/blob/master/src/utils/util_signal.py\n    Modified to work with scipy version 1.1.0 which does not have the fs parameter\n    \"\"\"\n    \n    # nyquist frequency is half the sample rate https://en.wikipedia.org/wiki/Nyquist_frequency\n    nyquist = 0.5 * sample_rate\n    norm_low_cutoff = low_cutoff / nyquist\n    \n    # Fault pattern usually exists in high frequency band. According to literature, the pattern is visible above 10^4 Hz.\n    # scipy version 1.2.0\n    #sos = butter(10, low_freq, btype='hp', fs=sample_fs, output='sos')\n    \n    # scipy version 1.1.0\n    sos = butter(10, Wn=[norm_low_cutoff], btype='highpass', output='sos')\n    filtered_sig = scipy.signal.sosfilt(sos, x)\n\n    return filtered_sig\n\ndef denoise_signal( x, wavelet='db4', level=1):\n    \"\"\"\n    1. Adapted from waveletSmooth function found here:\n    http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/\n    2. Threshold equation and using hard mode in threshold as mentioned\n    in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:\n    http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    \"\"\"\n    \n    # Decompose to get the wavelet coefficients\n    coeff = pywt.wavedec( x, wavelet, mode=\"per\" )\n    \n    # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf\n    # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation\n    sigma = (1/0.6745) * maddest( coeff[-level] )\n\n    # Calculte the univeral threshold\n    uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) )\n    coeff[1:] = ( pywt.threshold( i, value=uthresh, mode='hard' ) for i in coeff[1:] )\n    \n    # Reconstruct the signal using the thresholded coefficients\n    return pywt.waverec( coeff, wavelet, mode='per' )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7460e718a605803f1d9e4fbec61750a0deb02a47"},"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 prep_data(start, end):\n    # load a piece of data from file\n    praq_train = pq.read_pandas('../input/train.parquet', columns=[str(i) for i in range(start, end)]).to_pandas()\n\n    X = []\n    y = []\n    # using tdqm to evaluate processing time\n    # takes each index from df_train and iteract it from start to end\n    # it is divided by 3 because for each id_measurement there are 3 id_signal, and the start/end parameters are id_signal\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\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":{"trusted":true,"_uuid":"58a9c267dfd21bdcdeeed4efae2e7628b2c02b09"},"cell_type":"code","source":"def peaks(x):\n    y = denoise_signal(x, wavelet='haar', level=1)\n    peaks, properties = scipy.signal.find_peaks(y)\n    widths = scipy.signal.peak_widths(y, peaks)[0]\n    prominences = scipy.signal.peak_prominences(y, peaks)[0]\n    \n    count= peaks.size,\n    width_mean= widths.mean() if widths.size else -1.,\n    width_max= widths.max() if widths.size else -1.,\n    width_min= widths.min() if widths.size else -1.,\n    prominence_mean= prominences.mean() if prominences.size else -1.,\n    prominence_max= prominences.max() if prominences.size else -1.,\n    prominence_min= prominences.min() if prominences.size else -1.,\n    return count[0], width_mean[0], width_max[0], width_min[0], prominence_mean[0], prominence_max[0], prominence_min[0]\n\n\ndef denoised_std(x):\n    return np.std(denoise_signal(x, wavelet='haar', level=1))\n\n\ndef signal_entropy(x):\n\n    y = denoise_signal(x, wavelet='haar', level=1)\n\n    for i in range(3):\n        max_pos = y.argmax()\n        y[max_pos - 1000:max_pos + 1000] = 0.\n\n    return stats.entropy(np.histogram(y, 15)[0])\n\n\ndef detail_coeffs_entropy(x, wavelet='db1'):\n\n    c_a, c_d = pywt.dwt(x, wavelet)\n\n    return stats.entropy(np.histogram(c_d, 15)[0])\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"52dc826ab9ee1dd56c9fb29bd5c1b2d26b5928bf","scrolled":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 = prep_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":"723f1b9f68537d8c7964a2087627bf9f02a05193"},"cell_type":"code","source":"X.shape #1462 160 57 что-то такое","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dcacaf0e7a9553189d5385a446e8ca1892416fff"},"cell_type":"code","source":"print(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ebff0c50964f94c3101f85502a5895e006fd98da"},"cell_type":"code","source":"print(y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ad7c129da070baeb2245e39f921ea2ed74346f2f"},"cell_type":"code","source":"#    train_length = end\n#    print(praq_train.head())\n#    for i in range(train_length):\n#        signal_id = str(i)\n#        meta_row = df_train[df_train['signal_id'] == i].reset_index()\n#        measurement = str(meta_row['id_measurement'].values[0])\n#        signal_id = str(meta_row['signal_id'].values[0])\n#        phase = str(meta_row['phase'].values[0])\n#        subset_train_row = praq_train[signal_id]\n#        # Apply high pass filter with low cutoff of 10kHz, this will remove the low frequency 50Hz sinusoidal motion in the signal\n#        x_hp = high_pass_filter(subset_train_row, low_cutoff=10000, sample_rate=sample_rate)\n#        print('HP done')\n#        # Apply denoising\n#        x_dn = denoise_signal(x_hp, wavelet='haar', level=1)\n#        print('DN done')\n#        qqq = pd.Series(data = x_dn)\n#        praq_train[str(i)] = qqq\n#        del signal_id, meta_row, measurement, phase, subset_train_row, x_hp, x_dn,qqq\n#        gc.collect()\n#        print('Cycle done')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ee3f68b87976db3d518dee0e81cebde379c2205f"},"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":"51f450c5943ede860366a22070188d60814c5fe4"},"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":"5c5d5906b02955513e261e8e7f581d29fcacd9db"},"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),scipy.stats.kurtosis(x),denoised_std(x),signal_entropy(x),detail_coeffs_entropy(x, wavelet='db1')]+list(peaks(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c79f2d33c77bc4bc8e1d8c040250787995c61098"},"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":"ff829cda05a7ad9583c5e6d1d360d89e01c780d1"},"cell_type":"code","source":"features = np.array(features).reshape((len(features), 18))\nscaler = MinMaxScaler(feature_range=(0, 1))\nscaler.fit(features)\nfeatures = scaler.transform(features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2aab27f533058d68d34b16d9993238d052b4608c"},"cell_type":"code","source":"# model = Sequential([Dense(1, activation='sigmoid', input_shape=(7,))])\n# model.compile(loss='binary_crossentropy', optimizer='adam')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b147c8c0fdb69373dbcb25640a192242a7ab1598"},"cell_type":"code","source":"# model.fit(x, np.array(list(np.zeros((len(features_zero)))) + list(np.ones((len(features_one))))), epochs=1000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cc7eb66ecfeca0ada7b5a7ce680fd20e10b98421"},"cell_type":"code","source":"# Weights and bias for logistic regression of entropy features\n# weights = model.get_weights()[0]\n# biases = model.get_weights()[1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3b980b5a62c0db389f82bbe10bec7ed3731c9c6c"},"cell_type":"code","source":"# print(weights)\n# print(biases)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"69171fae6e8eec531735bbb1224cb8f4152b79a5"},"cell_type":"code","source":"# e = 2.718281\n\n# exp = make_function(function=lambda x: e**x, name='exp', arity=1)\n# tanh = make_function(function=lambda x: (e**(2*x) - 1)/(e**(2*x) + 1), name='tanh', arity=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d6c88e6595bb57b859f7ffe035ae1c8ce041f592"},"cell_type":"code","source":"# function_set = ['add', 'sub', 'mul', 'div',\n                # 'sqrt', 'log', 'abs', 'neg', 'inv',\n                # 'max', 'min', exp, tanh]\n\n# gp = SymbolicTransformer(generations=10, population_size=2000,\n                         # hall_of_fame=100, n_components=50,\n                         # function_set=function_set,\n                         # parsimony_coefficient=0.0005,\n                         # max_samples=0.9, verbose=1,\n                         # random_state=0, n_jobs=3)\n\n# gp.fit(features, np.array(list(np.zeros((len(features_zero)))) + list(np.ones((len(features_one))))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5ba22c728ece64aff2c93164370b02df33c516bc"},"cell_type":"code","source":"# engineered_features = gp._programs\n\n# for i in range(len(engineered_features)):\n    # for engineered_feature in engineered_features[i]:\n        # if engineered_feature != None:\n            # print(engineered_feature)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"51ad0e25b00536de6170168499923d82ae1d735f"},"cell_type":"code","source":"# The X shape here is very important. It is also important undertand a little how a LSTM works\n# X.shape[0] is the number of id_measuremts contained in train data\n# X.shape[1] is the number of chunks resultant of the transformation, each of this date enters in the LSTM serialized\n# This way the LSTM can understand the position of a data relative with other and activate a signal that needs\n# a serie of inputs in a specifc order.\n# X.shape[3] is the number of features multiplied by the number of phases (3)\nprint(X.shape, y.shape)\n# save data into file, a numpy specific format\nnp.save(\"X.npy\",X)\nnp.save(\"y.npy\",y)\nnp.save(\"features.npy\",features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"289bc7d1ab8048a60025801b457f8df1d848acbc"},"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 = Bidirectional(CuDNNLSTM(128, return_sequences=True), merge_mode='concat')(inp)\n    x = Dropout(0.23)(bi_lstm)\n    bi_gru = Bidirectional(CuDNNGRU(64, return_sequences=True), merge_mode='concat')(x)\n    \n    attention = Attention(input_shape[1])(bi_gru)\n    \n    x = concatenate([attention, feat], axis=1)\n    x = Dropout(0.2)(x)\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":{"trusted":true,"_uuid":"8d6f4ca319c383b1b4f671a37c5a324136e7a466"},"cell_type":"code","source":"# Here is where the training happens\n\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=70, 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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d28151fd0be9fd9762f3f55e307d82f89bfbd291"},"cell_type":"code","source":"# The output of this kernel must be binary (0 or 1), but the output of the NN Model is float (0 to 1).\n# So, find the best threshold to convert float to binary is crucial to the result\n# this piece of code is a function that evaluates all the possible thresholds from 0 to 1 by 0.01\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":{"trusted":true,"_uuid":"6fee7f722ed08bc1453a822a4371ed2d48e08abc"},"cell_type":"code","source":"optimal_values = threshold_search(y_val, preds_val)\nbest_threshold = optimal_values['threshold']\nbest_score = optimal_values['matthews_correlation']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4b96eed474f3712dfedae9cdd8fa6983cc6efc80"},"cell_type":"code","source":"print(\"Optimal Threshold : \" + str(best_threshold))\nprint(\"Best Matthews Correlation : \" + str(best_score))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ae9bd3fa9d8c0781c0708846bb7f2a9f9e6cbd3c"},"cell_type":"code","source":"%%time\n# Now load the test data\n# This first part is the meta data, not the main data, the measurements\nmeta_test = pd.read_csv('../input/metadata_test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3eb186d032f79c99ffba05dd1a7fabb77e13cec5"},"cell_type":"code","source":"meta_test = meta_test.set_index(['signal_id'])\nmeta_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6f8e94387f625bff0a9a6289e1ee038908bc5856"},"cell_type":"code","source":"%%time\n# First we daclarete a series of parameters to initiate the loading of the main data\n# it is too large, it is impossible to load in one time, so we are doing it in dividing in 10 parts\nfirst_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[str(i)]\n        subset_trans = transform_ts(subset_test_col)\n        X_test.append([i, id_measurement, phase, subset_trans])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"af9aa6b2b8f8a2beda1a02ff998e3072fcad8d06"},"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":"4f71465cd6d130e27c035e7596b3e3af397f7ed8"},"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":"1062376c23cab336a465775050e96db725f49f9f"},"cell_type":"code","source":"features_test = np.array(features_test).reshape((len(features_test), 18))\nfeatures_test = scaler.transform(features_test)\nnp.save(\"features_test.npy\",features_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cfd265d3e07c4cc1679d2c4d55fe7de631c813e7"},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv')\nprint(len(submission))\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2f7342296138f6bfd3e9cedd029e1035de3b98fc"},"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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9f76c471eaf983707d446c5081ab3d50c4e40ea5"},"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":{"trusted":true,"_uuid":"b35723f85d494b4b6ec630dd7c79135a110a4062"},"cell_type":"code","source":"submission['target'] = preds_test\nsubmission.to_csv('submission.csv', index=False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"0a6c2d2382d824d5287b08e7ba9932c8b5310259"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}