{"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":"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\nimport string\nimport math\nfrom IPython.display import clear_output\n\nRANDOM_SEED = 27\n\ndata_dir = '../input/vsb-power-line-fault-detection'","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-02T12:16:16.972346Z","iopub.execute_input":"2023-11-02T12:16:16.973303Z","iopub.status.idle":"2023-11-02T12:16:16.981691Z","shell.execute_reply.started":"2023-11-02T12:16:16.973255Z","shell.execute_reply":"2023-11-02T12:16:16.979705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 加载CSV文件\nmeta_train_df = pd.read_csv(data_dir + '/metadata_train.csv')\nmeta_test_df = pd.read_csv(data_dir + '/metadata_test.csv')\nmeta_train_df.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:16:16.984417Z","iopub.execute_input":"2023-11-02T12:16:16.984866Z","iopub.status.idle":"2023-11-02T12:16:17.025971Z","shell.execute_reply.started":"2023-11-02T12:16:16.984825Z","shell.execute_reply":"2023-11-02T12:16:17.024837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 加载数据集，选择样本\ntrain_df = pd.read_parquet(data_dir + '/train.parquet')\n\n# 选择五组样本\nnegative_signal_ids = meta_train_df[meta_train_df.id_measurement==2].signal_id.values #在此修改样本编码\npositive_signal_ids = meta_train_df[meta_train_df.id_measurement==1].signal_id.values\n\nnegative_sample = train_df.iloc[:, negative_signal_ids].values\npositive_sample = train_df.iloc[:, positive_signal_ids].values\n\nplt.figure(figsize=(12, 4))\nplt.title('Normal powerline')\nplt.plot(negative_sample, alpha=0.8);\nplt.savefig('Normal1.png')\n\nplt.figure(figsize=(12, 4))\nplt.title('Faulty powerline')\nplt.plot(positive_sample, alpha=0.8);\nplt.savefig('Faulty1.png')\n","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:16:17.027611Z","iopub.execute_input":"2023-11-02T12:16:17.028027Z","iopub.status.idle":"2023-11-02T12:17:26.228030Z","shell.execute_reply.started":"2023-11-02T12:16:17.027979Z","shell.execute_reply":"2023-11-02T12:17:26.226869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 相位对齐\ndef zero_phase_dft(signal):\n    # 计算DFT\n    dft_result = np.fft.fft(signal)\n    \n    # 计算相位角\n    phase = np.angle(dft_result)\n    \n    # 将第一个时间步长的相位调整为零度\n    phase[1] -= phase[1]\n    \n    # 重建信号\n    processed_signal = np.fft.ifft(np.abs(dft_result) * np.exp(1j * phase))\n    \n    return np.real(processed_signal)\n\n# 使用EMA残差进行信号平坦化\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\n    \nnegative_sample = zero_phase_dft(negative_sample)\npositive_sample = zero_phase_dft(positive_sample)\n\n\nflat_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])\nplt.figure(figsize=(12, 4))\nplt.title('Normal powerline')\nplt.plot(flat_negative_sample, alpha=0.8);\nplt.savefig('flat_Normal1.png')\n\nplt.figure(figsize=(12, 4))\nplt.title('Faulty powerline')\nplt.plot(flat_positive_sample, alpha=0.8);\nplt.savefig('flat_Faulty1.png')","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:17:26.229456Z","iopub.execute_input":"2023-11-02T12:17:26.230900Z","iopub.status.idle":"2023-11-02T12:17:38.901285Z","shell.execute_reply.started":"2023-11-02T12:17:26.230865Z","shell.execute_reply":"2023-11-02T12:17:38.900193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 噪声水平估计函数\ndef noise_level_estimation(s, Nnoise, Lnoise, Ncover, Cmax, Cstep):\n    # Step 1: Sample Nnoise sections from s\n    samples = []\n    segment_length = len(s) // Nnoise\n    for i in range(Nnoise):\n        sample_start = i * segment_length\n        sample_end = (i + 1) * segment_length\n        sample = s[sample_start:sample_end]\n        samples.append(sample)\n        \n    # Step 2: Calculate max absolute values for each sampled section\n    max_values = [max(abs(sample)) for sample in samples]\n    \n    # Step 3: Loop through possible noise levels\n    i=Cmax\n    while i<=Cmax:\n        k = 0\n        for j in range(Nnoise):\n            m = max_values[j]\n            if i> m >= (i - 1):\n                k += 1\n                if k > Ncover:\n                    a = i\n                    print(a)\n                    return a\n        i+=Cstep\n        \n    # Step 4: If no suitable noise level found, return default\n    return None\n\n#缩放信号\ndef scale_signal(signal, target_threshold=5):\n    # 获取阈值\n    a = noise_level_estimation(signal, Nnoise=1000, Lnoise=1000, Ncover=80, Cmax=15, Cstep=-0.5)\n    \n    # 缩放信号\n    scaled_signal = signal * (target_threshold / a)\n    \n    return scaled_signal\n\n# 噪声估计A\nflat_negative_sample[:,0] = scale_signal(flat_negative_sample[:,0])\nflat_positive_sample[:,0] = scale_signal(flat_positive_sample[:,0])\n    \n# 绘制并保存\nplt.figure(figsize=(12, 4))\nplt.title(f'Normal powerline Phase A')\nplt.plot(flat_negative_sample[:,0], alpha=0.8)\nplt.axhline(y=5, color='r', linestyle='--', label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\nplt.savefig(f'nle_Normal1_A.png')\n\nplt.figure(figsize=(12, 4))\nplt.title(f'Faulty powerline Phase A')\nplt.plot(flat_positive_sample[:,0], alpha=0.8)\nplt.axhline(y=5, color='r', linestyle='--',label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\nplt.savefig(f'nle_Faulty1_A.png')\n\n# 噪声估计B\nflat_negative_sample[:,1] = scale_signal(flat_negative_sample[:,1])\nflat_positive_sample[:,1] = scale_signal(flat_positive_sample[:,1])\n    \n# 绘制并保存\nplt.figure(figsize=(12, 4))\nplt.title(f'Normal powerline Phase B')\nplt.plot(flat_negative_sample[:,1],alpha=0.8,color='orange')\nplt.axhline(y=5, color='r', linestyle='--', label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\nplt.savefig(f'nle_Normal1_B.png')\n\nplt.figure(figsize=(12, 4))\nplt.title(f'Faulty powerline Phase B')\nplt.plot(flat_positive_sample[:,1],alpha=0.8,color='orange')\nplt.axhline(y=5, color='r', linestyle='--',label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\n\nplt.legend()\nplt.savefig(f'nle_Faulty1_B.png')\n\n# 噪声估计C\nflat_negative_sample[:,2] = scale_signal(flat_negative_sample[:,2])\nflat_positive_sample[:,2] = scale_signal(flat_positive_sample[:,2])\n    \n# 绘制并保存\nplt.figure(figsize=(12, 4))\nplt.title(f'Normal powerline Phase C')\nplt.plot(flat_negative_sample[:,2],alpha=0.8,color='g')\nplt.axhline(y=5, color='r', linestyle='--', label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\nplt.savefig(f'nle_Normal1_C.png')\n\nplt.figure(figsize=(12, 4))\nplt.title(f'Faulty powerline Phase C')\nplt.plot(flat_positive_sample[:,2],alpha=0.8,color='g')\nplt.axhline(y=5, color='r', linestyle='--',label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\nplt.savefig(f'nle_Faulty1_C.png')","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:17:38.904316Z","iopub.execute_input":"2023-11-02T12:17:38.905003Z","iopub.status.idle":"2023-11-02T12:17:48.220591Z","shell.execute_reply.started":"2023-11-02T12:17:38.904962Z","shell.execute_reply":"2023-11-02T12:17:48.219234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\"import os\nimport shutil\nfrom IPython.display import FileLink\n# 创建临时文件夹（如果不存在）\nos.makedirs('/kaggle/working/tem', exist_ok=True)\n\n# 遍历当前工作目录中的所有文件和文件夹\nfor root, dirs, files in os.walk('/kaggle/working'):\n    if root != '/kaggle/working/tem':\n        for file in files:\n            if file.endswith(('.png', '.jpg', '.jpeg', '.gif')):\n                # 构建源文件路径和目标文件路径\n                source_path = os.path.join(root, file)\n                target_path = os.path.join('/kaggle/working/tem', file)\n                # 复制文件\n                shutil.copy2(source_path, target_path)\nshutil.make_archive('/kaggle/working/tem', 'zip', '/kaggle/working/tem')\n\n# 创建一个链接以下载文件\nFileLink(r'/kaggle/working/tem.zip')\"\"\"","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:17:48.222353Z","iopub.execute_input":"2023-11-02T12:17:48.222802Z","iopub.status.idle":"2023-11-02T12:17:48.231785Z","shell.execute_reply.started":"2023-11-02T12:17:48.222763Z","shell.execute_reply":"2023-11-02T12:17:48.230660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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\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\ndef local_maxima_1d_window(x, w=1,threshold=5):\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    # 仅保留大于阈值的峰值\n    m = m[x[m] > threshold]\n\n    return m\n\n# 拐点检测函数\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\n\n\ndef get_peaks(x_flatten, window=25,visualise=False,visualise_color=None,):\n    \n    # 寻找局部最大峰值，并记录。\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    # 排序\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    # 平滑处理计算斜率\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    # 计算拐点\n    knee_x = plateau_detection(grad, -0.01, plateau_length=1000)\n    knee_x -= conv_length//2\n    \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    ii = np.argsort(peaks_x)\n    \n    peaks_x = peaks_x[ii]\n    peaks_y = peaks_y[ii]\n        \n    return peaks_x, peaks_y\n\n# 噪声估计A\nN_px, N_py = get_peaks(flat_negative_sample[:,0])\nF_px, F_py = get_peaks(flat_positive_sample[:,0])\n    \n# 绘制并保存\nplt.figure(figsize=(12, 4))\nplt.title(f'Normal powerline Phase A')\nplt.plot(flat_negative_sample[:,0], alpha=0.8)\nplt.axhline(y=5, color='r', linestyle='--', label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\n\n# 绘制峰值\nplt.scatter(N_px, flat_negative_sample[N_px,0], color=\"red\")\nplt.savefig(f'nle_Normal4_A.png')\n\nplt.figure(figsize=(12, 4))\nplt.title(f'Faulty powerline Phase A')\nplt.plot(flat_positive_sample[:,0], alpha=0.8)\nplt.axhline(y=5, color='r', linestyle='--',label='Threshold')\nplt.axhline(y=-5, color='r', linestyle='--')\nplt.legend()\n\n# 绘制峰值\nplt.scatter(F_px, flat_positive_sample[F_px,0], color=\"red\")\nplt.savefig(f'nle_Faulty4_A.png')\n","metadata":{"execution":{"iopub.status.busy":"2023-11-02T12:17:48.233629Z","iopub.execute_input":"2023-11-02T12:17:48.234052Z","iopub.status.idle":"2023-11-02T12:17:52.805490Z","shell.execute_reply.started":"2023-11-02T12:17:48.234013Z","shell.execute_reply":"2023-11-02T12:17:52.804208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}