{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":193361348,"sourceType":"kernelVersion"}],"dockerImageVersionId":30746,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# The issue\n\nAs @ambrosm has show in [this discussion](https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/527167), the transit time for the planets to pass the stars is not the same for each sample. Though, calculation and feature engineering on the light curves requires exact knowledge of these cutoffs.\n\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F7917824%2F9b51760603a8f72d3604305ed76b8841%2Fsl.png?generation=1723307476416972&alt=media)\n\n# The solution\n\nHere, I am presenting a very robust solution to identify the transit zone that works on ALL train samples, with both sensor types. Discussion about this solution here: https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/529533\n\n![](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F2675447%2Fc8dcc7b6800036ce45a352d298c99249%2F__results___6_1.png?generation=1724246305869502&alt=media)","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-09-02T07:14:21.274114Z","iopub.execute_input":"2024-09-02T07:14:21.274495Z","iopub.status.idle":"2024-09-02T07:14:22.76262Z","shell.execute_reply.started":"2024-09-02T07:14:21.274463Z","shell.execute_reply":"2024-09-02T07:14:22.761339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ROOT = \"/kaggle/input/ariel-data-challenge-2024/\"\n\ntrain_adc_info = pd.read_csv(ROOT + \"/train_adc_info.csv\", index_col=\"planet_id\")\ntrain_labels = pd.read_csv(ROOT + \"/train_labels.csv\", index_col=\"planet_id\")","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:14:22.764874Z","iopub.execute_input":"2024-09-02T07:14:22.765518Z","iopub.status.idle":"2024-09-02T07:14:22.903892Z","shell.execute_reply.started":"2024-09-02T07:14:22.765478Z","shell.execute_reply":"2024-09-02T07:14:22.901818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:14:22.906035Z","iopub.execute_input":"2024-09-02T07:14:22.906535Z","iopub.status.idle":"2024-09-02T07:14:22.952678Z","shell.execute_reply.started":"2024-09-02T07:14:22.906498Z","shell.execute_reply":"2024-09-02T07:14:22.951506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_raw_train = np.load(\"/kaggle/input/ariel24-data-prep/fgs_v4.npy\")\na_raw_train = np.load(\"/kaggle/input/ariel24-data-prep/airs_v4.npy\")","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:14:22.955798Z","iopub.execute_input":"2024-09-02T07:14:22.956311Z","iopub.status.idle":"2024-09-02T07:14:28.281343Z","shell.execute_reply.started":"2024-09-02T07:14:22.956251Z","shell.execute_reply":"2024-09-02T07:14:28.280186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_raw_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:14:28.282748Z","iopub.execute_input":"2024-09-02T07:14:28.283196Z","iopub.status.idle":"2024-09-02T07:14:28.291683Z","shell.execute_reply.started":"2024-09-02T07:14:28.283157Z","shell.execute_reply":"2024-09-02T07:14:28.290303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.signal import savgol_filter\n\n\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  # window size 51, polynomial order 3\n\n\ndef optimize_breakpoint(data, initial_breakpoint, window_size=500, buffer_size=50, smooth_window=121):\n    best_breakpoint = initial_breakpoint\n    best_score = float(\"-inf\")\n    midpoint = len(data) // 2\n    smoothed_data = smooth_data(data, smooth_window)\n\n    for i in range(-window_size, window_size):\n        new_breakpoint = initial_breakpoint + i\n        region1 = data[: new_breakpoint - buffer_size]\n        region2 = data[\n            new_breakpoint\n            + buffer_size : len(data)\n            - new_breakpoint\n            - buffer_size\n        ]\n        region3 = data[len(data) - new_breakpoint + buffer_size :]\n\n        # calc on smoothed data\n        breakpoint_region1 = smoothed_data[new_breakpoint - buffer_size: new_breakpoint + buffer_size]\n        breakpoint_region2 = smoothed_data[-(new_breakpoint + buffer_size): -(new_breakpoint - buffer_size)]\n\n        mean_diff = abs(np.mean(region1) - np.mean(region2)) + abs(\n            np.mean(region2) - np.mean(region3)\n        )\n        var_sum = np.var(region1) + np.var(region2) + np.var(region3)\n        range_at_breakpoint1 = (np.max(breakpoint_region1) - np.min(breakpoint_region1))\n        range_at_breakpoint2 = (np.max(breakpoint_region2) - np.min(breakpoint_region2))\n\n        mean_range_at_breakpoint = (range_at_breakpoint1 + range_at_breakpoint2) / 2\n\n        score = mean_diff - 0.5 * var_sum + mean_range_at_breakpoint\n\n        if score > best_score:\n            best_score = score\n            best_breakpoint = new_breakpoint\n\n    return best_breakpoint\n\n\n\ndef find_and_plot_breakpoints(f_raw, a_raw, IDX, verbose=False):\n    # Create a figure with two subplots\n    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 6))\n    planet_id = train_labels.index[IDX]\n\n    for sensor_name in [\"AIRS-CH0\", \"FGS1\"]:\n        if sensor_name == \"AIRS-CH0\":\n            buffer_size = 120\n            smooth_window = 121\n            window_size = 500\n            default_breakpoint = 1821\n            data = a_raw\n        elif sensor_name == \"FGS1\":\n            buffer_size = 1500\n            smooth_window = 1001\n            window_size = 3500\n            default_breakpoint = 21500\n            data = f_raw\n\n        initial_breakpoint = default_breakpoint\n\n        if verbose:\n            print(f\"Initial breakpoint: {initial_breakpoint}\")\n        optimized_breakpoint = optimize_breakpoint(\n            data,\n            initial_breakpoint,\n            window_size=window_size,\n            buffer_size=buffer_size,\n            smooth_window=smooth_window,\n        )\n\n        midpoint = len(data) // 2\n        breakpoints = [optimized_breakpoint, 2 * midpoint - optimized_breakpoint]\n\n        if sensor_name == \"FGS1\":\n            ax = ax2\n        elif sensor_name == \"AIRS-CH0\":\n            ax = ax1\n\n        # Plot the results\n        ax.plot(data, color='#4E79A7', alpha=0.7, label=\"Original Data\")\n        ax.plot(smooth_data(data, smooth_window), label=\"Smoothed Data\", color='#F28E2B')\n\n        for bp in breakpoints:\n            ax.axvline(x=bp - buffer_size, color=\"r\", linestyle=\"--\")\n            ax.axvline(x=bp + buffer_size, color=\"r\", linestyle=\"--\")\n            ax.axvspan(bp - buffer_size, bp + buffer_size, color=\"gray\", alpha=0.3)\n\n        ax.set_title(f\"{sensor_name}\")\n        ax.set_xlabel(\"Time step\")\n        ax.set_ylabel(\"Value\")\n        ax.legend()\n    fig.tight_layout()\n    fig.suptitle(f\"Planet ID: {planet_id}, IDX {IDX}\", fontsize=16)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:16:44.933053Z","iopub.execute_input":"2024-09-02T07:16:44.933486Z","iopub.status.idle":"2024-09-02T07:16:44.955097Z","shell.execute_reply.started":"2024-09-02T07:16:44.933453Z","shell.execute_reply":"2024-09-02T07:16:44.953681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for IDX in range(20):\n    f_raw = f_raw_train[IDX]\n    a_raw = a_raw_train[IDX]\n    find_and_plot_breakpoints(f_raw, a_raw, IDX)","metadata":{"execution":{"iopub.status.busy":"2024-09-02T07:16:45.988541Z","iopub.execute_input":"2024-09-02T07:16:45.98909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}