{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-07T10:46:19.913867Z","iopub.execute_input":"2022-01-07T10:46:19.914348Z","iopub.status.idle":"2022-01-07T10:46:19.921096Z","shell.execute_reply.started":"2022-01-07T10:46:19.914262Z","shell.execute_reply":"2022-01-07T10:46:19.920343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mục đích của Notebook.\n\nNotebook này sẽ thực hiện nhiệm vụ biến đổi dữ liệu đầu vào thành các đặc trưng khác nhau. Sau đó thực hiện tính feature importance để giảm số lượng feature phải lấy.\n\n## Định nghĩa các số liệu đầu vào","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Thư viện để trực quan hóa tiến độ.\nfrom tqdm import tqdm_notebook as tqdm\n\n# Tổng số dòng dữ liệu\ntotal = 629145481 \n\n# Khi sử dụng chunksize càng lớn, thì dữ liệu nạp càng nhanh.\n# 2700000 là số lớn nhất chia hết cho cả 18 và 150 000. Mà 18 lại là ước của 4194 là số lần load. \n# Do đó, 2700000 được chọn làm chunksize.\nchunksize = 2700000\n\n# Tải dữ liệu sử dụng chunk, do nếu tải hết dữ liệu vào thì phần Editing sẽ báo lỗi.\nchunks = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64}, chunksize = chunksize)\n\nc = 1\n\nsegments = total // chunksize\n# for debug\n# segments = 1\n    \n# Segment là số lượng đoạn được lấy từ input.\n# Các chunk được tải sẽ được nối lại bằng làm append.\n\ntrain = None\n\nfor i in tqdm(range(segments)):\n    chunk = next(chunks)\n    \n    if train is None:\n        train = chunk\n    else:\n        train = pd.concat([train, chunk])\n    del chunk","metadata":{"execution":{"iopub.status.busy":"2022-01-07T10:46:20.597965Z","iopub.execute_input":"2022-01-07T10:46:20.598637Z","iopub.status.idle":"2022-01-07T11:01:39.895699Z","shell.execute_reply.started":"2022-01-07T10:46:20.598581Z","shell.execute_reply":"2022-01-07T11:01:39.89456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Hiển thị đến độ chính xác 15 chữ số sau dấu phẩy.\npd.options.display.precision = 15","metadata":{"execution":{"iopub.status.busy":"2022-01-07T11:01:39.899558Z","iopub.execute_input":"2022-01-07T11:01:39.899821Z","iopub.status.idle":"2022-01-07T11:01:39.908547Z","shell.execute_reply.started":"2022-01-07T11:01:39.899776Z","shell.execute_reply":"2022-01-07T11:01:39.907294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# In ra 5 dòng đầu tiên của acoustic data\ntrain[:5]","metadata":{"execution":{"iopub.status.busy":"2022-01-07T11:01:39.910119Z","iopub.execute_input":"2022-01-07T11:01:39.910641Z","iopub.status.idle":"2022-01-07T11:01:39.952094Z","shell.execute_reply.started":"2022-01-07T11:01:39.910574Z","shell.execute_reply":"2022-01-07T11:01:39.951404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature engineering","metadata":{}},{"cell_type":"markdown","source":"## Chuẩn bị các thư viện sẽ sử dụng.","metadata":{}},{"cell_type":"code","source":"%%time\n# Load các thư viện liên quan đến biến đổi tín hiệu và thống kê sẽ được sủ dụng.\nfrom scipy import stats\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\n\n# Load thư viện về hồi quy tuyến tính đơn giản để tìm xu hướng của dữ liệu.\nfrom sklearn.linear_model import LinearRegression\n\nimport numpy as np","metadata":{"execution":{"iopub.status.busy":"2022-01-07T11:01:39.953296Z","iopub.execute_input":"2022-01-07T11:01:39.953736Z","iopub.status.idle":"2022-01-07T11:01:39.96226Z","shell.execute_reply.started":"2022-01-07T11:01:39.953678Z","shell.execute_reply":"2022-01-07T11:01:39.961315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n# Lọc bỏ các warning để tập trung hơn vào tiến độ xử lý dữ liệu.\nimport warnings\nwarnings.filterwarnings('ignore')\n\ndef classic_sta_lta(x, length_sta, length_lta):\n    \"\"\"\n    Input: x là một mảng số.\n    STA: Short Term Avg.\n    LTA: Long Term Avg.\n    \n    - Term Average: được tính bằng thương độ khác nhau của các phần tử cách nhau một khoảng cố định trên khoảng cố định đó.\n    - Short Term Avg: Term Average với khoảng nhỏ.\n    - Long Term Avg: Term Average với khoảng lớn.\n    \"\"\"\n    # Tổng cộng dồn các phần tử bình phương.\n    sta = np.cumsum(x ** 2)\n\n    # Chuyển đổi sang kiểu float.\n    sta = np.require(sta, dtype=np.float)\n\n    # Tương tự với mảng lta\n    lta = sta.copy()\n\n    # Tính STA là LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n\n    # Tiếp theo ta tiến hành chia sta cho lta. Các giá trị nan sẽ được biến thành 0.\n    sta[:length_lta - 1] = 0\n\n    # Định nghĩa dtiny là số nhỏ nhất mà máy tính có thể chứa.\n    dtiny = np.finfo(0.0).tiny\n    \n    # Với mỗi lta nhỏ hơn dtiny, sẽ đặt bằng dtiny.\n    idx = lta < dtiny\n    lta[idx] = dtiny\n    \n    return sta / lta\n    \n\ndef get_feature(ad):\n    \"\"\"\n    Đầu vào: ad (Acoustic Data): pd.Series: chứa 150 000 dòng dữ liệu.\n    Đầu ra:\n    - Các thông số thống kê của ad. Được tính cho cả ad và trị tuyệt đối của ad.\n    \"\"\"\n    features = dict()\n    ## Trung bình, đồ lệch chuẩn, max, min và khoảng cách giữa max và min của ad.\n    features['ave'] = ad.mean()\n    features['std'] = ad.std() \n    features['max'] = ad.max()\n    features['min'] = ad.min() \n    features['min_max_diff'] = features['max'] - features['min']\n\n    features[ 'abs_ave'] = np.abs(ad).mean()\n    features['abs_std'] = np.abs(ad).std()\n    features['abs_max'] = np.abs(ad).max() \n    features['abs_min'] = np.abs(ad).min()\n    features['abs_min_max_diff'] = features['abs_max'] - features['abs_min']\n\n    ## Số lượng thời điểm ad có giá trị > 400.\n    features['count_abnormal'] = len(ad[np.abs(ad) > 400])\n    ## Tổng các ad.\n    features['sum'] = ad.sum()\n\n    # Tìm xu hướng của ad sử dụng hồi quy tuyến tính.\n    # Thuật toán hồi quy tuyến tính sẽ được nói rõ trong notebook về thực nghiệm.\n    model = LinearRegression().fit(np.arange(len(ad)).reshape(-1, 1), ad)\n\n    features['coef_'] = model.coef_[0]\n    \n    model = LinearRegression().fit(np.arange(len(ad)).reshape(-1, 1), np.abs(ad))\n\n    features['abs_coef_'] = model.coef_[0]\n    \n    # Chia ad thành các phần nhỏ hơn và tính mean, std, max, min và khoảng cách giữa max và min.\n    div = 5\n    for i in range(div):\n        features['ave_div_' + str(i)] = ad[rows // div * i: rows // div * (i + 1)].mean()\n        features['std_div_' + str(i)] = ad[rows // div * i: rows // div * (i + 1)].std()\n        features['max_div_' + str(i)] = ad[rows // div * i: rows // div * (i + 1)].max()\n        features['min_div_' + str(i)] = ad[rows // div * i: rows // div * (i + 1)].min()\n        features['min_max_diff_div_' + str(i)] = features['max_div_' + str(i)] - features['min_div_' + str(i)]\n\n        features['abs_ave_div_' + str(i)] = np.abs(ad[rows // div * i: rows // div * (i + 1)]).mean()\n        features['abs_std_div_' + str(i)] = np.abs(ad[rows // div * i: rows // div * (i + 1)]).std()\n        features['abs_max_div_' + str(i)] = np.abs(ad[rows // div * i: rows // div * (i + 1)]).max() \n        features['abs_min_div_' + str(i)] = np.abs(ad[rows // div * i: rows // div * (i + 1)]).min()\n        features['abs_min_max_diff_div_' + str(i)] = features['abs_max_div_' + str(i)] - features['abs_min_div_' + str(i)]  \n\n    # trung bình điều hòa. x_i phải dương.\n    # H = n / (1/x_1 + 1/x_2 + ... + 1/x_n)\n    features['hmean'] = stats.hmean(np.abs(ad[ad != 0]))\n    # Trung bình nhân. x_i phải dương.\n    # G = power(x_1 * x_2 * ... * x_n, 1/n)\n    features['gmean'] = stats.gmean(np.abs(ad[ad != 0]))\n    \n    # Quantile: lượng tử. Liên kết với một tỉ lệ k.\n    # Là điểm mà các vị trí dưới nó có tổng tỉ lệ xuất hiện với xác suất k.\n    \n    # Percentile là quantile sử dụng các tỉ lệ được chia theo các số nguyên của 100%.\n    percentages = [0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 0.99]\n    for percentage in percentages:\n        features['percentile_' + str(percentage)] = np.percentile(ad, percentage)\n        \n    for quantile in percentages:\n        features['quantile_' + str(quantile)] = np.quantile(ad, quantile)\n        \n    for percentage in percentages:\n        features['abs_percentile_' + str(percentage)] = np.percentile(np.abs(ad), percentage)\n        \n    for quantile in percentages:\n        features['abs_quantile_' + str(quantile)] = np.quantile(np.abs(ad), quantile)\n        \n    # Median absolute deviation: Trung bình trị tuyệt đối của độ lệch.\n    features['mad'] = np.mean(np.abs(ad - np.mean(ad)))\n    # Độ lệch:\n    features['skew'] = stats.skew(ad)\n    # Độ nhọn\n    features['kurtosis'] = stats.kurtosis(ad)\n    # Trung vị,\n    features['median'] = np.median(ad)\n    \n    features['abs_mad'] = np.mean(np.abs(np.abs(ad) - np.mean(np.abs(ad))))\n    features['abs_skew'] = stats.skew(np.abs(ad))\n    features['abs_kurtosis'] = stats.kurtosis(np.abs(ad))\n    features['abs_median'] = np.median(np.abs(ad))\n    \n    features['hilbert_mean'] = np.abs(hilbert(ad)).mean()\n    \n    # Phép làm min sử dụng vector Hann.\n    # Tính tích chập giữa ad và vector hann, sau đó chia cho tổng của vector hann và tính mean.\n    windows = [10, 50, 100, 200, 500, 1000, 5000, 10000]\n    for window in windows:\n        features['hann_window_' + str(window)] = (convolve(ad, hann(window), mode='same') / sum(hann(window))).mean()\n    \n    # Tính các đặc tính theo tổng liên tiếp (rolling)\n    for window in windows:\n        # pd.Series.rolling nếu dùng hàm mặc định sẽ trả về tổng các phần tử liên tiếp trong một khoảng window.\n        roll_mean = ad.rolling(window).mean().dropna()\n        roll_std = ad.rolling(window).std().dropna()\n        \n        # Tính mean, min, max, std, percentile trên trên tập roll_mean và roll_std.\n        \n        features['rolling_max_ave_' + str(window)] = roll_mean.max()\n        features['rolling_min_ave_' + str(window)] = roll_mean.min()\n        \n        features['rolling_max_std_' + str(window)] = roll_std.max()\n        features['rolling_min_std_' + str(window)] = roll_std.min()\n        \n        features['abs_rolling_max_std_' + str(window)] = np.abs(roll_std).max()\n        features['abs_rolling_min_std_' + str(window)] = np.abs(roll_std).min()\n        \n        features['abs_rolling_max_mean_' + str(window)] = np.abs(roll_mean).max()\n        features['abs_rolling_min_mean_' + str(window)] = np.abs(roll_mean).min()\n        \n        for percentage in percentages:\n            features['abs_rolling_mean_{}_percentage_{}'.format(window, percentage)] = np.percentile(roll_mean, percentage)\n            features['abs_rolling_std_{}_percentage_{}'.format(window, percentage)] = np.percentile(roll_std, percentage)\n            \n    # tính các đặc tính theo hiệu liên tiếp\n    dif = ad.diff().dropna().values\n    \n    features['dif_mean'] = dif.mean()\n    \n    for percentage in percentages:\n        features['dif_percentage_{}'.format(percentage)] = np.percentile(dif, percentage)\n        features['dif_percentage_{}'.format(percentage)] = np.percentile(dif, percentage)\n        \n    # Tính trunh bình mức độ thay đổi giữa các giá trị liên tiếp trong ad.\n    change = (dif / ad.values[:-1])\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    \n    features['change_mean'] = change.mean()\n    # sử dụng short term average và long term average\n    \n    features['classic_sta_lta_mean_1'] = classic_sta_lta(ad, 50, 10000).mean()\n    features['classic_sta_lta_mean_1'] = classic_sta_lta(ad, 200, 10000).mean()\n    features['classic_sta_lta_mean_1'] = classic_sta_lta(ad, 500, 10000).mean()\n    features['classic_sta_lta_mean_1'] = classic_sta_lta(ad, 1000, 10000).mean()\n    \n    # sử dụng exponential Moving Average\n    # Là phương pháp làm mịn phổ biến của xử lý tín hiệu\n    for window in windows:\n        features['EMA_{}'.format(window)] = ad.ewm(span=window).mean(skipna=True).mean(skipna=True)\n    \n    # Sử dụng Interquartile range\n    # Là khoảng cách giữa percentile .75 và percentile .25.\n    features['iqr'] = np.subtract(*np.percentile(ad, [75, 25]))\n    \n    # xóa biến để giải phóng bộ nhớ\n    del change\n    del dif\n    \n    return features\n\n\nrows = 150000\nsegments = len(train) // rows\n\nX = []\ny = []\n\nfor segment in tqdm(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n\n    ad = seg['acoustic_data']\n    ttf = seg['time_to_failure'].values[-1]\n    \n    y.append(ttf)\n    \n    ## Tính với dữ liệu chính, dữ liệu thực và phức của fft.\n    \n    ft = np.fft.fft(ad)\n    real_ft = pd.Series(np.real(ft))\n    imag_ft = pd.Series(np.imag(ft))\n    \n    # Nhận một dict() các đặc trưng.\n    normal_feature = get_feature(ad)\n    real_feature = get_feature(real_ft)\n    img_feature = get_feature(imag_ft)\n    \n    # Nối các dict() lại với nhau.\n    for key, value in real_feature.items():\n        normal_feature['real_' + key] = value\n    for key, value in img_feature.items():\n        normal_feature['imag_' + key] = value\n    \n    X.append(normal_feature)","metadata":{"execution":{"iopub.status.busy":"2022-01-07T11:24:39.128277Z","iopub.execute_input":"2022-01-07T11:24:39.132461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = pd.DataFrame(X)\ny = pd.Series(y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Chạy mô hình catboost cơ bản để tính feature importance","metadata":{}},{"cell_type":"markdown","source":"Tiền xử lý dữ liệu sử dụng StandardScaler.\n\nStandardScaler sẽ đưa dữ liệu về mean 0 và std 1.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import train_test_split\n\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X)\n\nX_train, X_valid, y_train, y_valid = train_test_split(X_scaled, y, train_size=0.7)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T10:16:07.930324Z","iopub.execute_input":"2022-01-04T10:16:07.930828Z","iopub.status.idle":"2022-01-04T10:16:07.942629Z","shell.execute_reply.started":"2022-01-04T10:16:07.930757Z","shell.execute_reply":"2022-01-04T10:16:07.941541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Định nghĩa một mô hình catboost cơ bản.","metadata":{}},{"cell_type":"code","source":"%%time\nimport catboost as cb\n\ntrain_data = cb.Pool(X_train, y_train)\nvalid_data = cb.Pool(X_valid, y_valid)\n\nmodel = cb.CatBoostRegressor(\n    loss_function='MAE', \n    logging_level=\"Silent\",\n    iterations=300,\n    learning_rate=0.1,\n    depth=8,\n    l2_leaf_reg=0.2,\n    random_seed=500\n)\nmodel.fit(\n    train_data,             \n    verbose=False\n)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T10:16:07.944072Z","iopub.execute_input":"2022-01-04T10:16:07.944592Z","iopub.status.idle":"2022-01-04T10:16:11.010567Z","shell.execute_reply.started":"2022-01-04T10:16:07.944541Z","shell.execute_reply":"2022-01-04T10:16:11.009566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Tính và trực quan hóa feature importance.\n\nGiữ lại 120 feature có độ quan trọng cao nhất.","metadata":{}},{"cell_type":"code","source":"feature_importance = model.get_feature_importance(\n    train_data\n)\n\nfeature_importance = pd.DataFrame({\n    \"feature_name\": X.columns,\n    \"score\": feature_importance\n})\n\nfeature_importance.sort_values([\"score\"], axis=0, ascending=False, inplace=True)\n\nkeep_feature = feature_importance[feature_importance['score'] > 0.00001]\nkeep_feature = keep_feature.iloc[:120]\nkeep_feature","metadata":{"execution":{"iopub.status.busy":"2022-01-04T10:16:11.011837Z","iopub.execute_input":"2022-01-04T10:16:11.012063Z","iopub.status.idle":"2022-01-04T10:16:11.052325Z","shell.execute_reply.started":"2022-01-04T10:16:11.012027Z","shell.execute_reply":"2022-01-04T10:16:11.051374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for feature_name in keep_feature['feature_name']:\n    print(feature_name)\n    \nkeep_feature['feature_name'].values","metadata":{"execution":{"iopub.status.busy":"2022-01-04T10:16:11.053629Z","iopub.execute_input":"2022-01-04T10:16:11.053909Z","iopub.status.idle":"2022-01-04T10:16:11.066755Z","shell.execute_reply.started":"2022-01-04T10:16:11.053864Z","shell.execute_reply":"2022-01-04T10:16:11.065114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Các thuộc tính này sẽ được giữ lại để làm đặc trưng. Cho các notebook sau.","metadata":{}}]}