{"cells":[{"metadata":{"_uuid":"8c535b17163fb3d8c4081e4c9b0d6187362510fe"},"cell_type":"markdown","source":"# Copy from \nhttps://www.kaggle.com/junkoda/handmade-features\nI used my psuedo label to improve the scores of this model.  Both lb score and private score got improvement. \n\nAn example of handmade features + CatBoostClassifier, for those who wants to try something other than RNN."},{"metadata":{"trusted":true,"_uuid":"ecbcc82c6dc8d80cfde32c918168f6bc74394456"},"cell_type":"code","source":"import os\nos.listdir('../input')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":true},"cell_type":"code","source":"# Load training data\nimport numpy as np\nimport pandas as pd\nimport pyarrow.parquet as pq\nfrom tqdm import tqdm_notebook as tqdm\n\nv_raw_train = pq.read_pandas('../input/vsb-power-line-fault-detection/train.parquet').to_pandas().values\nmeta_train = np.loadtxt('../input/vsb-power-line-fault-detection/metadata_train.csv', skiprows=1, delimiter=',')\ny_train = meta_train[:, 3].astype(bool)\n\nprint(v_raw_train.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b6a16e5c566ea82c0b276e069e8e9f22269d71cc"},"cell_type":"markdown","source":"# Feature engineering\n\n## Compute *spectra*\n\nI first compute *spectra*, `percentiles - mean` for every `m=1000` chunk, a common data compression technique in this competition."},{"metadata":{"trusted":true,"_uuid":"5fac122ee6be83c8876aa9f1146d4399cdbacbd4"},"cell_type":"code","source":"def compute_spectra(v_raw, *, m = 1000):\n    \"\"\"\n    compute mean and percentile - mean for every chunk of m data\n    \n    Args:\n      v_raw (array): 800,000 x n_sample; the input\n      m (int): the chunk size\n    \n    Returns: d (dict)\n      d['mean']: mean in each chunks\n      d['percentile']: percentile - mean\n    \"\"\"\n    percentile = (100, 99, 95, 0, 1, 5)\n    \n    n = v_raw.shape[1] # number of samples\n    length = v_raw.shape[0] // m # 800,000 -> 800\n    n_spectra = len(percentile)\n    \n    mean_signal = np.zeros((n, length), dtype='float32') # mean in each chunk\n    percentile_spectra = np.zeros((n, length, n_spectra), dtype='float32')\n    \n    # compute spectra\n    print('computing spectra...', flush=True)\n    for i in tqdm(range(n)):\n        v = v_raw[:, i].astype('float32').reshape(-1, m) / 128.0\n        \n        mean = np.mean(v, axis=1)        \n        s = np.abs(np.percentile(v, percentile, axis=1) - mean)\n        \n        # subtract baseline\n        h = np.percentile(s, 5.0)\n        s = np.maximum(0.0, s - h)\n\n        mean_signal[i, :] = mean\n        percentile_spectra[i, :, :] = s.T\n            \n    d = {}\n    d['mean'] = mean_signal\n    d['percentile'] = percentile_spectra\n    \n    return d\n\nspec_train = compute_spectra(v_raw_train)\nprint('done.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3e4242fbd37a835715f96f9ee3a69b5097a850ed"},"cell_type":"code","source":"%matplotlib inline\nimport matplotlib.pyplot as plt\n\ni_vis = 0\n\nplt.title('mean signal')\nplt.ylabel('mean')\n        \nfor iphase in range(3):\n    plt.plot(spec_train['mean'][i_vis + iphase, :], label='phase %d' % iphase)\n\nplt.legend(loc=1)\nplt.show()\n\n#\n# spectra\n#\nnrow = 2\nncol = 3\nplt.figure(figsize=(3*ncol, 3*nrow))\nplt.suptitle('Spectra -- deviation from the mean')\n\nfor icol in range(3):\n    plt.subplot(nrow, ncol, icol + 1)\n    plt.title('phase %d' % icol)\n    plt.ylim(0, 0.3)\n    \n    if icol == 0:\n        plt.ylabel('max - mean')\n        \n    # max\n    plt.plot(spec_train['percentile'][i_vis + icol, :, 0])\n    \n    # 99% percentile\n    plt.subplot(nrow, ncol, icol + ncol + 1)\n    plt.ylim(0, 0.025)\n    \n    if icol == 0:\n        plt.ylabel('99% - mean')\n    plt.plot(spec_train['percentile'][i_vis + icol, :, 1])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2c7fde7001b21227800f0823ceb5a6374ae636fb"},"cell_type":"markdown","source":"## Peak interval\n\nThen, within the 800 chunks of the spectra, I find the *peak interval* of `width=150` that contains the maximum deviation in the max - mean spectrum."},{"metadata":{"trusted":true,"_uuid":"0529f8186f830aba9227dcaaf063e174d2358765"},"cell_type":"code","source":"import tensorflow as tf\n\ndef max_windowed(spec, *, width=150, stride=10):\n    \"\"\"\n    Smooth the spectrum with a tophat window function and find the\n    peak inteval that maximises the smoothed spectrum.\n    \n    Returns: d(dict)\n      d['w'] (array): smoothed max - mean spectrum\n      d['ibegin'] (array): the left edge index of the peak interval\n    \"\"\"\n    n = spec.shape[0]\n    length = spec.shape[1] # 800\n    nspec = spec.shape[2] # 6 spectra\n\n    n_triplet = n // 3\n\n    # Reorganize the max spectrum from 8712 data to 2904 triplets with 3 phases\n    max_spec3 = np.empty((n_triplet, length, 3))\n    for i_triplet in range(n_triplet):\n        max_spec3[i_triplet, :, 0] = spec[3*i_triplet, :, 0] # phase 0\n        max_spec3[i_triplet, :, 1] = spec[3*i_triplet + 1, :, 0] # phase 1\n        max_spec3[i_triplet, :, 2] = spec[3*i_triplet + 2, :, 0] # phase 2\n\n    x = tf.placeholder(tf.float32, [None, length, 3]) # input spectra before smoothing\n    # 800 -> 80: static convolaution\n    # convolution but not CNN, the kernel is static\n    # smoothing/convolution kernel\n    # tophat window function\n    # shape (3, 1) adds up 3 phases to one output\n    K = np.ones((width, 3, 1), dtype='float32') / width\n\n    W_conv1 = tf.constant(K)\n    h_conv1 = tf.nn.conv1d(x, W_conv1, stride=stride, padding='VALID')\n    \n    with tf.Session() as sess:\n        w = sess.run(h_conv1, feed_dict={x:max_spec3})\n\n    imax = np.argmax(w[:, :, 0], axis=1) # index of maximum smoothed spectrum\n    \n    d = {}\n    d['w'] = w # smoothed max spectrum\n    d['ibegin'] = imax*stride\n    \n    return d\n\npeaks = max_windowed(spec_train['percentile'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"4c44c66cb047b1cf5cd496a15020d577b53582c6"},"cell_type":"code","source":"i_vis = 0\nnrow = 1\nncol = 3\nplt.figure(figsize=(3*ncol, 3*nrow))\nplt.suptitle('spectra')\n\nfor j in range(3):\n    plt.subplot(nrow, ncol, j + 1)\n    plt.ylim(0, 0.3)\n    \n    ibegin = peaks['ibegin'][i_vis]\n    plt.axvspan(ibegin, ibegin + 150, alpha=0.2, color='red')\n    \n    # max - mean spectrum\n    plt.plot(spec_train['percentile'][i_vis + j, :, 0])\n    \n    # smoothed max - mean spectrum\n    xx = np.arange(peaks['w'].shape[1])*10 + 150//2\n    plt.plot(xx, peaks['w'][i_vis])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"48184f695d6b61dc7d2170f312e38f9f97af0fe2"},"cell_type":"markdown","source":"#### Fig 3: Peak interval\n* The *orange line* is the smooted `max - mean` spectrum (sum of three phases)\n* *Red band* is the `peak interval`, the interval that gives the maximum smoothed line (orange)[](http://)\n\nNow, I compute features in the peak interval.\n\nMean and Max of `percentile - mean` spectra in the interval."},{"metadata":{"trusted":true,"_uuid":"18e11978d856cd6fd05b528dd8a5cceb13465619"},"cell_type":"code","source":"def compute_features(v_raw, spec=None):\n    \"\"\"\n    Args:\n      v_raw (array): The original 800,000 x 8712 training data\n      spec (dict): The result of compute_spectra() if already computed.\n                   If it is None, it will be computed automatically.\n    \n    Returns:\n       X (array): Feature vector of shape (2904, 57)\n                  2904 triplets, 57 features\n    \"\"\"\n    if spec is None:\n        spec = compute_spectra(v_raw)\n    \n    v_spec = spec['percentile']\n    shape = v_spec.shape\n    n = shape[0] # number of data\n    length = shape[1]\n    nspec = shape[2]\n    \n    n_triplet = n // 3\n    \n    # Reorder to i_triplet, phase\n    spec3 = np.empty((n_triplet, length, nspec, 3))\n    \n    for i_triplet in range(n_triplet):\n        spec3[i_triplet, :, :, 0] = v_spec[3*i_triplet, :, :] # phase 0\n        spec3[i_triplet, :, :, 1] = v_spec[3*i_triplet + 1, :, :] # phase 1\n        spec3[i_triplet, :, :, 2] = v_spec[3*i_triplet + 2, :, :] # phase 2\n\n    # extract \"max-windowed\" from the spectra\n    width = 150\n    peaks = max_windowed(v_spec, width=width)\n    \n    # Feature vector\n    n_feature4 = 3\n    X = np.empty((n_triplet, n_feature4*nspec*3 + 3))\n    \n    # features for each percentile and phase\n    X4 = np.empty((n_triplet, n_feature4, nspec, 3)) # triplet, figure, spec type, phase\n        \n    for i_triplet in range(n_triplet):       \n        # Maximum of the spectra in the full range\n        # 18 features (6 percentiles x 3 phases)\n        X4[i_triplet, 0, :, :] = np.max(spec3[i_triplet, :, :, :], axis=0)\n        \n        # Peak interval\n        ibegin = peaks['ibegin'][i_triplet]\n        iend = ibegin + width\n        imid = ibegin + width // 2\n    \n        # Mean of the spectra in the peak inteval 18 features\n        X4[i_triplet, 1, :, :] = np.mean(spec3[i_triplet, ibegin:iend, :, :], axis=0)\n        \n        # Max of the spectra in the peak inteval (18 features)\n        X4[i_triplet, 2, :, :] = np.max(spec3[i_triplet, ibegin:iend, :, :], axis=0)\n        \n        # Mean signal at the midpoint of the interval (3 features)\n        X[i_triplet, 0] = spec['mean'][3*i_triplet,     imid]\n        X[i_triplet, 1] = spec['mean'][3*i_triplet + 1, imid]\n        X[i_triplet, 2] = spec['mean'][3*i_triplet + 2, imid]\n    \n    shape = X4.shape\n    \n    # Flatten the X4 tensor\n    # 3 + 18x3 = 57 features\n    X[:, 3:] = X4.reshape(shape[0], shape[1]*shape[2]*shape[3])\n    \n    return X\n\nX_all3 = compute_features(v_raw_train, spec_train)\n\n# The label for the triple\n# True iff two or more labels in 3 phases are True\ny_all3 = np.sum(y_train.reshape(-1, 3), axis=1) >= 2\n\nprint('Three phases are combined into one training data; the shapes are, therefore,')\nprint(X_all3.shape, y_all3.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3c555fa8beb77c649587efbc977cf8fc9276f425"},"cell_type":"code","source":"plt.figure()\nplt.title('mean max vs 99 percentile spectra in the peak interval')\nplt.xlabel('max')\nplt.ylabel('99 percentile')\nplt.ylim(0, 0.11)\nplt.xlim(0, 0.2)\n\n\nj_max = 3 + 18 + 0\nj_99 = 3 + 18 + 3\n\nidx_pos = y_all3 # positive samples\nidx_neg = np.logical_not(y_all3) # negative samples\n\nplt.plot(X_all3[idx_neg, j_max], X_all3[idx_neg, j_99], ',', color='gray', label='False')\nplt.plot(X_all3[idx_pos, j_max], X_all3[idx_pos, j_99], ',', color='red', label='True')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1fa05a3391c318f9ff541855f2d95204b25a456c"},"cell_type":"markdown","source":"#### Fig. 4\n\nA plot of 2 features:\n\n* x axis: mean \"max - mean spectra\" in the peak interval\n* y asix: mean \"99% percentile - mean\" spectra in the peak intervbal\n\nPositive samples tends to have small `99% - mean` feature for given `max - mean` feature,\n"},{"metadata":{"_uuid":"655a76d8d6890013a058eab66ddb7595b3f5ce3b"},"cell_type":"markdown","source":"# Training and Cross Validation\n\n1. Split to training and CV data width 5-fold cross validation"},{"metadata":{"trusted":true,"_uuid":"1a799f21c8d67b610ade53efbebbf03a99dce2e0"},"cell_type":"code","source":"#Loading the psuedo labels\n# knowledge_data = pd.read_csv('../input/vsb-knowledge-0744/sub_ens_v19.csv')\nknowledge_data = pd.read_csv('../input/vsb-knowledge-0744/sub_ens_v16.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"541704d0a6ab46fdcb5ab1368957c94099620678"},"cell_type":"code","source":"# Release RAM of the training data \nif 'v_raw_train' in globals():\n    del v_raw_train\n\n# Load test data\nid_test = np.loadtxt('../input/vsb-power-line-fault-detection/metadata_test.csv', skiprows=1, delimiter=',')[:, 0].astype(int)\nn_test = len(id_test)\n\nX_tests = []\n\n# Load test data and compute the feature vector\n# The test data is split into 4 to fit it into RAM\nn_subset = 4\nnread = 0\n\nfor i_subset in range(n_subset):\n    # signal_id range in the test data; 8712 is the first data in the test.parquet\n    ibegin = 8712 + 3*int(n_test // 3 * (i_subset/n_subset))\n    iend = 8712 + 3*int(n_test // 3 * ((i_subset + 1)/n_subset))\n    \n    print('Loading %d/%d; signal_id %d - %d...' % (i_subset, n_subset, ibegin, iend))\n    v_raw_test = pq.read_pandas('../input/vsb-power-line-fault-detection/test.parquet',\n                                columns=[str(i) for i in range(ibegin, iend)]).to_pandas().values\n    \n    nread += v_raw_test.shape[1]\n    X = compute_features(v_raw_test)\n    X_tests.append(X)\n    print('%d/%d test data processed.' % (nread, n_test))\n\n    del v_raw_test\n\nX_test = np.concatenate(X_tests, axis=0)\nassert(X_test.shape[0] == id_test.shape[0] // 3)\n\ndel X_tests\n\nprint('X_test computation done. shape', X_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6948832f6017a3f5fb06f5582fa85a2f6cad1bf9"},"cell_type":"code","source":"len_train=len(X_all3)\n\ny_list=knowledge_data['target'].values\ny_test=[]\nfor j in range(0,len(y_list),3):\n    y_test.append(y_list[j])\ny_test=np.asarray(y_test)\ndel knowledge_data\n# gc.collect()\nprint(X_all3.shape,y_all3.shape)\nprint(X_test.shape,y_test.shape)\n#\nX_all3=np.concatenate([X_all3,X_test])\ny_all3=np.concatenate([y_all3,y_test]) \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bf27b645f31177307796cc0800a902fa7b9b0d98"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e822a190a168b22e12ff76fc6deb875940e60f40"},"cell_type":"code","source":"import sklearn.model_selection\nimport sklearn.metrics\nfrom catboost import CatBoostClassifier\n\n# split into train and CV data\nn_splits = 5\n\nmodels = []\nscores = np.zeros(n_splits)\n\nprint('Training...')\nprint('MCC training & cv') # MCC = Matthews Correlation Coefficient\n\n\nseeds=[0,42,1204,2019]\nfor seed in seeds:\n    splits = list(sklearn.model_selection.StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=seed).split(X_all3, y_all3))\n\n    for i, (idx_train, idx_cv) in enumerate(splits):\n        X_train = X_all3[idx_train, :]\n        y_train = y_all3[idx_train]\n\n        X_cv = X_all3[idx_cv, :]\n        y_cv = y_all3[idx_cv]\n\n        # Learning rate is important; large values overfit the data\n        learning_rate = 0.006\n        model = CatBoostClassifier(learning_rate=learning_rate, od_type='IncToDec',\n                                loss_function='Logloss', use_best_model=True, eval_metric='MCC')\n\n        model.fit(X_train, y_train.astype(float), silent=True,\n                  eval_set=(X_cv, y_cv.astype(float)))\n\n        y_predict_train = model.predict(X_train)\n        y_predict_cv = model.predict(X_cv)\n\n    #     score_train = sklearn.metrics.matthews_corrcoef(y_train, y_predict_train)\n        score_cv = sklearn.metrics.matthews_corrcoef(y_cv, y_predict_cv)\n\n        models.append(model)\n        scores[i] = score_cv\n    \n#     print('%d %.3f %.3f' % (i, score_train, score_cv))\n\nprint('CV scores %.3f ± %.3f' % (np.mean(scores), np.std(scores)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b36e9105b6ce14309fd41714e568d261f71a56ee"},"cell_type":"markdown","source":"# Submit predictions\n\n# Load the test data and compute the feature vector\n\n```\nv_raw_test = pq.read_pandas('../input/test.parquet').to_pandas().values\nX_test = compute_features(v_raw_test)\n```\n\nTwo lines, but extra work is necessary to put the test data on RAM of the Kernel 17.2GB."},{"metadata":{"trusted":true,"_uuid":"ef90fe85eb1f560d0fdd4ecd7f3bf55abf5c2a98","_kg_hide-output":true},"cell_type":"code","source":"# # Release RAM of the training data \n# if 'v_raw_train' in globals():\n#     del v_raw_train\n\n# # Load test data\n# id_test = np.loadtxt('../input/vsb-power-line-fault-detection/metadata_test.csv', skiprows=1, delimiter=',')[:, 0].astype(int)\n# n_test = len(id_test)\n\n# X_tests = []\n\n# # Load test data and compute the feature vector\n# # The test data is split into 4 to fit it into RAM\n# n_subset = 4\n# nread = 0\n\n# for i_subset in range(n_subset):\n#     # signal_id range in the test data; 8712 is the first data in the test.parquet\n#     ibegin = 8712 + 3*int(n_test // 3 * (i_subset/n_subset))\n#     iend = 8712 + 3*int(n_test // 3 * ((i_subset + 1)/n_subset))\n    \n#     print('Loading %d/%d; signal_id %d - %d...' % (i_subset, n_subset, ibegin, iend))\n#     v_raw_test = pq.read_pandas('../input/vsb-power-line-fault-detection/test.parquet',\n#                                 columns=[str(i) for i in range(ibegin, iend)]).to_pandas().values\n    \n#     nread += v_raw_test.shape[1]\n#     X = compute_features(v_raw_test)\n#     X_tests.append(X)\n#     print('%d/%d test data processed.' % (nread, n_test))\n\n#     del v_raw_test\n\n# X_test = np.concatenate(X_tests, axis=0)\n# assert(X_test.shape[0] == id_test.shape[0] // 3)\n\n# del X_tests\n\n# print('X_test computation done. shape', X_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c84edf65688ef6df7e965aebb1fa2ca66afd3506"},"cell_type":"code","source":"# Take average of 5 predictions\ny_test_probas = np.empty((X_test.shape[0], n_splits*len(seeds)))\n\n# assert(len(models) == n_splits)\nfor i, model in enumerate(models):\n    y_test_probas[:, i] = model.predict_proba(X_test)[:, 1]\n\ny_test_proba = np.mean(y_test_probas, axis=1)\n\n# Convert to 0 1 with a threshold 0.25, then replicate 3 copies for 3 phases\ny_submit = np.repeat(y_test_proba > 0.25, 3)\n\nprint('Positive fraction %d/%d = %.3f' % (\n    np.sum(y_submit), len(y_submit), np.sum(y_submit)/len(y_submit)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a819ef7e2f88dd34e757d4209ac07f750369d92f"},"cell_type":"code","source":"# Write predictions\nassert(len(id_test) == len(y_submit))\nn = len(id_test)\n\nwith open('submission.csv', 'w') as f:\n    f.write('signal_id,target\\n')\n    for i in range(n):\n        f.write('%d,%d\\n' % (id_test[i], int(y_submit[i])))\n\nprint('submission.csv written')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b26c9035e2a3df8a1a8f222bb52895549fcb112f"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}