{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"},{"sourceId":192766898,"sourceType":"kernelVersion"}],"dockerImageVersionId":30747,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\n# Set to \"1\" to directly force test df, this is much quicker to commit\nos.environ[\"KAGGLE_IS_COMPETITION_RERUN\"] = \"1\"\n\nMODELS_LOAD = False  # True\nPATH_MODELS = '/kaggle/input/ariel-2025-02-01'\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    MODELS_LOAD = False\n","metadata":{"execution":{"iopub.status.busy":"2025-07-08T14:03:08.433536Z","iopub.execute_input":"2025-07-08T14:03:08.433783Z","iopub.status.idle":"2025-07-08T14:03:22.118401Z","shell.execute_reply.started":"2025-07-08T14:03:08.433755Z","shell.execute_reply":"2025-07-08T14:03:22.117337Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile preprocess.py\n\nimport os\nimport pandas as pd\nimport numpy as np\n#import itertools\nfrom tqdm import tqdm\nimport multiprocessing as mp\nfrom astropy.stats import sigma_clip\nimport torch\nimport torch.nn.functional as F\n\nROOT = \"/kaggle/input/ariel-data-challenge-2025/\"\nBINNING = 15\n# 16 center pixels, rest is just noise\ncl = 8\ncr = 24\nsensor_sizes_dict = {\n    \"AIRS-CH0\": [[11250, 32, 356], [32, 356]],\n    \"FGS1\": [[135000, 32, 32], [32, 32]],\n}  # input, mask\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    MODE = \"test\"\n    MODELS_LOAD = False\nelse:\n    MODE = \"train\"\n\ndef get_gain_offset():\n    \"\"\"\n    Get the gain and offset for a given planet and sensor\n\n    Unlike last year's challenge, all planets use the same adc_info.\n    We can just hard code it.\n    \"\"\"\n    gain = 0.4369\n    offset = -1000.0\n    return gain, offset\n\ndef read_data(planet_id, sensor, mode):\n    \"\"\"\n    Read the data for a given planet and sensor\n    \"\"\"\n    # get all noise correction frames and signal\n    signal = pd.read_parquet(\n        f\"{ROOT}/{mode}/{planet_id}/{sensor}_signal_0.parquet\",\n        engine=\"pyarrow\",\n    )\n    dark_frame = pd.read_parquet(\n        f\"{ROOT}/{mode}/{planet_id}/{sensor}_calibration_0/dark.parquet\",\n        engine=\"pyarrow\",\n    )\n    dead_frame = pd.read_parquet(\n        f\"{ROOT}/{mode}/{planet_id}/{sensor}_calibration_0/dead.parquet\",\n        engine=\"pyarrow\",\n    )\n    linear_corr_frame = pd.read_parquet(\n        f\"{ROOT}/{mode}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\",\n        engine=\"pyarrow\",\n    )\n    flat_frame = pd.read_parquet(\n        f\"{ROOT}/{mode}/{planet_id}/{sensor}_calibration_0/flat.parquet\",\n        engine=\"pyarrow\",\n    )\n    # read_frame = pd.read_parquet(\n    #     f\"{ROOT}/{mode}/{planet_id}/{sensor}_calibration/read.parquet\",\n    #     engine=\"pyarrow\",\n    # )\n    # reshape to sensor shape and cast to float64\n    signal = signal.values.astype(np.float64).reshape(sensor_sizes_dict[sensor][0])[\n        :, cl:cr, :\n    ]\n    dark_frame = dark_frame.values.astype(np.float64).reshape(\n        sensor_sizes_dict[sensor][1]\n    )[cl:cr, :]\n    dead_frame = dead_frame.values.reshape(sensor_sizes_dict[sensor][1])[cl:cr, :]\n    flat_frame = flat_frame.values.astype(np.float64).reshape(\n        sensor_sizes_dict[sensor][1]\n    )[cl:cr, :]\n    # read_frame = read_frame.values.reshape(sensor_sizes_dict[sensor][1])\n    linear_corr = linear_corr_frame.values.astype(np.float64).reshape(\n        [6] + sensor_sizes_dict[sensor][1]\n    )[:, cl:cr, :]\n    return (\n        signal,\n        dark_frame,\n        dead_frame,\n        linear_corr,\n        flat_frame,\n        # read_frame,\n    )\n\ndef ADC_convert(signal, gain, offset):\n    \"\"\"\n    Step 1: Analog-to-Digital Conversion (ADC) correction\n    The Analog-to-Digital Conversion (adc) is performed by the detector to convert the\n    pixel voltage into an integer number. We revert this operation by using the gain\n    and offset for the calibration files 'train_adc_info.csv'.\n    \"\"\"\n    return signal / gain + offset\n\ndef mask_hot_dead(signal, dead, dark):\n    \"\"\"\n    Step 2: Mask hot/dead pixel\n    The dead pixels map is a map of the pixels that do not respond to light and, thus,\n    can't be accounted for any calculation. In all these frames the dead pixels are\n    masked using python masked arrays. The bad pixels are thus masked but left\n    uncorrected. Some methods can be used to correct bad-pixels but this task,\n    if needed, is left to the participants.\n    \"\"\"\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    # Set values to np.nan where dead or hot pixels are found\n    signal[dead] = np.nan\n    signal[hot] = np.nan\n    return signal\n\ndef apply_linear_corr(c, signal):\n    \"\"\"\n    Step 3: linearity Correction\n    The non-linearity of the pixels' response can be explained as capacitive leakage\n    on the readout electronics of each pixel during the integration time. The number\n    of electrons in the well is proportional to the number of photons that hit the\n    pixel, with a quantum efficiency coefficient. However, the response of the pixel\n    is not linear with the number of electrons in the well. This effect can be\n    described by a polynomial function of the number of electrons actually in the well.\n    The data is provided with calibration files linear_corr.parquet that are the\n    coefficients of the inverse polynomial function and can be used to correct this\n    non-linearity effect.\n    Using horner's method to evaluate the polynomial\n    \"\"\"\n    assert c.shape[0] == 6  # Ensure the polynomial is of degree 5\n    return (\n        (((c[5] * signal + c[4]) * signal + c[3]) * signal + c[2]) * signal + c[1]\n    ) * signal + c[0]\n\ndef clean_dark(signal, dark, dt):\n    \"\"\"\n    Step 4: dark current subtraction\n    The data provided include calibration for dark current estimation, which can be\n    used to pre-process the observations. Dark current represents a constant signal\n    that accumulates in each pixel during the integration time, independent of the\n    incoming light. To obtain the corrected image, the following conventional approach\n    is applied: The data provided include calibration files such as dark frames or\n    dead pixels' maps. They can be used to pre-process the observations. The dark frame\n    is a map of the detector response to a very short exposure time, to correct for the\n    dark current of the detector.\n    image - (dark * dt)\n    The corrected image is conventionally obtained via the following: where the dark\n    current map is first corrected for the dead pixel.\n    \"\"\"\n    dark = torch.tile(dark, (signal.shape[0], 1, 1))\n    signal -= dark * dt[:, None, None]\n    return signal\n\ndef get_cds(signal):\n    \"\"\"\n    Step 5: Get Correlated Double Sampling (CDS)\n    The science frames are alternating between the start of the exposure and the end of\n    the exposure. The lecture scheme is a ramp with a double sampling, called\n    Correlated Double Sampling (CDS), the detector is read twice, once at the start\n    of the exposure and once at the end of the exposure. The final CDS is the\n    difference (End of exposure) - (Start of exposure).\n    \"\"\"\n    return torch.subtract(signal[1::2, :, :], signal[::2, :, :])\n\ndef bin_obs(signal, binning):\n    \"\"\"\n    Step 5.1: Bin Observations\n    The data provided are binned in the time dimension. The binning is performed by\n    summing the signal over the time dimension.\n    \"\"\"\n    assert signal.shape[0] % binning == 0  # Ensure the binning is possible\n    # cds_transposed = signal.transpose(0, 2, 1)\n    cds_binned = torch.zeros((signal.shape[0] // binning, signal.shape[1], signal.shape[2]), device=\"cuda:0\")\n    for i in range(signal.shape[0] // binning):\n        cds_binned[i, :, :] = torch.sum(signal[i * binning : (i + 1) * binning, :, :], axis=0)\n    return cds_binned\n\ndef correct_flat_field(flat, signal):\n    \"\"\"\n    Step 6: Flat Field Correction\n    The flat field is a map of the detector response to uniform illumination, to\n    correct for the pixel-to-pixel variations of the detector, for example the\n    different quantum efficiencies of each pixel.\n    \"\"\"\n    return signal / flat\n\ndef nan_interpolation(tensor):\n    # Assume tensor is of shape (batch, height, width)\n    nan_mask = torch.isnan(tensor)\n    # Replace NaNs with zero temporarily\n    tensor_filled = torch.where(nan_mask, torch.tensor(0.0, device=tensor.device), tensor)\n    # Create a binary mask (0 where NaNs were and 1 elsewhere)\n    ones = torch.ones_like(tensor, device=tensor.device)\n    weight = torch.where(nan_mask, torch.tensor(0.0, device=tensor.device), ones)\n    # Perform interpolation by convolving with a kernel\n    # using bilinear interpolation\n    kernel = torch.ones(1, 1, 1, 3, device=tensor.device, dtype=tensor.dtype)\n    # Apply padding to the tensor and weight to prevent boundary issues\n    tensor_padded = F.pad(tensor_filled.unsqueeze(1), (1, 1, 0, 0), mode=\"replicate\").squeeze(1)\n    weight_padded = F.pad(weight.unsqueeze(1), (1, 1, 0, 0), mode=\"replicate\").squeeze(1)\n    # Convolve the filled tensor and the weight mask\n    tensor_conv = F.conv2d(tensor_padded.unsqueeze(1), kernel, stride=1)\n    weight_conv = F.conv2d(weight_padded.unsqueeze(1), kernel, stride=1)\n    # Compute interpolated values (normalized by weights)\n    interpolated_tensor = tensor_conv / weight_conv\n    # Apply the interpolated values only to the positions of NaNs\n    result = torch.where(nan_mask, interpolated_tensor.squeeze(1), tensor)\n    return result\n\ndef process_planet(planet_id):\n    \"\"\"\n    Process a single planet's data\n    \"\"\"\n    axis_info = pd.read_parquet(ROOT + \"axis_info.parquet\")  # (135000, 4)\n    dt_airs = axis_info[\"AIRS-CH0-integration_time\"].dropna().values  # (11250,)\n    for sensor in [\"FGS1\", \"AIRS-CH0\"]:\n        # load all data for this planet and sensor\n        signal, dark_frame, dead_frame, linear_corr, flat_frame = read_data(planet_id, sensor, mode=MODE)  \n        # FGS1=[signal=(135000, 16, 32), dark_frame=(16, 32), dead_frame=(16, 32), linear_corr=(6, 16, 32), flat_frame=(16, 32)] # AIRS0=((11250, 16, 356), (16, 356), (16, 356), (6, 16, 356), (16, 356))\n        gain, offset = get_gain_offset()\n        # Step 1: ADC correction\n        signal = ADC_convert(signal, gain, offset)  # FGS1=(135000, 16, 32) / AIRS0=(11250, 16, 356)\n        # Step 2: Mask hot/dead pixel\n        signal = mask_hot_dead(signal, dead_frame, dark_frame)  # FGS1=(135000, 16, 32) / AIRS0=(11250, 16, 356)\n        # clip at 0\n        signal = signal.clip(0)\n        # Step 3: linearity Correction\n        signal = apply_linear_corr(torch.tensor(linear_corr).to(\"cuda:0\"), torch.tensor(signal).to(\"cuda:0\"))  # FGS1=(135000, 16, 32) / AIRS0=[11250, 16, 356]\n        # Step 4: dark current subtraction\n        if sensor == \"FGS1\":\n            dt = torch.ones(len(signal), device=\"cuda:0\") * 0.1  # FGS1=[135000]\n            dt[1::2] += 4.5\n        elif sensor == \"AIRS-CH0\":\n            dt = torch.tensor(dt_airs).to(\"cuda:0\")  # AIRS0=[11250]\n            dt[1::2] += 0.1\n        signal = clean_dark(signal, torch.tensor(dark_frame).to(\"cuda:0\"), dt)  # FGS1=[135000, 16, 32] / AIRS0=[11250, 16, 356]\n        # Step 5: Get Correlated Double Sampling (CDS)\n        signal = get_cds(signal)  # FGS1=[67500, 16, 32] / AIRS0=[5625, 16, 356]\n        # Step 5.1: Bin Observations\n        if sensor == \"FGS1\":\n            signal = bin_obs(signal, binning=BINNING * 12)  # FGS1=[375, 16, 32]\n        elif sensor == \"AIRS-CH0\":\n            signal = bin_obs(signal, binning=BINNING)  # AIRS0=[375, 16, 356]\n        # Step 6: Flat Field Correction\n        signal = correct_flat_field(torch.tensor(flat_frame).to(\"cuda:0\"), signal)  # FGS1=[375, 16, 32] / AIRS0=[375, 16, 356]\n        # Step 7: Interpolate NaNs (twice!)\n        signal = nan_interpolation(signal)  # FGS1=[375, 16, 32] / AIRS0=[375, 16, 356]\n        signal = nan_interpolation(signal)  # FGS1=[375, 16, 32] / AIRS0=[375, 16, 356]\n        # Step 8: Sum over spatial axis\n        if sensor == \"FGS1\":\n            signal = torch.nanmean(signal, axis=[1, 2]).cpu().numpy()  # FGS1=(375,)\n        elif sensor == \"AIRS-CH0\":\n            signal = torch.nanmean(signal, axis=1).cpu().numpy()       # AIRS0=[375, 356]\n        # save the processed signal\n        np.save(f\"{planet_id}_{sensor}_signal.npz\", signal.astype(np.float64))\n\nif __name__ == \"__main__\":\n    star_info = pd.read_csv(ROOT + f\"/{MODE}_star_info.csv\")\n    star_info[\"planet_id\"] = star_info[\"planet_id\"].astype(int)\n    star_info = star_info.set_index(\"planet_id\")\n    planet_ids = star_info.index.tolist()\n    # process_planet(planet_ids[0])\n    with mp.Pool(processes=4) as pool:\n        list(tqdm(pool.imap(process_planet, planet_ids), total=len(planet_ids)))\n    signal_train = []\n    for i, planet_id in enumerate(planet_ids):\n        f_raw = np.load(f\"{planet_id}_FGS1_signal.npz.npy\")\n        a_raw = np.load(f\"{planet_id}_AIRS-CH0_signal.npz.npy\")\n        # flip a_raw\n        signal = np.concatenate([f_raw[:, None], a_raw[:, ::-1]], axis=1)\n        signal_train.append(signal)\n    signal_train = np.array(signal_train)\n    np.save(\"/kaggle/working/signal.npy\", signal_train, allow_pickle=False)\n    print(\"Processing complete!\")\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2025-07-08T14:04:34.388352Z","iopub.execute_input":"2025-07-08T14:04:34.389096Z","iopub.status.idle":"2025-07-08T14:04:34.400681Z","shell.execute_reply.started":"2025-07-08T14:04:34.389067Z","shell.execute_reply":"2025-07-08T14:04:34.39964Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not MODELS_LOAD:\n    ! python preprocess.py","metadata":{"execution":{"iopub.status.busy":"2025-07-08T14:04:48.215321Z","iopub.execute_input":"2025-07-08T14:04:48.215606Z","iopub.status.idle":"2025-07-08T14:04:49.251098Z","shell.execute_reply.started":"2025-07-08T14:04:48.215581Z","shell.execute_reply":"2025-07-08T14:04:49.250164Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"! rm -rf *FGS1_signal*\n! rm -rf *AIRS-CH0_signal*\n! ls","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\n#from sklearn.linear_model import Ridge\n#from sklearn.linear_model import LinearRegression\nfrom scipy.signal import savgol_filter\nfrom scipy.optimize import curve_fit\nfrom tqdm import tqdm\n\nSMOOTH_WINDOW = 19\nPRE_BINNED_TIME = 15\nBINNING = 15\nbuffer_size_poly = 150 // PRE_BINNED_TIME  # 10\nROOT = \"/kaggle/input/ariel-data-challenge-2025/\"\n\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    MODE = \"test\"\nelse:\n    MODE = \"train\"\n\nstar_info = pd.read_csv(ROOT + f\"/{MODE}_star_info.csv\")\nstar_info[\"planet_id\"] = star_info[\"planet_id\"].astype(int)\nstar_info = star_info.set_index(\"planet_id\")\nwavelengths = pd.read_csv(ROOT + \"/wavelengths.csv\")\n\nsignal_train = np.load(f\"{PATH_MODELS}/signal.npy\" if MODELS_LOAD else \"/kaggle/working/signal.npy\")\ncut_inf, cut_sup = 36, 318\nsignal_train = np.concatenate([signal_train[:, :, 0][:, :, None], signal_train[:, :, cut_inf:cut_sup]], axis=2)\nprint(signal_train.shape)\n# signal_train = signal_train.mean(axis=2)\n\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  # window size 51, polynomial order 3\n\n# find transit zones\ndef phase_detector(signal_orig, binning=15, smooth_window=11, verbose=False):\n    signal = signal_orig.reshape(-1, binning).mean(-1)  # collapse by 15; 375\n    signal = savgol_filter(signal, smooth_window, 2)  # smooth\n    first_derivative = np.gradient(signal)\n    phase1 = np.argmin(first_derivative)\n    phase2 = np.argmax(first_derivative)\n    if verbose:\n        plt.plot(signal_orig, color=\"grey\", alpha=0.5, label=\"original\")\n        plt.plot(signal, color=\"blue\", alpha=0.9, label=\"smoothed\")\n        plt.axvline(phase1, color=\"r\")\n        plt.axvline(phase2, color=\"r\")\n        plt.show()\n        plt.plot(first_derivative, color=\"green\", alpha=0.9, label=\"first derivative\")\n        plt.show()\n    #assert phase1 < phase2\n    #assert phase1 >= 0\n    #assert phase2 <= signal.shape[0]\n    return phase1 * binning, phase2 * binning\n\ndef get_breakpoints(x, pre_binned_time, verbose=False):\n    bp = np.zeros(x.shape[0], dtype=np.int32)  # (1100,)\n    bp2 = np.zeros(x.shape[0], dtype=np.int32)  # (1100,)\n    for i in range(x.shape[0]):\n        signal = x[i].mean(-1)   # signal=(1100, 375, 283) => (375,)\n        p1, p2 = phase_detector(signal, binning=BINNING // pre_binned_time, smooth_window=SMOOTH_WINDOW, verbose=verbose)\n        bp[i] = p1\n        bp2[i] = p2\n    return [bp, bp2]\n\ndef poly_exp_fit(data, optimized_breakpoints, buffer_size, degree=3):\n    # Define the three regions\n    low_break = max(1, optimized_breakpoints[0] - buffer_size)\n    high_break = min(len(data)-1, optimized_breakpoints[1] + buffer_size)\n    if low_break >= high_break:\n        raise Exception(\"e001\")\n    if optimized_breakpoints[0] + buffer_size >= optimized_breakpoints[1] - buffer_size:\n        raise Exception(\"e002\")\n    x1 = np.arange(low_break)  # (120,)\n    y1 = data[: low_break]  # (120,)\n    x2 = np.arange(optimized_breakpoints[0] + buffer_size, optimized_breakpoints[1] - buffer_size)  # (122,)\n    y2 = data[optimized_breakpoints[0] + buffer_size : optimized_breakpoints[1] - buffer_size]  # (122,)\n    x3 = np.arange(high_break, len(data))  # (93,)\n    y3 = data[high_break :]  # (93,)\n    # Concatenate the x-values and y-values for regions 1 and 3\n    x_combined = np.concatenate([x1, x3])  # (213,)\n    y_combined = np.concatenate([y1, y3])  # (213,)\n\n    def fit_function(x, *params):\n        poly_params = params[: degree + 1]\n        y_fit = np.polyval(poly_params, x)\n        return y_fit\n\n    # Define the polynomial fit function with an additional shift parameter for region 2\n    def fit_function_with_shift(x, shift, *poly_params):\n        x1_adjusted = x[: len(x1)]\n        x2_adjusted = x[len(x1) : len(x1) + len(x2)]\n        x3_adjusted = x[len(x1) + len(x2) :]\n        y1_fit = np.polyval(poly_params, x1_adjusted)\n        y2_fit = np.polyval(poly_params, x2_adjusted) * shift\n        y3_fit = np.polyval(poly_params, x3_adjusted)\n        return np.concatenate([y1_fit, y2_fit, y3_fit])\n\n    # Define the combined x-values (including region 2)\n    x_combined_with_region2 = np.concatenate([x1, x2, x3])  # (335,)\n    y_combined_with_region2 = np.concatenate([y1, y2, y3])  # (335,)\n    # Initial guesses for the polynomial coefficients and shift\n    try:\n        poly_guess = np.polyfit(x_combined, y_combined, degree)  # (3,)\n    except Exception as error:\n        print(f\"{error}\")\n        pass  # x1.shape, y1.shape, x2.shape, y2.shape, x3.shape, y3.shape,\n              # ((0,), (365,), (329,), (329,), (16,), (16,))\n              # x_combined.shape, y_combined.shape, degree => ((16,), (381,), 2)\n              # optimized_breakpoints => [0, 349]\n    p0 = list(poly_guess)\n    initial_shift_guess = 1.0\n    p0 = [initial_shift_guess] + list(p0)\n    # Fit the polynomial and the shift using curve_fit\n    popt, _ = curve_fit(\n        fit_function_with_shift,\n        x_combined_with_region2,\n        y_combined_with_region2,\n        p0=p0,\n        maxfev=10000,\n    )  # (4,)\n    # Extract the optimized shift and polynomial coefficients\n    optimized_shift = popt[0]\n    #assert optimized_shift > 0.8\n    if optimized_shift <= 0.8:\n        optimized_shift = 0.8\n    optimized_poly_params = popt[1:]  # (3,)\n    return fit_function, optimized_poly_params, optimized_shift\n\ndef feature_engineering(signal_train):\n    \"\"\"Create a dataframe with two features from the raw data.\n    Parameters:\n    f_raw: ndarray of shape (n_planets, 67500)\n    a_raw: ndarray of shape (n_planets, 5625)\n    Return value:\n    df: DataFrame of shape (n_planets, 2)\n    \"\"\"\n    y_shifts = []\n    for IDX in tqdm(range(len(signal_train))):  # (1100, 375, 283)\n        try:\n            data = signal_train[IDX]  # (375, 283)\n            optimized_breakpoints = [all_bp[IDX].item(), all_bp2[IDX].item()]  # [130, 272]\n            fit_func, params, y_shift = poly_exp_fit(\n                data[:, 1:].mean(1) / data[:, 1:].mean(1).mean(),  # (375,)\n                optimized_breakpoints,\n                buffer_size_poly,\n                degree=2,\n            )\n            y_shifts.append(y_shift)\n        except Exception as error:\n            print(f\"Feature engineering error {IDX}: {error}\")\n            y_shifts.append(0.991)\n    y_shifts = np.array(y_shifts)\n    df = pd.DataFrame(1 - y_shifts, index=star_info.index,)\n    return df\n\ndef postprocessing(pred_array, index, sigma_pred):\n    \"\"\"Create a submission dataframe from its components\n    Parameters:\n    pred_array: ndarray of shape (n_samples, 283)\n    index: pandas.Index of length n_samples with name 'planet_id'\n    sigma_pred: series of length n_samples or float\n    Return value:\n    df: DataFrame of shape (n_samples, 566) with planet_id as index\n    \"\"\"\n    if isinstance(sigma_pred, float):\n        expanded_sigmas = np.ones(len(pred_array)) * sigma_pred\n    else:\n        expanded_sigmas = sigma_pred\n    expanded_sigmas = np.repeat(expanded_sigmas[:, np.newaxis], 283, axis=1)\n    if pred_array.shape[1] == 1:\n        pred_array = np.repeat(pred_array, 283, axis=1)\n    return pd.concat(\n        [\n            pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns),\n            pd.DataFrame(expanded_sigmas, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)]),\n        ],\n        axis=1,\n    )\n\n# breakpoint detection\nall_bp, all_bp2 = get_breakpoints(signal_train, PRE_BINNED_TIME, verbose=False)  # signal_train => (1100, 375, 283)\npd.DataFrame(dict(all_bp=all_bp, all_bp2=all_bp2)).to_csv('all_bp.csv', index=False)\n\ndf = feature_engineering(signal_train)\ndisplay(df)\npredictions = df.values\nprint(predictions.shape)\nstar_info = pd.read_csv(ROOT + f\"/{MODE}_star_info.csv\")\nstar_info[\"planet_id\"] = star_info[\"planet_id\"].astype(int)\nstar_info = star_info.set_index(\"planet_id\")\ndisplay(star_info)\n\nsub_df = postprocessing(predictions, star_info.index, sigma_pred=0.0008)\ndisplay(sub_df)\nsub_df.to_csv('submission.csv')\npd.read_csv(\"submission.csv\")\n! rm signal_v0.npy preprocess.py\n! ls\n","metadata":{"execution":{"iopub.status.busy":"2025-07-08T14:04:49.252868Z","iopub.execute_input":"2025-07-08T14:04:49.253158Z","iopub.status.idle":"2025-07-08T14:04:49.330098Z","shell.execute_reply.started":"2025-07-08T14:04:49.253135Z","shell.execute_reply":"2025-07-08T14:04:49.329371Z"},"trusted":true},"outputs":[],"execution_count":null}]}