{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Introduction\n### In this kernel, I employ a new type of feature generation technique not used my most kernels in this competition.\n### I will do EDA on new, fresh features and analyze the visualizations with detailed comments and descriptions.\n### I will do training and inference with a stacked BiLSTM model (on these features) in a future version (you can try it out yourself)\n### Thanks to Bruno Aquino, whose work this method is based on. Check out [this kernel](https://www.kaggle.com/braquino/5-fold-lstm-attention-fully-commented-0-694) of his. The entropy and fractal dimension functions are from Raphael Vallat's [entropy](https://github.com/raphaelvallat/entropy/blob/master/entropy) repository\n### PLEASE UPVOTE IF YOU LIKE THIS KERNEL\n"},{"metadata":{},"cell_type":"markdown","source":"<center><img src=\"https://i.imgur.com/hBPv3fh.png\" width=\"750px\"></center>"},{"metadata":{},"cell_type":"markdown","source":"### Import necessary libraries"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":false},"cell_type":"code","source":"import os\nimport gc\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\ntqdm.pandas()\n\nfrom numba import jit\nfrom math import log, floor\nfrom sklearn.neighbors import KDTree\nfrom scipy.signal import periodogram, welch\n\nfrom keras.layers import *\nfrom keras.models import *\nfrom tqdm import tqdm\nfrom sklearn.model_selection import train_test_split \nfrom keras import backend as K\nfrom keras import optimizers\nfrom sklearn.model_selection import GridSearchCV, KFold\nfrom keras.callbacks import *\nfrom keras import activations\nfrom keras import regularizers\nfrom keras import initializers\nfrom keras import constraints\nfrom keras.engine import Layer\nfrom keras.engine import InputSpec\nfrom keras.objectives import categorical_crossentropy\nfrom keras.objectives import sparse_categorical_crossentropy\nfrom keras.utils import plot_model\nfrom keras.utils.vis_utils import model_to_dot\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom IPython.display import SVG\n\nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Initialize necessay constants"},{"metadata":{"trusted":false},"cell_type":"code","source":"SIGNAL_LEN = 150000\nMIN_NUM = -27\nMAX_NUM = 28","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Download seismic signal data along with targets (time left for occurance of laboratory earthquake)"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"seismic_signals = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Extract the acoustic data and targets from the dataframe\nNote : I delete the original dataframe to save memory"},{"metadata":{"trusted":false},"cell_type":"code","source":"acoustic_data = seismic_signals.acoustic_data\ntime_to_failure = seismic_signals.time_to_failure\ndata_len = len(seismic_signals)\ndel seismic_signals\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Break the data down into parts\nWe have one long array of seismic data. We will break it down into chunks of size 150k (**SIGNAL_LEN**) and each chunk will be one signal in our data (this is because each segment in the test data has length 150k). The **time_to_failure** at the last time step of each segment becomes the target associated with that segment."},{"metadata":{"trusted":false},"cell_type":"code","source":"signals = []\ntargets = []\n\nfor i in range(data_len//SIGNAL_LEN):\n    min_lim = SIGNAL_LEN * i\n    max_lim = min([SIGNAL_LEN * (i + 1), data_len])\n    \n    signals.append(list(acoustic_data[min_lim : max_lim]))\n    targets.append(time_to_failure[max_lim])\n    \ndel acoustic_data\ndel time_to_failure\ngc.collect()\n    \nsignals = np.array(signals)\ntargets = np.array(targets)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Functions for preparing signal features"},{"metadata":{},"cell_type":"markdown","source":"### Scaling the signals\nThis function scales the seismic signals from its original range (-27 to 28 : this where 99% of the data lies) to the range (-1, 1)"},{"metadata":{"trusted":false},"cell_type":"code","source":"def min_max_transfer(ts, min_value, max_value, range_needed=(-1,1)):\n    ts_std = (ts - min_value) / (max_value - min_value)\n\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":{},"cell_type":"markdown","source":"### Extracting features from each part of the segment\nThe original long seismic signal has already been broken down into several segments. The segments are scaled using the **min_max_transfer** function. Then, we break down each segment into several parts. Usual features such as mean, standard deviation, range, percentiles etc are calculated over each part of the segment and now, each part of the segment is represented by its own list of such features. Finally, the representations of all the small parts of the segment are strung together into a time series. This time series becomes the representation of that segment."},{"metadata":{"trusted":false},"cell_type":"code","source":"def transform_ts(ts, n_dim=160, min_max=(-1,1)):\n    ts_std = min_max_transfer(ts, min_value=MIN_NUM, max_value=MAX_NUM)\n    bucket_size = int(SIGNAL_LEN / n_dim)\n    new_ts = []\n    for i in range(0, SIGNAL_LEN, bucket_size):\n        ts_range = ts_std[i:i + bucket_size]\n        mean = ts_range.mean()\n        std = ts_range.std()\n        std_top = mean + std\n        std_bot = mean - std\n        percentil_calc = ts_range.quantile([0, 0.01, 0.25, 0.50, 0.75, 0.99, 1])\n        max_range = ts_range.quantile(1) - ts_range.quantile(0)\n        relative_percentile = percentil_calc - mean\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":{},"cell_type":"markdown","source":"### Prepare the final signal features\nThe time series representations of all segments of the signal are calulated and concatenated. This results in a 3D tensor containing the time series representations of all segments in the signal."},{"metadata":{"trusted":false},"cell_type":"code","source":"def prepare_data(start, end):\n    train = pd.DataFrame(np.transpose(signals[int(start):int(end)]))\n    X = []\n    for id_measurement in tqdm(train.index[int(start):int(end)]):\n        X_signal = transform_ts(train[id_measurement])\n        X.append(X_signal)\n    X = np.asarray(X)\n    return X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Implement the feature generation process"},{"metadata":{"trusted":false},"cell_type":"code","source":"X = []\n\ndef load_all():\n    total_size = len(signals)\n    for start, end in [(0, int(total_size))]:\n        X_temp = prepare_data(start, end)\n        X.append(X_temp)\n        \nload_all()\nX = np.concatenate(X)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here is the shape of X. There are a total of 4194 segments. Each segment is divided into 161 parts and each part is represented by a list of 19 features (mean, stddev etc). Therefore, X is a 3D tensor with shape (4194, 161, 19)."},{"metadata":{"trusted":false},"cell_type":"code","source":"X.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Flattening the features and doing basic EDA with seaborn"},{"metadata":{},"cell_type":"markdown","source":"Now, we flatten the 2D tensors associated with each segment into 1D arrays. Now, each data point (segment) is represented by a 1D array.\n\nHere are some flattened 1D arrays (with sparse selection) visualized with **matplotlib**."},{"metadata":{"trusted":false},"cell_type":"code","source":"shape = X.shape\nnew_signals = X.reshape((shape[0], shape[1]*shape[2]))\n\nsparse_signals = []\nfor i in range(3):\n    sparse_signal = []\n    for j in range(len(new_signals[i])):\n        if j % 3 == 0:\n            sparse_signal.append(new_signals[i][j])\n    sparse_signals.append(sparse_signal)\n\nplt.plot(sparse_signals[0], 'purple')\nplt.show()\nplt.plot(sparse_signals[1], 'mediumvioletred')\nplt.show()\nplt.plot(sparse_signals[2], 'crimson')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Permutation Entropy\nThe permutation entropy is a complexity measure for time-series first introduced by Bandt and Pompe in 2002. It represents the information contained\nin comparing *n* consecutive values of the time series. It is a measure of entropy or disorderliness in a time series."},{"metadata":{"trusted":false},"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","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"perm_entropies = np.array([perm_entropy(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=perm_entropies, y=targets, kind='kde', color='orangered')\nplot.set_axis_labels('perm_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) along a line with positive slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=perm_entropies, y=targets, kind='hex', color='orangered')\nplot.set_axis_labels('perm_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a positively-sloped line."},{"metadata":{},"cell_type":"markdown","source":"#### Scatterplot with line of best fit"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=perm_entropies, y=targets, kind='reg', color='orangered')\nplot.set_axis_labels('perm_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear positive slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **positive correlation** between the permutation entropy of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### Approximate Entropy\nApproximate entropy is a technique used to quantify the amount of\nregularity and the unpredictability of fluctuations over time-series data.\nSmaller values indicates that the data is more regular and predictable."},{"metadata":{"trusted":false},"cell_type":"code","source":"def _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])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"app_entropies = np.array([app_entropy(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=app_entropies, y=targets, kind='kde', color='magenta')\nplot.set_axis_labels('app_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) around a line with negative slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=app_entropies, y=targets, kind='hex', color='magenta')\nplot.set_axis_labels('app_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a negatively-sloped line."},{"metadata":{},"cell_type":"markdown","source":"#### Scatterplot with line of best fit"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=app_entropies, y=targets, kind='reg', color='magenta')\nplot.set_axis_labels('app_entropy', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear negative slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **negative correlation** between the approximate entropy of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### Higuchi Fractal Dimension\nThe Higuchi fractal dimension is a method to calculate the [fractal dimension](https://en.wikipedia.org/wiki/Fractal_dimension) of any two-dimensional curve. Generally, curves with higher fractal dimension are \"rougher\" or more \"complex\" (higher entropy)."},{"metadata":{"trusted":false},"cell_type":"code","source":"@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)\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"higuchi_fds = np.array([higuchi_fd(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate KDE distribution plot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=higuchi_fds, y=targets, kind='kde', color='crimson')\nplot.set_axis_labels('higuchi_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) around a line with negative slope."},{"metadata":{},"cell_type":"markdown","source":"#### Bivariate hexplot"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=higuchi_fds, y=targets, kind='hex', color='crimson')\nplot.set_axis_labels('higuchi_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot is also darkest around a negatively-sloped line."},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=higuchi_fds, y=targets, kind='reg', color='crimson')\nplot.set_axis_labels('higuchi_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit in the scatterplot has a clear negative slope."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see a somewhat **negative correlation** between the Higuchi fractal dimension of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"### Katz Fractal Dimension\nThe Katz fractal dimension is yet another way to calculate the fractal dimension of a two-dimensional curve."},{"metadata":{"trusted":false},"cell_type":"code","source":"def 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))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"katz_fds = np.array([katz_fd(new_signal) for new_signal in new_signals])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=katz_fds, y=targets, kind='kde', color='forestgreen')\nplot.set_axis_labels('katz_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The KDE plot has highest density (darkness) around an almost vertical line."},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=katz_fds, y=targets, kind='hex', color='forestgreen')\nplot.set_axis_labels('katz_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The hexplot also has highest density around an almost vertical line."},{"metadata":{"trusted":false},"cell_type":"code","source":"plot = sns.jointplot(x=katz_fds, y=targets, kind='reg', color='forestgreen')\nplot.set_axis_labels('katz_fd', 'time_to_failure', fontsize=16)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The line of best fit has a very slight positive correlation."},{"metadata":{},"cell_type":"markdown","source":"From the above three plots we can see an **very slight positive correlation (or maybe no correlation)** between the Katz fractal dimension of the flattened feature array and the time left for the laboratory earthquake to occur."},{"metadata":{},"cell_type":"markdown","source":"## Example model to use for these features"},{"metadata":{},"cell_type":"markdown","source":"### Attention layer\nAttention mechanisms in neural networks serve to orient perception as well as memory access (you might even say perception is just a very short-term subset of all memory). Attention filters the perceptions that can be stored in memory, and filters them again on a second pass when they are to be retrieved from memory. Attention can be aimed at the present and the past."},{"metadata":{"trusted":false},"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_type":"markdown","source":"### Example neural network architecture (stacked BiLSTM with Attention)\n\n#### Create a model that :\n\n1). Passes the 3D input through a Bidirectional LSTM layer with 128 units\n\n2). Passes the output of the previous layer through another Bidirectional LSTM layer with 64 units\n\n3). Then passes the result through an Attention layer\n\n4). The Attention layer's is passed through a Dense layer with 64 units\n\n5). Finally, the output is obtained by passing the previous layer's output through a Dense layer with one neuron"},{"metadata":{"trusted":false},"cell_type":"code","source":"def model_lstm(input_shape):\n    inp = Input(shape=(input_shape[1], input_shape[2],))\n\n    bi_lstm = Bidirectional(LSTM(128, return_sequences=True))(inp)\n    bi_gru = Bidirectional(LSTM(64, return_sequences=True))(bi_lstm)\n    \n    attention = Attention(input_shape[1])(bi_gru)\n    \n    x = Dense(64, activation=\"relu\")(attention)\n    x = Dense(1)(x)\n    \n    model = Model(inputs=inp, outputs=x)\n    model.compile(loss='mae', optimizer='adam')\n    \n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"sample_model = model_lstm(X.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Model summary"},{"metadata":{"trusted":false},"cell_type":"code","source":"sample_model.summary()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Model plot visualization "},{"metadata":{"trusted":false},"cell_type":"code","source":"SVG(model_to_dot(sample_model).create(prog='dot', format='svg'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### I will do training and inference with this model in a future version (you can try it out yourself)"},{"metadata":{},"cell_type":"markdown","source":"### That's it ! Thanks for reading this kernel :)"},{"metadata":{},"cell_type":"markdown","source":"### Please post your feedback and suggestions in the comments."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.4"}},"nbformat":4,"nbformat_minor":1}