{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook contains the winning solution to the VSB Power Line Fault Detection competition.\n\nKey Information:\n\n    * Simple LightGBM model\n    * Standard 5-fold CV\n    * Uses just 9 features\n    * Features all calculated using the peaks in the signals\n    * Run times:\n        * Generating features: ~15 minutes\n        * Training: ~2 minutes\n        * Prediction: ~10 seconds"},{"metadata":{"trusted":true,"_uuid":"36d68d77d86ea10bc97451305d3f1b67c3f17903"},"cell_type":"code","source":"# Set RECALCULATE_PEAK_FEATURES to True to recalculate the per peak features\nRECALCULATE_PEAK_FEATURES = True\n# The PREPROCESSED_DIR is used to load the previous commit's per peak features if \n# RECALCULATE_PEAK_FEATURES is set to False. Note, if the kernel is forked this \n# variable needs to be updated to point to the forked kernel's directory\nPREPROCESSED_DIR = '../input/vsb-1st-place-solution'","execution_count":1,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8dc7d6b7cb5f924b6172d12e0cea21836dceb6ee"},"cell_type":"code","source":"USE_SIMPLIFIED_VERSION = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"79edf41624070a6cf531e8122838bc12fbf3922e"},"cell_type":"code","source":"if RECALCULATE_PEAK_FEATURES:\n    # If performing preprocessing to calculate the per peak features need a newer version of\n    # numba than the latest kaggle docker image\n    !conda update -y numba","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import copy\nimport gc\nimport os\nimport sys\nimport warnings\n\nfrom IPython.core.display import display, HTML\nimport lightgbm as lgb\nfrom matplotlib import pyplot as plt\nimport numba\nimport numpy as np\nimport pandas as pd\nimport pyarrow\nimport pyarrow.parquet as pq\nimport scipy.stats\nimport seaborn as sns \nimport sklearn\nfrom sklearn.metrics import matthews_corrcoef, precision_recall_fscore_support","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3dc8d667d7684041e3fb53d5fa0b743e64d9eb09"},"cell_type":"code","source":"data_dir = '../input/vsb-power-line-fault-detection'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"292caf2daad402c77f222eadd5dd8e3bf0596433"},"cell_type":"code","source":"for p in [scipy, np, numba, sklearn, lgb, pyarrow, pd]:\n    print('{}=={}'.format(p.__name__, p.__version__))\n    \nprint()\n\nprint(sys.version)\n\n# Tested using versions:\n#scipy==1.1.0\n#numpy==1.16.2\n#numba==0.43.0\n#sklearn==0.20.3\n#lightgbm==2.2.3\n#pyarrow==0.10.0\n#pandas==0.23.4","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"meta_train_df = pd.read_csv(data_dir + '/metadata_train.csv')\nmeta_test_df = pd.read_csv(data_dir + '/metadata_test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"df9d1d8f523ee1eefd8922cee233544f4c61a617"},"cell_type":"code","source":"meta_train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"22c0cb7116b8d29e7bdcfc89132e08cf2c05c645"},"cell_type":"markdown","source":"# Functions"},{"metadata":{"trusted":true,"_uuid":"d26cbe1cf03e563e0a8aec2475681d8a5c9ee150"},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef flatiron(x, alpha=100., beta=1):\n    \"\"\"\n    Flatten signal\n    \n    Creator: Michael Kazachok\n    Source: https://www.kaggle.com/miklgr500/flatiron\n    \"\"\"\n    new_x = np.zeros_like(x)\n    zero = x[0]\n    for i in range(1, len(x)):\n        zero = zero*(alpha-beta)/alpha + beta*x[i]/alpha\n        new_x[i] = x[i] - zero\n    return new_x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1e4161bee610416e9fe21a99186b101a3de5b556"},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef drop_missing(intersect,sample):\n    \"\"\"\n    Find intersection of sorted numpy arrays\n    \n    Since intersect1d sort arrays each time, it's effectively inefficient.\n    Here you have to sweep intersection and each sample together to build\n    the new intersection, which can be done in linear time, maintaining order. \n\n    Source: https://stackoverflow.com/questions/46572308/intersection-of-sorted-numpy-arrays\n    Creator: B. M.\n    \"\"\"\n    i=j=k=0\n    new_intersect=np.empty_like(intersect)\n    while i< intersect.size and j < sample.size:\n        if intersect[i]==sample[j]: # the 99% case\n            new_intersect[k]=intersect[i]\n            k+=1\n            i+=1\n            j+=1\n        elif intersect[i]<sample[j]:\n            i+=1\n        else : \n            j+=1\n    return new_intersect[:k]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"20a88e1924050b1b72b0a723649386c1bcda8123"},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef _local_maxima_1d_window_single_pass(x, w):\n    \n    midpoints = np.empty(x.shape[0] // 2, dtype=np.intp)\n    left_edges = np.empty(x.shape[0] // 2, dtype=np.intp)\n    right_edges = np.empty(x.shape[0] // 2, dtype=np.intp)\n    m = 0  # Pointer to the end of valid area in allocated arrays\n\n    i = 1  # Pointer to current sample, first one can't be maxima\n    i_max = x.shape[0] - 1  # Last sample can't be maxima\n    while i < i_max:\n        # Test if previous sample is smaller\n        if x[i - 1] < x[i]:\n            i_ahead = i + 1  # Index to look ahead of current sample\n\n            # Find next sample that is unequal to x[i]\n            while i_ahead < i_max and x[i_ahead] == x[i]:\n                i_ahead += 1\n                    \n            i_right = i_ahead - 1\n            \n            f = False\n            i_window_end = i_right + w\n            while i_ahead < i_max and i_ahead < i_window_end:\n                if x[i_ahead] > x[i]:\n                    f = True\n                    break\n                i_ahead += 1\n                \n            # Maxima is found if next unequal sample is smaller than x[i]\n            if x[i_ahead] < x[i]:\n                left_edges[m] = i\n                right_edges[m] = i_right\n                midpoints[m] = (left_edges[m] + right_edges[m]) // 2\n                m += 1\n                \n            # Skip samples that can't be maximum\n            i = i_ahead - 1\n        i += 1\n\n    # Keep only valid part of array memory.\n    midpoints = midpoints[:m]\n    left_edges = left_edges[:m]\n    right_edges = right_edges[:m]\n    \n    return midpoints, left_edges, right_edges\n\n@numba.jit(nopython=True)\ndef local_maxima_1d_window(x, w=1):\n    \"\"\"\n    Find local maxima in a 1D array.\n    This function finds all local maxima in a 1D array and returns the indices\n    for their midpoints (rounded down for even plateau sizes).\n    It is a modified version of scipy.signal._peak_finding_utils._local_maxima_1d\n    to include the use of a window to define how many points on each side to use in\n    the test for a point being a local maxima.\n    Parameters\n    ----------\n    x : ndarray\n        The array to search for local maxima.\n    w : np.int\n        How many points on each side to use for the comparison to be True\n    Returns\n    -------\n    midpoints : ndarray\n        Indices of midpoints of local maxima in `x`.\n    Notes\n    -----\n    - Compared to `argrelmax` this function is significantly faster and can\n      detect maxima that are more than one sample wide. However this comes at\n      the cost of being only applicable to 1D arrays.\n    \"\"\"    \n        \n    fm, fl, fr = _local_maxima_1d_window_single_pass(x, w)\n    bm, bl, br = _local_maxima_1d_window_single_pass(x[::-1], w)\n    bm = np.abs(bm - x.shape[0] + 1)[::-1]\n    bl = np.abs(bl - x.shape[0] + 1)[::-1]\n    br = np.abs(br - x.shape[0] + 1)[::-1]\n\n    m = drop_missing(fm, bm)\n\n    return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a9f61289bb745e194988bc272971b42b4d4e1b2e"},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef plateau_detection(grad, threshold, plateau_length=5):\n    \"\"\"Detect the point when the gradient has reach a plateau\"\"\"\n    \n    count = 0\n    loc = 0\n    for i in range(grad.shape[0]):\n        if grad[i] > threshold:\n            count += 1\n        \n        if count == plateau_length:\n            loc = i - plateau_length\n            break\n            \n    return loc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"70a8b049dc199743ae325e418669dba3b55f01fd"},"cell_type":"code","source":"#@numba.jit(nopython=True)\ndef get_peaks(\n    x, \n    window=25,\n    visualise=False,\n    visualise_color=None,\n):\n    \"\"\"\n    Find the peaks in a signal trace.\n    Parameters\n    ----------\n    x : ndarray\n        The array to search.\n    window : np.int\n        How many points on each side to use for the local maxima test\n    Returns\n    -------\n    peaks_x : ndarray\n        Indices of midpoints of peaks in `x`.\n    peaks_y : ndarray\n        Absolute heights of peaks in `x`.\n    x_hp : ndarray\n        An absolute flattened version of `x`.\n    \"\"\"\n    \n    x_hp = flatiron(x, 100, 1)\n\n    x_dn = np.abs(x_hp)\n    \n    peaks = local_maxima_1d_window(x_dn, window)\n    \n    heights = x_dn[peaks]\n    \n    ii = np.argsort(heights)[::-1]\n    \n    peaks = peaks[ii]\n    heights = heights[ii]\n    \n    ky = heights\n    kx = np.arange(1, heights.shape[0]+1)\n    \n    conv_length = 9\n\n    grad = np.diff(ky, 1)/np.diff(kx, 1)\n    grad = np.convolve(grad, np.ones(conv_length)/conv_length)#, mode='valid')\n    grad = grad[conv_length-1:-conv_length+1]\n    \n    knee_x = plateau_detection(grad, -0.01, plateau_length=1000)\n    knee_x -= conv_length//2\n    \n    if visualise:\n        plt.plot(grad, color=visualise_color)\n        plt.axvline(knee_x, ls=\"--\", color=visualise_color)\n    \n    peaks_x = peaks[:knee_x]\n    peaks_y = heights[:knee_x]\n    \n    ii = np.argsort(peaks_x)\n    peaks_x = peaks_x[ii]\n    peaks_y = peaks_y[ii]\n        \n    return peaks_x, peaks_y, x_hp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"53db94c29df4edecaa28bfb52413bb45bbd5047d"},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef clip(v, l, u):\n    \"\"\"Numba helper function to clip a value\"\"\"\n    \n    if v < l:\n        v = l\n    elif v > u:\n        v = u\n        \n    return v","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9a8084eb0eaaa0ee1515e55690ab83c6020d2bca"},"cell_type":"code","source":"peak_features_names = [\n    'ratio_next',\n    'ratio_prev',\n    'small_dist_to_min',\n    'sawtooth_rmse',\n]\n\nnum_peak_features = len(peak_features_names)\n\n@numba.jit(nopython=True)\ndef create_sawtooth_template(sawtooth_length, pre_length, post_length):\n    \"\"\"Generate sawtooth template\"\"\"\n    \n    l = pre_length+post_length+1\n    \n    st = np.zeros(l)\n    for i in range(sawtooth_length+1):\n        \n        j = pre_length+i\n        if j < l:\n            st[j] = 1 - ((2./sawtooth_length) * i)\n        \n    return st\n\n@numba.jit(nopython=True)\ndef calculate_peak_features(px, x_hp0, ws=5, wl=25):\n    \"\"\"\n    Calculate features for peaks.\n    Parameters\n    ----------\n    px : ndarray\n        Indices of peaks.\n    x_hp0 : ndarray\n        The array to search.\n    ws : np.int\n        How many points on each side to use for small window features\n    wl : np.int\n        How many points on each side to use for large window features\n    Returns\n    -------\n    features : ndarray\n        Features calculate for each peak in `x_hp0`.\n    \"\"\"\n    \n    features = np.ones((px.shape[0], num_peak_features), dtype=np.float64) * np.nan\n    \n    for i in range(px.shape[0]):\n        \n        feature_number = 0\n        \n        x = px[i]\n        x_next = x+1\n        x_prev = x-1\n        \n        h0 = x_hp0[x]\n\n        ws_s = clip(x-ws, 0, 800000-1)\n        ws_e = clip(x+ws, 0, 800000-1)\n        wl_s = clip(x-wl, 0, 800000-1)\n        wl_e = clip(x+wl, 0, 800000-1)\n        \n        ws_pre = x - ws_s\n        ws_post = ws_e - x\n        \n        wl_pre = x - wl_s\n        wl_post = wl_e - x\n        \n        if x_next < 800000:\n            h0_next = x_hp0[x_next]\n            features[i, feature_number] = np.abs(h0_next)/np.abs(h0)\n        feature_number += 1\n            \n        if x_prev >= 0:\n            h0_prev = x_hp0[x_prev]\n            features[i, feature_number] = np.abs(h0_prev)/np.abs(h0)\n        feature_number += 1\n            \n        x_hp_ws0 = x_hp0[ws_s:ws_e+1]\n        x_hp_wl0 = x_hp0[wl_s:wl_e+1]\n        x_hp_wl0_norm = (x_hp_wl0/np.abs(h0))\n        x_hp_ws0_norm = (x_hp_ws0/np.abs(h0))\n        x_hp_abs_wl0 = np.abs(x_hp_wl0)\n        wl_max_0 = np.max(x_hp_abs_wl0)\n        \n        ws_opp_peak_i = np.argmin(x_hp_ws0*np.sign(h0))\n        \n        features[i, feature_number] = ws_opp_peak_i - ws\n        feature_number += 1\n        \n        x_hp_wl0_norm_sign = x_hp_wl0_norm * np.sign(h0)\n        \n        sawtooth_length = 3\n        st = create_sawtooth_template(sawtooth_length, wl_pre, wl_post)\n        assert np.argmax(st) == np.argmax(x_hp_wl0_norm_sign)\n        assert st.shape[0] == x_hp_wl0_norm_sign.shape[0]\n        features[i, feature_number] = np.mean(np.power(x_hp_wl0_norm_sign - st, 2))\n        feature_number += 1\n        \n        if i == 0:\n            assert feature_number == num_peak_features\n        \n    return features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2d6a2ffc2c829cb70f01cbc3bb92fd9125c6cd93"},"cell_type":"code","source":"def process_signal(\n    data,\n    window=25,\n):\n    \"\"\"\n    Process a signal trace to find the peaks and calculate features for each peak.\n    Parameters\n    ----------\n    data : ndarray\n        The array to search.\n    window : np.int\n        How many points on each side to use for the local maxima test\n    Returns\n    -------\n    px0 : ndarray\n        Indices for each peak in `data`.\n    height0 : ndarray\n        Absolute heaight for each peak in `data`.\n    f0 : ndarray\n        Features calculate for each peak in `data`.\n    \"\"\"\n    \n    px0, height0, x_hp0 = get_peaks(\n        data.astype(np.float),\n        window=window, \n    )\n            \n    f0 = calculate_peak_features(px0, x_hp0)\n    \n    return px0, height0, f0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c9f2c506c728f4b3bddb50c2149326b09c2ab6f3"},"cell_type":"code","source":"def process_measurement_peaks(data, signal_ids):\n    \"\"\"\n    Process three signal traces in measurment to find the peaks\n    and calculate features for each peak.\n    Parameters\n    ----------\n    data : ndarray\n        Signal traces.\n    signal_ids : ndarray\n        Signal IDs for each of the signal traces in measurment\n    Returns\n    -------\n    res : ndarray\n        Data for each peak in the three traces in `data`.\n    sigid_res : ndarray\n        Signal ID for each row in `res`.\n    \"\"\"\n    res = []\n    sigid_res = []\n    \n    assert data.shape[1] % 3 == 0\n    N = data.shape[1]//3\n    \n    for i in range(N):\n        \n        sigids = signal_ids[i*3:(i+1)*3]\n        x = data[:, i*3:(i+1)*3].astype(np.float)\n        \n        px0, height0, f0 = process_signal(\n            x[:, 0]\n        )\n        \n        px1, height1, f1 = process_signal(\n            x[:, 1]\n        )\n        \n        px2, height2, f2 = process_signal(\n            x[:, 2]\n        )\n        \n        if px0.shape[0] != 0:\n            res.append(np.hstack([\n                px0[:, np.newaxis], \n                height0[:, np.newaxis],\n                f0,\n            ]))\n            \n            sigid_res.append(np.ones(px0.shape[0], dtype=np.int) * sigids[0])\n        \n        if px1.shape[0] != 0:\n            res.append(np.hstack([\n                px1[:, np.newaxis], \n                height1[:, np.newaxis],\n                f1,\n            ]))\n\n            sigid_res.append(np.ones(px1.shape[0], dtype=np.int) * sigids[1])\n        \n        if px2.shape[0] != 0:\n            res.append(np.hstack([\n                px2[:, np.newaxis], \n                height2[:, np.newaxis],\n                f2,\n            ]))\n\n            sigid_res.append(np.ones(px2.shape[0], dtype=np.int) * sigids[2])\n            \n    return res, sigid_res","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bfe983e9594a9d0a2a5c3b62dd9330bd3340a21d"},"cell_type":"code","source":"def process_measurement(data_df, meta_df, fft_data):\n    \"\"\"\n    Process three signal traces in measurment to find the peaks\n    and calculate features for each peak.\n    Parameters\n    ----------\n    data_df : pandas.DataFrame\n        Signal traces.\n    meta_df : pandas.DataFrame\n        Meta data for measurement\n    fft_data : ndarray\n        50Hz fourier coefficient for three traces\n    Returns\n    -------\n    peaks : pandas.DataFrame\n        Data for each peak in the three traces in `data`.\n    \"\"\"\n    peaks, sigids = process_measurement_peaks(\n        data_df.values, # [:, :100*3], \n        meta_df['signal_id'].values, # [:100*3]\n    )\n    \n    peaks = np.concatenate(peaks)\n\n    peaks = pd.DataFrame(\n        peaks,\n        columns=['px', 'height'] + peak_features_names\n    )\n    peaks['signal_id'] = np.concatenate(sigids)\n\n    # Calculate the phase resolved location of each peak\n\n    phase_50hz = np.angle(fft_data, deg=False) # fft_data[:, 1]\n\n    phase_50hz = pd.DataFrame(\n        phase_50hz,\n        columns=['phase_50hz']\n    )\n    phase_50hz['signal_id'] = meta_df['signal_id'].values\n\n    peaks = pd.merge(peaks, phase_50hz, on='signal_id', how='left')\n\n    dt = (20e-3/(800000))\n    f1 = 50\n    w1 = 2*np.pi*f1\n    peaks['phase_aligned_x'] = (np.degrees(\n        (w1*peaks['px'].values*dt) + peaks['phase_50hz'].values\n    ) + 90) % 360\n    \n    # Calculate the phase resolved quarter for each peak\n    peaks['Q'] = pd.cut(peaks['phase_aligned_x'], [0, 90, 180, 270, 360], labels=[0, 1, 2, 3])\n    \n    return peaks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"628b07cea495ea037af0a851ca1603fe8f351c44"},"cell_type":"code","source":"@numba.jit(nopython=True, parallel=True)\ndef calculate_50hz_fourier_coefficient(data):\n    \"\"\"Calculate the 50Hz Fourier coefficient of a signal.\n    Assumes the signal is 800000 data points long and covering 20ms.\n    \"\"\"\n\n    n = 800000\n    assert data.shape[0] == n\n    \n    omegas = np.exp(-2j * np.pi * np.arange(n) / n).reshape(n, 1)\n    m_ = omegas ** np.arange(1, 2)\n    \n    m = m_.flatten()\n\n    res = np.zeros(data.shape[1], dtype=m.dtype)\n    for i in numba.prange(data.shape[1]):\n        res[i] = m.dot(data[:, i].astype(m.dtype))\n            \n    return res","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f8b42f97ae88f41709c622f9674a4cf3ee055b57"},"cell_type":"markdown","source":"# Preprocessing Overview\n\nEach signal trace is preprocessed to identify the peaks and calculate features.\n\nThe steps performed to do this are:\n\n    1. Flatten the trace using the flatiron function\n    2. Identify the local maxima peaks in each trace using the local_maxima_1d_window function\n    3. Filter the peaks identified in step 2 to separate the signal from the noise\n    4. Calculate features for each peak identified in step 3 using calculate_peak_features"},{"metadata":{"_uuid":"55d238643ee258df8ce21cac06ff3d6fcd220fe0"},"cell_type":"markdown","source":"Step 2. To identify the local maxima the function local_maxima_1d_window is used. This function takes a window length argument, which is the number of points on each side to use for the comparison. An example of the behaviour of this function can be seen below:"},{"metadata":{"trusted":true,"_uuid":"d3c0620b07e745753a46b1eb5accdde6fa2193c6"},"cell_type":"code","source":"a = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 1, 0, 0, 3, 0, 0, 0, 0, 0, 1, 0, 2, 0, 3, 0, 2, 0, 1, 0])\n\np1 = local_maxima_1d_window(a, w=1)\np3 = local_maxima_1d_window(a, w=3)\np4 = local_maxima_1d_window(a, w=4)\n\nplt.plot(a, marker='o')\nplt.scatter(p1, a[p1]+0.2, color='red', label='1')\nplt.scatter(p3, a[p3]+0.4, color='orange', marker='x', label='3')\nplt.scatter(p4, a[p4]+0.6, color='grey', marker='^', label='4')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0b3f796a77aa19b8e7812828284f14c6fcb1cefd"},"cell_type":"markdown","source":"Step 3. Once all the peaks in a trace have been identified, the peaks caused by the noise in the signal need to be removed. This is performed in the get_peaks function. When the peaks are ordered by height, knee detection is performed to identify the point when the height of the peaks stops changing due to the noise floor being reached. The steps are:\n    \n    1. Order the peaks by their height\n    2. Calculate the gradient between each consecutive pair of peaks\n    3. Smooth the gradients using a convolution operation\n    4. Find the noise floor using the plateau_detection function\n    \nThe location of the threshold identified for a number of signals can be seen below:"},{"metadata":{"trusted":true,"_uuid":"1550200e68125d1ff0cd8c5b1c47b5af17a270e2"},"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n\nsigids = [2323, 10, 4200, 4225]\ncolours = ['blue', 'red', 'orange', 'grey']\nfor i, sigid in enumerate(sigids):\n    d = pq.read_pandas(\n        data_dir + '/train.parquet',\n        columns=[str(sigid)]\n    ).to_pandas().values[:, 0].astype(np.float)\n    get_peaks(d, visualise=True, visualise_color=colours[i])\n\nplt.xlim([0, 4000])\nplt.axhline(-0.01, color='black', ls='--')\nplt.yscale(\"symlog\")\nplt.xscale(\"symlog\")\n\nplt.xlabel('Sorted peak index')\nplt.ylabel('Gradient')\nplt.suptitle('Example of peak filtering')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"17c3623d58099b02f0e7ef1609e6c9c98311bd01"},"cell_type":"markdown","source":"Some examples of the peaks detected using the steps outlined above can be seen below:"},{"metadata":{"trusted":true,"_uuid":"a22ee7046d8ba12ff8bbfa5d554b41cd9ed98fa4"},"cell_type":"code","source":"sigids = [100, 230, 2323, 10, 4200, 4225]\n\nfor sigid in sigids:\n    d = pq.read_pandas(\n        data_dir + '/train.parquet',\n        columns=[str(sigid)]\n    ).to_pandas().values[:, 0].astype(np.float)\n\n    fig, ax = plt.subplots(1, 1, figsize=(12, 4))\n    plt.plot(d, alpha=0.75)\n\n    px, py, _ = get_peaks(d)\n    \n    plt.scatter(px, d[px], color=\"red\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b41c1c815501ebdc0a6c3a6860ced7fbf489f81f"},"cell_type":"markdown","source":"Step 4. Once the noisy peaks have been removed, features are calcuated for each of the remaining peaks. This is performed in calculate_peak_features. Of the features the most interesting is the sawtooth_rmse feature. This is the RMSE between the window of values including the peak and 25 data points either side of the peak and a sawtooth like template. An example of this can be seen below:"},{"metadata":{"trusted":true,"_uuid":"f066ab6089debca8a98887711fe6ca14773634c4"},"cell_type":"code","source":"def plot_sawtooth_example(measid, px, w):\n    \n    display(HTML(meta_train_df[meta_train_df['id_measurement'] == measid].to_html()))\n    \n    sigids = [measid * 3 + i for i in range(3)]\n    \n    fig, ax = plt.subplots(1, 1, figsize=(12, 4))\n\n    for i, sigid in enumerate(sigids):\n        d = pq.read_pandas(\n            data_dir + '/train.parquet',\n            columns=[str(sigid)]\n        ).to_pandas().values[:, 0].astype(np.float)\n\n        d = flatiron(d)\n        d = d[px-w:px+w+1]\n        plt.plot(d, marker='o', label='{}) {}'.format(i, sigid))\n\n        if i == 1:\n            s = d[w]\n            ft = create_sawtooth_template(3, w, w) * s\n            plt.plot(ft, color='black', ls='--')\n\n    plt.axvline(w, color='black', ls='--')\n    plt.legend()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c7803274875981ccd0736beab5a8ec045420c81e"},"cell_type":"code","source":"measid = 705\npx = 520128 # 93911 # 123144\nw = 10\n\nplot_sawtooth_example(measid, px, w)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c632a9b21e52c36c8e86d5ac05d949d5eeccb92d"},"cell_type":"code","source":"measid = 1091\npx = 593144\nw = 10\n\nplot_sawtooth_example(measid, px, w)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f10b050729d291d0de236a7d12fafb42d727505c"},"cell_type":"code","source":"measid = 2760\npx = 269640\nw = 10\n\nplot_sawtooth_example(measid, px, w)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"902e449045c1a8ae9620ba5a8da03827d6d20c96"},"cell_type":"markdown","source":"# Preprocess"},{"metadata":{},"cell_type":"markdown","source":"### Preprocess Train Data"},{"metadata":{"trusted":true,"_uuid":"21d283825e5c8236a309922088782d5868e4170f"},"cell_type":"code","source":"%%time\n\nif RECALCULATE_PEAK_FEATURES:\n\n    train_df = pq.read_pandas(\n        data_dir + '/train.parquet'\n    ).to_pandas()\n    \n    train_fft = calculate_50hz_fourier_coefficient(train_df.values)\n    \n    train_peaks = process_measurement(\n        train_df, \n        meta_train_df, \n        train_fft\n    )\n\n    del train_df, train_fft\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e633b3d3cd231af87d73b522910649e50a6014bd"},"cell_type":"code","source":"if not RECALCULATE_PEAK_FEATURES:\n    train_peaks = pd.read_pickle(PREPROCESSED_DIR + '/train_peaks.pkl')\n    \ntrain_peaks.to_pickle(\"train_peaks.pkl\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e20043531178e06c0911e08ba14717909397ce91"},"cell_type":"code","source":"train_peaks = pd.merge(train_peaks, meta_train_df[['signal_id', 'id_measurement', 'target']], on='signal_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3c781ad4f13156d4fec61d59071d3f5b21684f48"},"cell_type":"code","source":"train_peaks.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"88792d4a0c98895da1cceb7f1b264e7fb185db8d"},"cell_type":"code","source":"train_peaks.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Preprocess Test Data"},{"metadata":{"trusted":true,"_uuid":"583e407ed4f6261244ea31dc96add2afe671c9fe"},"cell_type":"code","source":"%%time\n\nif RECALCULATE_PEAK_FEATURES:\n    \n    NUM_TEST_CHUNKS = 10\n    \n    test_chunk_size = int(np.ceil((meta_test_df.shape[0]/3.)/float(NUM_TEST_CHUNKS))*3.)\n    \n    test_peaks = []\n\n    for j in range(NUM_TEST_CHUNKS):\n\n        j_start = j*test_chunk_size\n        j_end = (j+1)*test_chunk_size\n        \n        signal_ids = meta_test_df['signal_id'].values[j_start:j_end]\n\n        test_df = pq.read_pandas(\n            data_dir + '/test.parquet',\n            columns=[str(c) for c in signal_ids]\n        ).to_pandas()\n\n        test_fft = calculate_50hz_fourier_coefficient(test_df.values)\n        \n        p = process_measurement(\n            test_df, \n            meta_test_df.iloc[j_start:j_end], \n            test_fft\n        )\n\n        test_peaks.append(p)\n        \n        print(j)\n        \n        del test_df\n        gc.collect()\n        \n    del test_fft\n    gc.collect()\n        \n    test_peaks = pd.concat(test_peaks)\n    \n# Wall time: 9min 45s","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7e117ad6050ccacc07162349ed7bb7fc26803a13"},"cell_type":"code","source":"if not RECALCULATE_PEAK_FEATURES:\n    test_peaks = pd.read_pickle(PREPROCESSED_DIR + '/test_peaks.pkl')\n    \ntest_peaks.to_pickle(\"test_peaks.pkl\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f638e2ec57620c8a99890659fcb3ab3653834c6f"},"cell_type":"code","source":"test_peaks = pd.merge(test_peaks, meta_test_df[['signal_id', 'id_measurement']], on='signal_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"85b49d1acb963ecc67c2516fc7c056e981a0e950"},"cell_type":"code","source":"test_peaks.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c91c0f92609e6b2584f63c56c4f299e62b68fd67"},"cell_type":"code","source":"test_peaks.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"75a44026079a055e8043160861f2fe1e9e5962ed"},"cell_type":"code","source":"meta_test_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2c02c459d14e9dda0dee4a84c404fa89f2967dac"},"cell_type":"code","source":"test_peaks.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"64f7d449917916b670620b9ede6e222eee0a57ba"},"cell_type":"markdown","source":"# Features"},{"metadata":{"trusted":true,"_uuid":"a20dee063467672f86e388592cf721cde9b7f543"},"cell_type":"code","source":"def save_importances(importances_, filename='importances.png', print_results=False, figsize=(8, 6)):\n    mean_gain = importances_[['gain', 'feature']].groupby('feature').mean()\n    if print_results:\n        display(HTML(mean_gain.sort_values('gain').to_html()))\n    importances_['mean_gain'] = importances_['feature'].map(mean_gain['gain'])\n    plt.figure(figsize=figsize)\n    sns.barplot(x='gain', y='feature', data=importances_.sort_values('mean_gain', ascending=False))\n    plt.tight_layout()\n    plt.savefig(filename)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"96038328c89ab19cdceafb162b6d4a433d5c644f"},"cell_type":"code","source":"def process(peaks_df, meta_df, use_improved=True):\n\n    results = pd.DataFrame(index=meta_df['id_measurement'].unique())\n    results.index.rename('id_measurement', inplace=True)\n    \n    ################################################################################\n    \n    if not USE_SIMPLIFIED_VERSION:\n        # Filter peaks using ratio_next and height features\n        # Note: may not be all that important\n        peaks_df = peaks_df[~(\n            (peaks_df['ratio_next'] > 0.33333)\n            & (peaks_df['height'] > 50)\n        )]\n    \n    ################################################################################\n\n    # Count peaks in phase resolved quarters 0 and 2\n    p = peaks_df[peaks_df['Q'].isin([0, 2])].copy()\n    res = p.groupby('id_measurement').agg(\n    {\n        'px': ['count'],\n    })\n    res.columns = [\"peak_count_Q02\"]\n    results = pd.merge(results, res, on='id_measurement', how='left')\n        \n    ################################################################################\n    \n    # Count total peaks for each measurement id\n    res = peaks_df.groupby('id_measurement').agg(\n    {\n        'px': ['count'],\n    })\n    res.columns = [\"peak_count_total\"]\n    results = pd.merge(results, res, on='id_measurement', how='left')\n\n    ################################################################################\n\n    # Count peaks in phase resolved quarters 1 and 3\n    p = peaks_df[peaks_df['Q'].isin([1, 3])].copy()\n    res = p.groupby('id_measurement').agg(\n    {\n        'px': ['count'],\n    })\n    res.columns = ['peak_count_Q13']\n    results = pd.merge(results, res, on='id_measurement', how='left')\n    \n    ################################################################################\n    \n    # Calculate additional features using phase resolved quarters 0 and 2\n    \n    feature_quarters = [0, 2]\n    \n    p = peaks_df[peaks_df['Q'].isin(feature_quarters)].copy()\n    \n    p['abs_small_dist_to_min'] = np.abs(p['small_dist_to_min'])\n    \n    res = p.groupby('id_measurement').agg(\n    {\n        \n        'height': ['mean', 'std'],\n        'ratio_prev': ['mean'],\n        'ratio_next': ['mean'],\n        'abs_small_dist_to_min': ['mean'],        \n        'sawtooth_rmse': ['mean'],\n    })\n    res.columns = [\"_\".join(f) + '_Q02' for f in res.columns]     \n    results = pd.merge(results, res, on='id_measurement', how='left')\n        \n    return results","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8d516509a457ac118c43852c83c932a77c7a987a"},"cell_type":"code","source":"def create_features(\n    meta_df, \n    peaks_df,\n):\n    \n    res = process(peaks_df, meta_df)\n    \n    return res","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create Train features"},{"metadata":{"trusted":true,"_uuid":"960cd8e8da2212a37c8bb1f7fcc9952007f414d9"},"cell_type":"code","source":"%%time\nX_train_full = create_features( \n    meta_train_df, \n    train_peaks,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d1d3ba0948bf4bb0fda1452edd8504e753e7ae17"},"cell_type":"code","source":"y_train = (\n    meta_train_df.groupby('id_measurement')['target'].sum().round(0).astype(np.int)\n    != 0\n).astype(np.float)\n\nassert np.all(y_train.index.values == X_train_full.index.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dbf9832f0bab500dd31c5dffdc0f6b55d5be5636"},"cell_type":"code","source":"X_train_full.head(6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fc9c8b4345b6ba9325b24df4ae0d4a128e68bf50"},"cell_type":"code","source":"y_train.head(6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"653f652d0bfac22e619512950351c0e0ce665248"},"cell_type":"code","source":"y_train.sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ef9e37e1475da96f8b9394ce8f3e739523f98146"},"cell_type":"code","source":"X_train_full.shape, y_train.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create Test Features"},{"metadata":{"trusted":true,"_uuid":"5b1e07d459e3c17ca387fb06574abb480e4a682a"},"cell_type":"code","source":"%%time\nX_test_full = create_features( \n    meta_test_df, \n    test_peaks,\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4d153f7ecb932d9b4b5b5cb22eccbb856341c0a5"},"cell_type":"code","source":"X_test_full.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bd0a8834b5be918319a5229a04d9cadd5992c850"},"cell_type":"code","source":"X_test_full.tail()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d22cd660eff11c384dabfc3cb2d4d288faf68fc4"},"cell_type":"markdown","source":"# Drop Features"},{"metadata":{"trusted":true,"_uuid":"00f6cf66b3f72cbf4ec6c3946bed3f96ed17ba58"},"cell_type":"code","source":"default_features_to_drop = [\n    'id_measurement',\n    'phase',\n    'fft_phase',\n    'signal_id',\n]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a581519a47bda6e44622e743ec55b5854f5ae0b"},"cell_type":"code","source":"if USE_SIMPLIFIED_VERSION:\n    additional_features_to_drop = [\n        #'peak_count_Q13',\n        'abs_small_dist_to_min_mean_Q02',\n        'height_mean_Q02',\n        #'height_std_Q02',\n        #'sawtooth_rmse_mean_Q02',\n        #'peak_count_Q02',\n        'ratio_next_mean_Q02',\n        'ratio_prev_mean_Q02',\n        #'peak_count_total' \n    ]\nelse:\n    additional_features_to_drop = []","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b9b5ea7b62d43c774ec3438e676100ef95e097fa"},"cell_type":"code","source":"features_to_drop = (\n    default_features_to_drop\n    + additional_features_to_drop \n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8b5b757aa8a709678a6227cfd8919dabea6ada8f"},"cell_type":"code","source":"if USE_SIMPLIFIED_VERSION:\n    feature_names = sorted([c for c in X_train_full.columns if c not in features_to_drop])\nelse:\n    # to maintain consistent order with original version after feature renaming\n    feature_names = [\n        'peak_count_Q13',\n        'abs_small_dist_to_min_mean_Q02',\n        'height_mean_Q02',\n        'height_std_Q02',\n        'sawtooth_rmse_mean_Q02',\n        'peak_count_Q02',\n        'ratio_next_mean_Q02',\n        'ratio_prev_mean_Q02',\n        'peak_count_total'   \n    ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"283120e97789533e03c5235a8018732aabd69dc7"},"cell_type":"code","source":"X_train = X_train_full[feature_names]\nX_test = X_test_full[feature_names]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"61b6fb0832d49d44af741407f9dedaf3cf8cf8d7"},"cell_type":"code","source":"X_train.head(6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"18fdc58675849577700c26e582c1e8cc5e90d51b"},"cell_type":"code","source":"X_train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f355c17127e6f8e4eac73f104e0100e006970e3f"},"cell_type":"code","source":"feature_names","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"61ce738f08ad86e7fbe77ef7a6dd62772e210b8f"},"cell_type":"markdown","source":"# Training Code"},{"metadata":{"trusted":true,"_uuid":"57420da36aa36d13e03967d31daaf065b41a8fb1"},"cell_type":"code","source":"@numba.jit\ndef mcc(tp, tn, fp, fn):\n    sup = tp * tn - fp * fn\n    inf = (tp + fp) * (tp + fn) * (tn + fp) * (tn + fn)\n    if inf==0:\n        return 0\n    else:\n        return sup / np.sqrt(inf)\n\n#@numba.jit\ndef eval_mcc(y_true, y_prob, show=False):\n    \"\"\"\n    A fast implementation of Anokas mcc optimization code.\n\n    This code takes as input probabilities, and selects the threshold that \n    yields the best MCC score. It is efficient enough to be used as a \n    custom evaluation function in xgboost\n    \n    Source: https://www.kaggle.com/cpmpml/optimizing-probabilities-for-best-mcc\n    Source: https://www.kaggle.com/c/bosch-production-line-performance/forums/t/22917/optimising-probabilities-binary-prediction-script\n    Creator: CPMP\n    \"\"\"\n    idx = np.argsort(y_prob)\n    y_true_sort = y_true[idx]\n    n = y_true.shape[0]\n    nump = 1.0 * np.sum(y_true) # number of positive\n    numn = n - nump # number of negative\n    tp = nump\n    tn = 0.0\n    fp = numn\n    fn = 0.0\n    best_mcc = 0.0\n    best_id = -1\n    prev_proba = -1\n    best_proba = -1\n    mccs = np.zeros(n)\n    for i in range(n):\n        # all items with idx < i are predicted negative while others are predicted positive\n        # only evaluate mcc when probability changes\n        proba = y_prob[idx[i]]\n        if proba != prev_proba:\n            prev_proba = proba\n            new_mcc = mcc(tp, tn, fp, fn)\n            if new_mcc >= best_mcc:\n                best_mcc = new_mcc\n                best_id = i\n                best_proba = proba\n        mccs[i] = new_mcc\n        if y_true_sort[i] == 1:\n            tp -= 1.0\n            fn += 1.0\n        else:\n            fp -= 1.0\n            tn += 1.0\n    if show:\n        y_pred = (y_prob >= best_proba).astype(int)\n        score = matthews_corrcoef(y_true, y_pred)\n        plt.plot(mccs)\n        return best_proba, best_mcc, y_pred\n    else:\n        return best_proba,  best_mcc, None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"88fc4a2f8e4bafe79fa075aa2738efb89a78d7ad"},"cell_type":"code","source":"def train(\n    X_data, \n    y_data,\n    params,\n    feature_names, \n    verbose_eval=50,\n    kfold_random_state=23948,\n    predict=False,\n    calculate_mcc=False,\n    early_stopping_rounds=10,\n    num_folds=5,\n    num_iterations=10,\n):\n    \n    params = copy.deepcopy(params)\n    \n    models = []\n    cv_scores = []\n    val_cv_scores = []\n    if calculate_mcc:\n        mcc_scores = []\n        val_mcc_scores = []\n        thresholds=[]\n    else:\n        mcc_scores = None\n        val_mcc_scores = None\n        thresholds = None\n        \n    if predict:\n        yp_train = np.zeros(X_train.shape[0])\n        yp_val = np.zeros(X_train.shape[0])\n        yp_test = np.zeros(X_train.shape[0])\n    else:\n        yp_train = None\n        yp_val = None\n        yp_test = None\n        \n    for iter in range(num_iterations):\n        \n        np.random.seed(kfold_random_state + iter)\n        splits = np.zeros(X_data.shape[0], dtype=np.int)\n        m = y_data == 1\n        splits[m] = np.random.randint(0, 5, size=m.sum())\n        m = y_data == 0\n        splits[m] = np.random.randint(0, 5, size=m.sum())\n\n        for fold in range(num_folds):\n\n            val_fold = fold\n            test_fold = (fold + 1) % num_folds\n\n            train_folds = [f for f in range(num_folds) if f not in [val_fold, test_fold]]\n\n            train_indices = np.where(np.isin(splits, train_folds))[0]\n            val_indices = np.where(splits == val_fold)[0]\n            test_indices = np.where(splits == test_fold)[0]\n\n            trn = lgb.Dataset(\n                X_data.values[train_indices],\n                y_data[train_indices],\n                feature_name=feature_names,\n            )\n\n            val = lgb.Dataset(\n                X_data.values[val_indices],\n                y_data[val_indices],\n                feature_name=feature_names,\n            )\n\n            test = lgb.Dataset(\n                X_data.values[test_indices],\n                y_data[test_indices],\n                feature_name=feature_names,\n            )\n\n            # train model\n\n            with warnings.catch_warnings():\n                warnings.simplefilter(\"ignore\")\n                model = lgb.train(\n                    params, \n                    trn, \n                    #num_boost_round=10000, \n                    valid_sets=(trn, test, val), \n                    valid_names=(\"train\", \"test\", \"validation\"), \n                    fobj=None, \n                    feval=None,\n                    early_stopping_rounds=early_stopping_rounds,\n                    #evals_result=training_data,\n                    verbose_eval=verbose_eval,\n                )\n\n\n            yp = None\n\n            if predict:\n                yp = model.predict(X_train.values[train_indices])\n                yp_train[train_indices] += yp\n\n                yp_val_fold = model.predict(X_train.values[val_indices])\n                yp_val[val_indices] += yp_val_fold\n                \n                yp_test_fold = model.predict(X_train.values[test_indices])\n                yp_test[test_indices] += yp_test_fold\n\n            models.append(model)\n\n            s = model.best_score['test']\n            cv_scores.append(s[list(s.keys())[0]])\n\n            s = model.best_score['validation']\n            val_cv_scores.append(s[list(s.keys())[0]])\n\n            if calculate_mcc:\n\n                # find optimal probability threshold using training data\n                best_proba, _, _ = eval_mcc(\n                    y_train.values[train_indices].astype(np.float),\n                    yp, \n                    False\n                )\n\n                # calculate validation MCC score using optimal probability threshold\n                best_mcc_val = matthews_corrcoef(\n                    y_train.values[val_indices].astype(np.float), \n                    yp_val_fold.astype(np.float) > best_proba\n                )\n\n                # calculate test MCC score using optimal probability threshold\n                best_mcc = matthews_corrcoef(\n                    y_train.values[test_indices].astype(np.float), \n                    yp_test_fold.astype(np.float) > best_proba\n                )\n\n                val_mcc_scores.append(best_mcc_val)\n                mcc_scores.append(best_mcc)\n                thresholds.append(best_proba)\n\n    if yp_train is not None:\n        yp_train /= ((num_folds - 2) * num_iterations)\n        yp_val /= num_iterations\n        yp_test /= num_iterations\n            \n    cv_scores = np.array(cv_scores)\n    val_cv_scores = np.array(val_cv_scores)\n    if calculate_mcc:\n        mcc_scores = np.array(mcc_scores)\n        val_mcc_scores = np.array(val_mcc_scores)\n        thresholds = np.array(thresholds)\n        \n    return (\n        models,\n        cv_scores, val_cv_scores, \n        mcc_scores, val_mcc_scores, \n        yp_train, yp_val, yp_test, thresholds\n    )","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a10df827027e281e3aac130e26fd51228bff01d4"},"cell_type":"markdown","source":"# Train"},{"metadata":{"trusted":true,"_uuid":"37985a51764c3712cddf2f5049280081df441e83"},"cell_type":"code","source":"params = {\n    'objective': 'binary',\n    #'is_unbalance': True,\n    'boosting': 'gbdt',\n    'learning_rate': 0.01,\n    'num_leaves': 80, # 61,\n    'num_threads': 4,\n    'metric': 'binary_logloss',\n    'feature_fraction': 0.8,\n    'bagging_freq': 1,\n    'bagging_fraction': 0.8,\n    'seed': 23974,\n    'num_boost_round': 10000,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d368ffd58032cdd9629b21fbf8622c0766782b5e","scrolled":false},"cell_type":"code","source":"%%time\n\n(\n    models,\n    cv_scores, val_cv_scores, \n    mcc_scores, val_mcc_scores, \n    yp_train, yp_val, yp_test, thresholds\n) = train(\n    X_train, \n    y_train,\n    params,\n    feature_names, \n    verbose_eval=0, # 50,\n    kfold_random_state=123948,\n    predict=True,\n    calculate_mcc=True,\n    early_stopping_rounds=100,\n    num_iterations=25 # 100\n)\n\nprint(\"\")\nprint(\"CV Logloss: {:.4f} +/- {:.4f} ({:.4f})\".format(cv_scores.mean(), cv_scores.std()/np.sqrt(cv_scores.shape[0]), cv_scores.std()))\nprint(\"CV Val Logloss: {:.4f} +/- {:.4f} ({:.4f})\".format(val_cv_scores.mean(), val_cv_scores.std()/np.sqrt(val_cv_scores.shape[0]), val_cv_scores.std()))\nprint(\"CV MCC: {:.4f} +/- {:.4f} ({:.4f})\".format(mcc_scores.mean(), mcc_scores.std()/np.sqrt(mcc_scores.shape[0]), mcc_scores.std()))\nprint(\"CV VAL MCC: {:.4f} +/- {:.4f} ({:.4f})\".format(val_mcc_scores.mean(), val_mcc_scores.std()/np.sqrt(val_mcc_scores.shape[0]), val_mcc_scores.std()))\nprint(\"Threshold {:.4f} +/- {:.4f}\".format(thresholds.mean(), thresholds.std()/np.sqrt(thresholds.shape[0])))\n\n##################################\n\nbest_proba_meas, _, _ = eval_mcc(\n    y_train.values.astype(np.float),\n    yp_train, \n    False\n)\n\nbest_mcc_meas = matthews_corrcoef(\n    y_train.values.astype(np.float), \n    yp_val > best_proba_meas\n)\n\nbest_mcc_meas_test = matthews_corrcoef(\n    y_train.values.astype(np.float), \n    yp_test > best_proba_meas\n)\n\nplt.show()\n\n##################################\n\nyp_train_df = pd.DataFrame(\n    yp_train,\n    index=X_train_full.index\n)\n\nyp_val_df = pd.DataFrame(\n    yp_val,\n    index=X_train_full.index\n)\n\nyp_test_df = pd.DataFrame(\n    yp_test,\n    index=X_train_full.index\n)\n\ntrain_pred = meta_train_df[['id_measurement', 'signal_id', 'target']].copy()\ntrain_pred = pd.merge(train_pred, yp_train_df, on='id_measurement')\ntrain_pred.rename({0:'prediction'}, axis=1, inplace=True)\n\nval_pred = meta_train_df[['id_measurement', 'signal_id', 'target']].copy()\nval_pred = pd.merge(val_pred, yp_val_df, on='id_measurement')\nval_pred.rename({0:'prediction'}, axis=1, inplace=True)\n\ntest_pred = meta_train_df[['id_measurement', 'signal_id', 'target']].copy()\ntest_pred = pd.merge(test_pred, yp_test_df, on='id_measurement')\ntest_pred.rename({0:'prediction'}, axis=1, inplace=True)\n\nbest_proba, _, _ = eval_mcc(\n    train_pred['target'].values.astype(np.float),\n    train_pred['prediction'].values.astype(np.float),\n    False\n)\n\nbest_mcc = matthews_corrcoef(\n    val_pred['target'].values.astype(np.float), \n    val_pred['prediction'].values.astype(np.float) > best_proba\n)\n\nbest_mcc_test = matthews_corrcoef(\n    test_pred['target'].values.astype(np.float), \n    test_pred['prediction'].values.astype(np.float) > best_proba\n)\n\n##################################\n\nplt.show()\n\nprint(\"MCC (measurement): {:.3f}\".format(best_mcc_meas))\nprint(\"MCC Test (measurement): {:.3f}\".format(best_mcc_meas_test))\nprint(\"Best Probability Threshold(measurement): {:.3f}\".format(best_proba_meas))\n\nprint(\"MCC: {:.3f}\".format(best_mcc))\nprint(\"MCC Test: {:.3f}\".format(best_mcc_test))\nprint(\"Best Probability Threshold: {:.3f}\".format(best_proba))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f2c768d9f9d121d18c4b6522280e03b551b8aa9b"},"cell_type":"code","source":"_, bins, _ = plt.hist(mcc_scores, alpha=0.5, label='Test fold')\nplt.hist(val_mcc_scores, alpha=0.5, bins=bins, label='Validation fold')\nplt.legend()\nplt.suptitle(\"Distributions of validation and test fold MCC scores\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f4f8ed17dc5d5a99265e2f37bac33d2791ba797a"},"cell_type":"code","source":"thresholds = np.linspace(0, 1, 100)[1:]\n\nscores_train = []\nscores_val = []\nscores_test = []\n\nfor t in thresholds:\n    \n    s_train = matthews_corrcoef(\n        y_train.values.astype(np.float), \n        yp_train > t\n    )\n    \n    s_val = matthews_corrcoef(\n        y_train.values.astype(np.float), \n        yp_val > t\n    )\n    \n    s_test = matthews_corrcoef(\n        y_train.values.astype(np.float), \n        yp_test > t\n    )\n    \n    scores_train.append(s_train)\n    scores_val.append(s_val)\n    scores_test.append(s_test)\n    \nplt.plot(thresholds, scores_train)\nplt.plot(thresholds, scores_val)\nplt.plot(thresholds, scores_test)\nplt.axvline(best_proba)\nplt.axvline(thresholds[np.argmax(scores_val)], ls='--')\nplt.show()\n\nprint(np.max(scores_val), thresholds[np.argmax(scores_val)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fac3d1d15dbd18407d8a51b2bd38ef583c9be5c3"},"cell_type":"code","source":"train_pred['probability_thresholded'] = (train_pred['prediction'] > best_proba).astype(np.int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"96acdc644f65e36a16117b6d4172b047e17e7c59"},"cell_type":"code","source":"train_pred.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1fd8b5db500249a16db754d728e2490ec294300a"},"cell_type":"code","source":"train_pred.to_csv('train_results.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6a4b64bd998a66b6f112185bc8dfd14747262d3c"},"cell_type":"markdown","source":"# Analysis"},{"metadata":{"_uuid":"6758642eba4a8625c6408d1f04a6bfdf82ab78eb"},"cell_type":"markdown","source":"## Feature Importance"},{"metadata":{"trusted":true,"_uuid":"f4def295428dfc5e5e93cac825453a9975fe3b09"},"cell_type":"code","source":"importances = pd.DataFrame()\n\nfor fold_ in range(len(models)):\n    \n    model = models[fold_]\n    \n    imp_df = pd.DataFrame()\n    imp_df['feature'] = X_train.columns\n    imp_df['gain'] = model.feature_importance('gain')\n    imp_df['fold'] = fold_ + 1\n    importances = pd.concat([importances, imp_df], axis=0, sort=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"7d7826fd67dc5c2146f16362d7fe7d5e8bc08f95"},"cell_type":"code","source":"save_importances(importances_=importances, print_results=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d3d53edefa3f3cbb049f41f37627a2ed2afeb89a"},"cell_type":"markdown","source":"## Prediction Errors"},{"metadata":{"trusted":true,"_uuid":"8c9df4f1af1791f151eab2b512ef0e42bd30e6f6"},"cell_type":"code","source":"prediction_errors = y_train.copy().to_frame()\nprediction_errors.columns = ['target']\nprediction_errors['error'] = np.abs(y_train.values - yp_train)\n\nprediction_errors = prediction_errors.sort_values('error')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"554a8fa985e1a280baf9b4a61307637bad0b0fe5"},"cell_type":"code","source":"prediction_errors.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"b6a7c20f0ca6072fc2076e14d93f6c7cacaa020d"},"cell_type":"code","source":"prediction_errors.tail(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"82f9703d9dd4c3f8f91633b45c6dbdaa3495b826"},"cell_type":"code","source":"prediction_errors[prediction_errors['target'] == 0].tail(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25384529820317eed192cd8ce95afc4925518bd5"},"cell_type":"markdown","source":"# Test"},{"metadata":{"trusted":true,"_uuid":"e255ceaa75b1e3f46534d4c7e76da0b40a7dedbe"},"cell_type":"code","source":"X_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b8af5eaccb5cdeb22178ac66c49054753504823e"},"cell_type":"code","source":"X_test.shape, meta_test_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3a08a6e94a80d599ea07ee4e5090507a01dd20de"},"cell_type":"code","source":"len(models)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e65b5e3acf2eb0441b4241e9ca854f6f56d52a63"},"cell_type":"code","source":"%%time\n\nyp_test = np.zeros(X_test.shape[0])\n\nfor j in range(len(models)):\n    \n    model = models[j]\n\n    yp_test += model.predict(X_test.values)/len(models)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d727fcb15d1c53860ef76d8a0fc5e55cc5cfa359"},"cell_type":"code","source":"test_submission = pd.DataFrame(\n    yp_test,\n    index=X_test.index,\n    columns=['probability']\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ada62bc51fa8aca157f4ec31c5e286b83379ff98"},"cell_type":"code","source":"test_submission['probability_thresholded'] = (yp_test > best_proba).astype(np.int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"526814fff119a23bc61ec6dd0e03c32d83bfd464"},"cell_type":"code","source":"test_submission = pd.merge(\n    meta_test_df[['id_measurement', 'signal_id']],\n    test_submission,\n    on='id_measurement',\n    how='left'\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dcd059a6c11ab4cfba88bb95f75bcd3e09325880"},"cell_type":"code","source":"test_submission['probability_thresholded'].sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f7281eef92233274e89083b1ff12f77605d70762"},"cell_type":"code","source":"test_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2eb8baa3c69dcaf0f51c2036634191a953c60ba7"},"cell_type":"code","source":"test_submission.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b0ad9f62620fb5d936de68c1ee57534f9e319269"},"cell_type":"code","source":"test_submission.to_csv('test_results.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2ec56637878c89b5b71eea02512e5533af381b1e"},"cell_type":"code","source":"submission = pd.read_csv(data_dir + '/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a5ab48fd1e62456973e010124647c62bff47a684"},"cell_type":"code","source":"assert np.all(submission['signal_id'].values == test_submission['signal_id'].values)\n\nsubmission['target'] = test_submission['probability_thresholded'].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"79e04bbcf4c37009953525921a293ac6f3eef185"},"cell_type":"code","source":"submission.shape, X_test.shape[0]*3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b21e9eb326752370bf214f29fbcff3f028a56cbb"},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7af68e150e86e5db1b61e1c133db4dadd61ecd04"},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6b01c5f24c9d381c609189442ed01012d6c81f2b"},"cell_type":"code","source":"!head submission.csv","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"835f75d6e42ba5a58098452b9264f5c483e6d639"},"cell_type":"markdown","source":"# Analysis"},{"metadata":{},"cell_type":"markdown","source":"### Proportions of train and test data detected as faulty"},{"metadata":{"trusted":true,"_uuid":"c98b1b52a8a2b99bec944a3db3c838ebdcefc9dc"},"cell_type":"code","source":"submission['target'].sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"544cdd96de2d84a108d80de810044b3362687663"},"cell_type":"code","source":"meta_train_df['target'].mean(), (yp_train > best_proba).astype(np.float).mean(), submission['target'].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"46a71d517ff1f512992ce949e919be56253880ec"},"cell_type":"code","source":"(meta_train_df['target'].values.reshape(-1, 3).sum(axis=1) > 0).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bcecf6cc3133878a323de026104072e549a7bfa1"},"cell_type":"code","source":"(submission['target'].values.reshape(-1, 3).sum(axis=1) > 0).mean()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9cfd0245154cb9510cdf77c21278d43d86842f5d"},"cell_type":"markdown","source":"## Important Features Partial Plots"},{"metadata":{"trusted":true,"_uuid":"2a6b041f766625d881ad60911264bb5362189dee"},"cell_type":"code","source":"X_train_full.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"221b91d386110adf65553afd1f72c293ee2cb8aa"},"cell_type":"code","source":"y_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"6485f59f04de3f414e80e43ba2e47089f72f7fd9"},"cell_type":"code","source":"sns.set_context(\"paper\", font_scale=2)\n\nimportant_features = importances[['gain', 'feature']].groupby('feature').mean().sort_values('gain').index.values[::-1]\n\nfor f in important_features:\n    fig, ax = plt.subplots(1, 1, figsize=(8,6))\n    sns.regplot(\n        f,\n        'target',\n        pd.merge(X_train, y_train.to_frame(), on='id_measurement', how='left'),\n        logistic=True,\n        n_boot=100,\n        y_jitter=.1,\n        scatter_kws={'alpha':0.1, 'edgecolor':'none'},\n        ax=ax\n    )\n    plt.show()\n    #break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"03e74c9208d0ff82806c487a37025fca710b5b78"},"cell_type":"markdown","source":"# Save Models"},{"metadata":{"trusted":true,"_uuid":"839b746dbec3187a38771abd0449eb2d1a5a6f67"},"cell_type":"code","source":"for j in range(len(models)):\n    models[j].save_model('model_{}.txt'.format(j))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cbcffc3fb464562f087b63a55739b7f860bc9157"},"cell_type":"code","source":"!head model_0.txt","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}