{"cells":[{"metadata":{},"cell_type":"markdown","source":"**An example of a data science project.**\n\nHere I have essentially reproduced: **https://www.kaggle.com/mark4h/vsb-1st-place-solution** and added a few comments.\n\n\n\nOther kernels I have used: \n* https://www.kaggle.com/genericurl/basic-eda\n* https://www.kaggle.com/miklgr500/flatiron"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport numba   # JIT compiler for python\nimport matplotlib.pyplot as plt  # graphics\nimport lightgbm as lgb  # Gradient boosting\nimport scipy.stats  # stats\nimport gc   # Garabage collector\nfrom sklearn import metrics\nimport seaborn as sns\n\nRANDOM_SEED = 27\n\ndata_dir = '../input/vsb-power-line-fault-detection'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load and explore metadata"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","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},"cell_type":"code","source":"meta_train_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df.target.value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df.id_measurement.nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df.groupby('id_measurement').target.sum().value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df[meta_train_df.target==0].sample(5, random_state=RANDOM_SEED)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_train_df[meta_train_df.target==1].sample(5, random_state=RANDOM_SEED)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load and explore training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_parquet(data_dir + '/train.parquet')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"negative_signal_ids = meta_train_df[meta_train_df.id_measurement==1287].signal_id.values\npositive_signal_ids = meta_train_df[meta_train_df.id_measurement==2649].signal_id.values\n\nnegative_sample = train_df.iloc[:, negative_signal_ids].values\npositive_sample = train_df.iloc[:, positive_signal_ids].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(18, 4))\nplt.title('Normal powerline')\nplt.plot(negative_sample, alpha=0.8);\n\nplt.figure(figsize=(18, 4))\nplt.title('Faulty powerline')\nplt.plot(positive_sample, alpha=0.8);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Preprocessing overview"},{"metadata":{},"cell_type":"markdown","source":"In order to extract features each signal has to be processed first.\n\nThis is done in 4 steps:\n    \n1. Flatten signal using EMA residuals\n2. Identify local maxima\n3. Filter the peaks to separate signal from noise\n4. Transform scale"},{"metadata":{},"cell_type":"markdown","source":"## Step 1. Flatten signal"},{"metadata":{},"cell_type":"markdown","source":"Signals are flattened by calculating exponential moving average (EMA, https://en.wikipedia.org/wiki/Exponential_smoothing)\nand only keeping the difference between EMA and the actual signal."},{"metadata":{"trusted":true},"cell_type":"code","source":"assert numba.__version__ == '0.46.0'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef ema_residuals(x, alpha=0.01):\n    \"\"\"\n    Flatten signal\n    Based on: https://www.kaggle.com/miklgr500/flatiron\n    \"\"\"\n    new_x = np.zeros_like(x)\n    ema = x[0]\n    for i in range(1, len(x)):\n        ema = ema*(1-alpha) + alpha*x[i]\n        new_x[i] = x[i] - ema\n    return new_x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"@numba.jit(nopython=True)\ndef ema(x, alpha=0.01):\n    \"\"\"\n    Flatten signal\n    Based on: https://www.kaggle.com/miklgr500/flatiron\n    \"\"\"\n    new_x = np.zeros_like(x)\n    ema = x[0]\n    for i in range(1, len(x)):\n        ema = ema*(1-alpha) + alpha*x[i]\n        new_x[i] = ema\n    return new_x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"flat_negative_sample = np.zeros_like(negative_sample)\nflat_positive_sample = np.zeros_like(positive_sample)\n\nfor i in range(3):\n    flat_negative_sample[:,i] = ema_residuals(negative_sample[:,i])\n    flat_positive_sample[:,i] = ema_residuals(positive_sample[:,i])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(18, 4))\nplt.title('Normal powerline')\nplt.plot(flat_negative_sample, alpha=0.8);\n\nplt.figure(figsize=(18, 4))\nplt.title('Faulty powerline')\nplt.plot(flat_positive_sample, alpha=0.8);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"flat_negative_sample = np.zeros_like(negative_sample)\nflat_positive_sample = np.zeros_like(positive_sample)\n\nfor i in range(3):\n    flat_negative_sample[:,i] = ema(negative_sample[:,i])\n    flat_positive_sample[:,i] = ema(positive_sample[:,i])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(18, 4))\nplt.title('Normal powerline')\nplt.plot(flat_negative_sample, alpha=0.8);\n\nplt.figure(figsize=(18, 4))\nplt.title('Faulty powerline')\nplt.plot(flat_positive_sample, alpha=0.8);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Step 2. Identify local maxima"},{"metadata":{"trusted":true},"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]\n\n@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":{},"cell_type":"markdown","source":"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:\n\n(https://www.kaggle.com/mark4h/vsb-1st-place-solution)"},{"metadata":{"trusted":true},"cell_type":"code","source":"a = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 1, 0, 0, 3, 0, 5, 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":{},"cell_type":"markdown","source":"## Step 3. Filter the peaks to separate signal from noise"},{"metadata":{},"cell_type":"markdown","source":"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\n1. Order the peaks by their height\n2. Calculate the gradient between each consecutive pair of peaks\n3. Smooth the gradients using a convolution operation\n4. Find the noise floor using the plateau_detection function"},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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_flatten_abs : ndarray\n        An absolute flattened version of `x`.\n    \"\"\"\n    \n    x_flatten = ema_residuals(x)\n    x_flatten_abs = np.abs(x_flatten)\n    \n    peaks_indices = local_maxima_1d_window(x_flatten_abs, window)\n    heights = x_flatten_abs[peaks_indices]\n    \n    peaks_sorted_indices = np.argsort(heights)[::-1]\n    \n    peaks_indices = peaks_indices[peaks_sorted_indices]\n    heights = heights[peaks_sorted_indices]\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_indices[: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\n\n\n@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},"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 = train_df.iloc[:, sigid].values.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":{"trusted":true},"cell_type":"code","source":"sigids = [2323, 10, 4200, 4225]\n\nfor sigid in sigids:\n    d = train_df.iloc[:, sigid].values.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":{},"cell_type":"markdown","source":"## Step 4. Scale transform"},{"metadata":{},"cell_type":"markdown","source":"Since partial discharge peaks location seems to depend on the current phase (https://www.kaggle.com/c/vsb-power-line-fault-detection/discussion/77600) it is useful to compute their location not on the time, but on the \"angle\" scale"},{"metadata":{"trusted":true},"cell_type":"code","source":"@numba.jit(nopython=True, parallel=True)\ndef calculate_current_phase(data):\n    \"\"\"Calculate the current phase shift relative to sine wave.\n    Assumes the signal is 800000 data points long\n    \"\"\"\n    n = 800000\n    assert data.shape[0] == n\n    \n    # uses 50Hz Fourier coefficient\n    omegas = np.exp(-2j * np.pi * np.arange(n) / n)\n    res = np.zeros(data.shape[1], dtype=omegas.dtype)\n    for i in numba.prange(data.shape[1]):\n        res[i] = omegas.dot(data[:, i].astype(omegas.dtype))\n            \n    return np.angle(res, deg=False)\n\n\ndef to_angle(x, phase):\n    dt = 1/800000\n    return (np.degrees(2*np.pi*dt*x + phase) + 90) % 360","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sampl_id = 0\n\nx = train_df.iloc[:,sampl_id:sampl_id+1].values\nphase = calculate_current_phase(x)\nangles = to_angle(np.arange(800000), phase)\n\nplt.figure(figsize=(18, 4))\nplt.plot(x, alpha=.8)\n\nplt.figure(figsize=(18, 4))\nplt.plot(angles, x, alpha=.8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocess"},{"metadata":{"trusted":true},"cell_type":"code","source":"def process_measurement(data_df, meta_df):\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    Returns\n    -------\n    peaks : pandas.DataFrame\n        Data for each peak in the three traces in `data`.\n    \"\"\"\n    peaks = []\n    for i, sig_id in enumerate(meta_df.signal_id):\n        mat = []\n        signal = data_df.iloc[:, i].values.astype(np.float)\n        px, h = get_peaks(signal)\n        mat.append(px)\n        mat.append(h)\n        mat.append([sig_id]*len(px))\n        peaks.append(np.asarray(mat))    \n\n    peaks = pd.DataFrame(\n        np.concatenate(peaks, axis=1).T,\n        columns=['px', 'height', 'signal_id']\n    )\n\n    # Calculate the phase resolved location of each peak\n    phase_50hz = calculate_current_phase(data_df.values)\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    peaks = pd.merge(peaks, phase_50hz, on='signal_id', how='left')\n\n    peaks['phase_aligned_x'] = to_angle(peaks['px'], peaks['phase_50hz'])\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    return peaks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_peaks = process_measurement(train_df, meta_train_df)\ntrain_peaks = pd.merge(train_peaks, meta_train_df[['signal_id', 'id_measurement', 'target']], on='signal_id', how='left')\n\ndel train_df\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_peaks.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_peaks.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Features"},{"metadata":{},"cell_type":"markdown","source":"We compute basic statistics for each meaurement:\n    \n1. Total, count of peaks\n2. Count of peaks in 0, 2 and 1, 3 quarters\n3. Average height and standard deviation of peak heights in 0, 2 quarters"},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_features(peaks_df, meta_df):\n    results = pd.DataFrame(index=meta_df['id_measurement'].unique())\n    results.index.rename('id_measurement', inplace=True)\n\n\n    # Count total peaks for each measurement id\n    res = peaks_df.groupby('id_measurement').agg({\n        'px': 'count',\n    })\n    res.columns = [\"peak_count_total\"]\n    results = pd.merge(results, res, on='id_measurement', how='left')\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        'px': 'count',\n    })\n    res.columns = [\"peak_count_Q02\"]\n    results = pd.merge(results, res, on='id_measurement', how='left')\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        'px': 'count',\n    })\n    res.columns = ['peak_count_Q13']\n    results = pd.merge(results, res, on='id_measurement', how='left')\n\n\n    # Calculate height properties using 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        'height': ['mean', 'std'],\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},"cell_type":"code","source":"X_train = calculate_features(train_peaks, meta_train_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train[y_train==0].peak_count_total.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train[y_train==1].peak_count_total.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train = meta_train_df.groupby('id_measurement')['target'].sum() > 0\ny_train = y_train.astype(np.float)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train model"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_folds = 5\n\nnp.random.seed(13)\n\nsplits = np.zeros(X_train.shape[0], dtype=np.int)\nm = y_train == 1\nsplits[m] = np.random.randint(0, num_folds, size=m.sum())\nm = y_train == 0\nsplits[m] = np.random.randint(0, num_folds, size=m.sum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.Series(splits).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {\n    'objective': 'binary',\n    'metric': 'binary_logloss',\n    \n    'boosting': 'gbdt',\n    'learning_rate': 0.01,\n    'num_leaves': 80,\n    'num_boost_round': 10000,\n    \n    'feature_fraction': 0.8,\n    'bagging_fraction': 0.8,\n\n    'num_threads': 4,\n    'seed': 23974,\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"models = []\ncv_scores = []\nval_cv_scores = []\nfeature_names = X_train.columns.tolist()\n\nyp_train = np.zeros(X_train.shape[0])\nyp_val = np.zeros(X_train.shape[0])\nyp_test = np.zeros(X_train.shape[0])\n\nfor fold in range(num_folds):\n    val_fold = fold\n    test_fold = (fold + 1) % num_folds\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_train.values[train_indices],\n        y_train[train_indices],\n        feature_name=feature_names,\n    )\n    val = lgb.Dataset(\n        X_train.values[val_indices],\n        y_train[val_indices],\n        feature_name=feature_names,\n    )\n    test = lgb.Dataset(\n        X_train.values[test_indices],\n        y_train[test_indices],\n        feature_name=feature_names,\n    )\n\n    # train model\n    model = lgb.train(\n        params, \n        trn, \n        valid_sets=(trn, test, val), \n        valid_names=(\"train\", \"test\", \"validation\"), \n        early_stopping_rounds=10,\n        verbose_eval=50\n    )\n\n    # predict\n    yp = model.predict(X_train.values[train_indices])\n    yp_train[train_indices] += yp\n    yp_val_fold = model.predict(X_train.values[val_indices])\n    yp_val[val_indices] += yp_val_fold\n    yp_test_fold = model.predict(X_train.values[test_indices])\n    yp_test[test_indices] += yp_test_fold\n    \n    # save \n    models.append(model)\n    cv_scores.append(model.best_score['test']['binary_logloss'])\n    val_cv_scores.append(model.best_score['validation']['binary_logloss'])\n\nyp_train /= (num_folds - 2)\ncv_scores = np.asarray(cv_scores)\nval_cv_scores = np.asarray(val_cv_scores)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"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 Test Logloss: {:.4f} +/- {:.4f} ({:.4f})\".format(cv_scores.mean(), cv_scores.std()/np.sqrt(cv_scores.shape[0]), cv_scores.std()))\n\nprint(\"Train  accuracy: {:.4f}\".format(metrics.accuracy_score(y_train, yp_train > 0.5)))\nprint(\"CV Val accuracy: {:.4f}\".format(metrics.accuracy_score(y_train, yp_val > 0.5)))\nprint(\"CV Test accuracy: {:.4f}\".format(metrics.accuracy_score(y_train, yp_test > 0.5)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"thresholds = np.linspace(.01, .9, 90)\n\nscores_train = []\nscores_val = []\nscores_test = []\n\nfor t in thresholds:\n    s_train = metrics.f1_score(\n        y_train.values.astype(np.float), \n        yp_train > t\n    )\n    s_val = metrics.f1_score(\n        y_train.values.astype(np.float), \n        yp_val > t\n    )\n    s_test = metrics.f1_score(\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(thresholds[np.argmax(scores_val)], ls='--')\nplt.xlabel('Threshold')\nplt.ylabel('F1 Score')\nplt.show()\n\nprint(round(np.max(scores_val), 4), thresholds[np.argmax(scores_val)])\nbest_thresh = thresholds[np.argmax(scores_val)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"best_thresh = 0.4","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_problem = yp_test > best_thresh\npred_neg     = yp_test <= best_thresh\n\ntrue_problem = y_train > 0.5\ntrue_neg     = y_train <= 0.5","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(true_problem & pred_problem).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(pred_neg & true_neg).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(pred_neg & true_problem).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(pred_problem & true_neg).sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)\nimportances.groupby('feature').gain.mean().sort_values(ascending=True).plot(kind='barh');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Analysis"},{"metadata":{"trusted":true,"scrolled":false},"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    print(f)\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()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict on test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"del X_train, train_peaks\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NUM_TEST_CHUNKS = 10\n\ntest_chunk_size = int(np.ceil((meta_test_df.shape[0]/3.)/float(NUM_TEST_CHUNKS))*3.)\n\ntest_peaks = []\n\nfor 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 = pd.read_parquet(\n        data_dir + '/test.parquet',\n        columns=[str(c) for c in signal_ids]\n    )\n\n    p = process_measurement(\n        test_df, \n        meta_test_df.iloc[j_start:j_end], \n    )\n\n    test_peaks.append(p)\n\n    print(j)\n\n    del test_df\n    gc.collect()\n\n\ntest_peaks = pd.concat(test_peaks)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_peaks = pd.merge(test_peaks, meta_test_df[['signal_id', 'id_measurement']], on='signal_id', how='left')\ntest_peaks.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = calculate_features(test_peaks, meta_test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"yp_test = np.zeros(X_test.shape[0])\n\nfor j in range(len(models)):\n    model = models[j]\n    yp_test += model.predict(X_test.values)/len(models)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_submission = pd.DataFrame(\n    yp_test,\n    index=X_test.index,\n    columns=['probability']\n)\ntest_submission['target'] = (yp_test > best_thresh).astype(np.int)\n\ntest_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},"cell_type":"code","source":"test_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_submission.target.value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = test_submission[['signal_id', 'target']]\nsubmission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}