{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.13"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":1087.424982,"end_time":"2024-08-03T13:00:18.164826","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-08-03T12:42:10.739844","version":"2.5.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel Data Challenge 2024: Introductory model: training\n\nIn this notebook, we show how to train and cross-validate a model. At the end, we save the model so that it can be used for inference in the separate notebook [ADC24 Intro inference](https://www.kaggle.com/code/ambrosm/adc24-intro-inference).\n","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.009657,"end_time":"2024-08-03T12:42:14.074227","exception":false,"start_time":"2024-08-03T12:42:14.06457","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import pandas as pd\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\nimport pickle\nimport itertools\nimport time\nimport os\nimport logging\nimport glob \nfrom astropy.stats import sigma_clip\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error","metadata":{"_kg_hide-input":true,"papermill":{"duration":3.023083,"end_time":"2024-08-03T12:42:17.107326","exception":false,"start_time":"2024-08-03T12:42:14.084243","status":"completed"},"tags":[],"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:29:12.187213Z","iopub.execute_input":"2025-01-14T13:29:12.187576Z","iopub.status.idle":"2025-01-14T13:29:15.874401Z","shell.execute_reply.started":"2025-01-14T13:29:12.18753Z","shell.execute_reply":"2025-01-14T13:29:15.873131Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# A look at the data\n\nWe start by reading the metadata:","metadata":{"papermill":{"duration":0.008428,"end_time":"2024-08-03T12:42:17.180297","exception":false,"start_time":"2024-08-03T12:42:17.171869","status":"completed"},"tags":[]}},{"cell_type":"code","source":"path_folder = '/kaggle/input/ariel-data-challenge-2024/' # path to the folder containing the data\npath_out = '/kaggle/tmp/data_light_raw/' # path to the folder to store the light data\noutput_dir = '/kaggle/tmp/data_light_raw/' # path for the output directory","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:30:13.1908Z","iopub.execute_input":"2025-01-14T13:30:13.191194Z","iopub.status.idle":"2025-01-14T13:30:13.196101Z","shell.execute_reply.started":"2025-01-14T13:30:13.191164Z","shell.execute_reply":"2025-01-14T13:30:13.195098Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if not os.path.exists(path_out):\n    os.makedirs(path_out)\n    print(f\"Directory {path_out} created.\")\nelse:\n    print(f\"Directory {path_out} already exists.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:30:13.919441Z","iopub.execute_input":"2025-01-14T13:30:13.919826Z","iopub.status.idle":"2025-01-14T13:30:13.926221Z","shell.execute_reply.started":"2025-01-14T13:30:13.919796Z","shell.execute_reply":"2025-01-14T13:30:13.925155Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"CHUNKS_SIZE = 1","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:30:14.774156Z","iopub.execute_input":"2025-01-14T13:30:14.775182Z","iopub.status.idle":"2025-01-14T13:30:14.779463Z","shell.execute_reply.started":"2025-01-14T13:30:14.775148Z","shell.execute_reply":"2025-01-14T13:30:14.778358Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def adc_convert(signal, gain, offset):\n    \"\"\"\n    Reverts the ADC process by applying the gain and offset to the signal.\n\n    Parameters:\n    - signal (np.ndarray): The digital signal to be converted back to analog form.\n    - gain (float or np.ndarray): The gain factor from the calibration data.\n    - offset (float or np.ndarray): The offset value from the calibration data.\n\n    Returns:\n    - np.ndarray: The converted analog signal.\n    \"\"\"\n    signal = signal.astype(np.float64)  # Ensure precision\n    signal /= gain\n    signal += offset   # Apply offset correction and scale by gain\n    return signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:41.625534Z","iopub.execute_input":"2025-01-14T13:31:41.626187Z","iopub.status.idle":"2025-01-14T13:31:41.632554Z","shell.execute_reply.started":"2025-01-14T13:31:41.626149Z","shell.execute_reply":"2025-01-14T13:31:41.631118Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from scipy.ndimage import generic_filter\ndef mask_and_correct_bad_pixels(signal, dead, dark_frame):\n    \"\"\"\n    Masks and corrects hot and dead pixels in the input signal using provided maps and dark frame data.\n\n    Parameters:\n    - signal (np.ndarray): The input signal array (frames, height, width).\n    - dead_pixel_map (np.ndarray): Binary map of dead pixels (height, width).\n    - dark_frame (np.ndarray): Dark frame used to identify hot pixels (height, width).\n\n    Returns:\n    - np.ndarray: Signal with hot and dead pixels corrected.\n    \"\"\"\n    # Identify hot pixels using sigma clipping on the dark frame\n    hot = sigma_clip(dark_frame, sigma=5, maxiters=5).mask\n\n    # # Convert hot_pixel_mask to boolean type\n    # hot_pixel_mask = hot_pixel_mask.astype(bool)\n\n    # # Ensure dead_pixel_map is boolean (if it's not already)\n    # dead_pixel_mask = dead_pixel_map.astype(bool)\n\n    # Broadcast masks to match the dimensions of the signal\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal)\n\n    return signal\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:42.004402Z","iopub.execute_input":"2025-01-14T13:31:42.004769Z","iopub.status.idle":"2025-01-14T13:31:42.011849Z","shell.execute_reply.started":"2025-01-14T13:31:42.004744Z","shell.execute_reply":"2025-01-14T13:31:42.010512Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from scipy.interpolate import CubicSpline\n\ndef apply_linear_corr_spline(linear_corr, clean_signal):\n    \"\"\"\n    Applies non-linearity correction using spline interpolation.\n\n    Parameters:\n    - linear_corr (np.ndarray): Spline control points or correction coefficients \n                                (shape: degree + 1, height, width).\n    - clean_signal (np.ndarray): Signal array to correct (shape: frames, height, width).\n\n    Returns:\n    - np.ndarray: Corrected signal array with non-linearity addressed using spline interpolation.\n    \"\"\"\n    # Create an empty array to store the corrected signal\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(\n                range(clean_signal.shape[1]), range(clean_signal.shape[2])\n            ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:44.083489Z","iopub.execute_input":"2025-01-14T13:31:44.084622Z","iopub.status.idle":"2025-01-14T13:31:44.09075Z","shell.execute_reply.started":"2025-01-14T13:31:44.084588Z","shell.execute_reply":"2025-01-14T13:31:44.089524Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef clean_dark(signal, dead, dark, dt):\n    \"\"\"\n    Subtracts dark current from the signal while accounting for dead pixels.\n\n    Parameters:\n    - signal (np.ndarray): Observed signal array (shape: frames, height, width).\n    - dead (np.ndarray): Dead pixel mask (boolean array, shape: height, width).\n    - dark (np.ndarray): Dark current map (shape: height, width).\n    - dt (np.ndarray): Integration time array (shape: frames).\n\n    Returns:\n    - np.ndarray: Signal array with dark current subtracted.\n    \"\"\"\n    # Mask dark current values at dead pixel positions\n    masked_dark = np.ma.masked_where(dead, dark)\n    \n    # Expand dark map to match the shape of the signal array\n    expanded_dark = np.tile(masked_dark, (signal.shape[0], 1, 1))\n\n    # Ensure integration time array is broadcastable (frames, 1, 1)\n    expanded_dt = dt[:, np.newaxis, np.newaxis]\n\n    # Subtract the scaled dark current from the signal\n    corrected_signal = signal - expanded_dark * expanded_dt\n\n    return corrected_signal\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:44.545572Z","iopub.execute_input":"2025-01-14T13:31:44.545939Z","iopub.status.idle":"2025-01-14T13:31:44.552225Z","shell.execute_reply.started":"2025-01-14T13:31:44.545914Z","shell.execute_reply":"2025-01-14T13:31:44.551114Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_cds(signal):\n    \"\"\"\n    Computes the Correlated Double Sampling (CDS) by subtracting the start-of-exposure\n    frames from the end-of-exposure frames.\n\n    Parameters:\n    - signal (np.ndarray): Signal array (shape: frames, height, width), where frames\n                           alternate between the start and end of exposure.\n\n    Returns:\n    - np.ndarray: CDS array (shape: frames/2, height, width).\n    \"\"\"\n\n    \n    cds = signal[:, 1::2, :, :] - signal[:, ::2, :, :]\n\n    return cds\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:47.767907Z","iopub.execute_input":"2025-01-14T13:31:47.768318Z","iopub.status.idle":"2025-01-14T13:31:47.774308Z","shell.execute_reply.started":"2025-01-14T13:31:47.768288Z","shell.execute_reply":"2025-01-14T13:31:47.77304Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def bin_obs(cds_signal, binning):\n    \"\"\"\n    Bins the time series data at the specified frequency (binning).\n\n    Parameters:\n    - cds_signal (np.ndarray): Input CDS signal array (shape: frames, height, width)\n    - binning (int): The binning factor (number of frames to bin together)\n\n    Returns:\n    - np.ndarray: Binned CDS signal array (shape: frames//binning, height, width)\n    \"\"\"\n\n    # Transpose the signal to facilitate the binning along the time axis (frames)\n    cds_transposed = cds_signal.transpose(0, 1, 3, 2)\n\n    # Compute the new shape for the binned signal\n    new_frames = cds_transposed.shape[1] // binning\n    cds_binned = np.zeros((cds_transposed.shape[0], new_frames, cds_transposed.shape[2], cds_transposed.shape[3]))\n\n    # Perform binning (sum the frames within each bin)\n    for i in range(new_frames):\n        cds_binned[:, i, :, :] = np.sum(cds_transposed[:, i * binning:(i + 1) * binning, :, :], axis=1)\n\n    return cds_binned\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:31:50.765848Z","iopub.execute_input":"2025-01-14T13:31:50.76689Z","iopub.status.idle":"2025-01-14T13:31:50.773438Z","shell.execute_reply.started":"2025-01-14T13:31:50.766839Z","shell.execute_reply":"2025-01-14T13:31:50.772113Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef correct_flat_field(flat, dead, signal):\n    \"\"\"\n    Corrects the signal using flat field data by accounting for pixel-to-pixel variations in the detector.\n\n    Parameters:\n    - flat (np.ndarray): Flat field correction map (shape: height, width)\n    - dead (np.ndarray): Dead pixel map (shape: height, width)\n    - signal (np.ndarray): Input signal to correct (shape: frames, height, width)\n\n    Returns:\n    - np.ndarray: Corrected signal array (shape: frames, height, width)\n    \"\"\"\n    # Ensure the flat field and dead pixel map are in the correct orientation (height, width)\n    flat = flat.transpose(1, 0) # Transpose to align with the expected signal shape\n    dead = dead.transpose(1, 0)\n    \n    # Mask dead pixels in the flat field to avoid dividing by zero or erroneous values\n    flat = np.ma.masked_where(dead, flat)\n\n    # Expand flat field to match the signal's frame count\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    # Avoid division by zero by masking flat field pixels that are zero or NaN\n    with np.errstate(divide='ignore', invalid='ignore'):\n        corrected_signal = signal / flat\n\n    # Mask the output signal where the flat field was masked (dead pixels)\n    corrected_signal = np.ma.masked_where(np.ma.getmask(flat), corrected_signal)\n\n    return corrected_signal\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:32:04.692441Z","iopub.execute_input":"2025-01-14T13:32:04.692805Z","iopub.status.idle":"2025-01-14T13:32:04.699701Z","shell.execute_reply.started":"2025-01-14T13:32:04.692777Z","shell.execute_reply":"2025-01-14T13:32:04.698533Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"## we will start by getting the index of the training data:\ndef get_index(files,CHUNKS_SIZE ):\n    index = []\n    for file in files :\n        file_name = file.split('/')[-1]\n        if file_name.split('_')[0] == 'AIRS-CH0' and file_name.split('_')[1] == 'signal.parquet':\n            file_index = os.path.basename(os.path.dirname(file))\n            index.append(int(file_index))\n    index = np.array(index)\n    index = np.sort(index) \n\n    index=np.array_split(index, len(index)//CHUNKS_SIZE)\n    \n    return index","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:32:11.250498Z","iopub.execute_input":"2025-01-14T13:32:11.250872Z","iopub.status.idle":"2025-01-14T13:32:11.257583Z","shell.execute_reply.started":"2025-01-14T13:32:11.250843Z","shell.execute_reply":"2025-01-14T13:32:11.256387Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"files = glob.glob(os.path.join(path_folder + 'train/', '*/*'))\nlen(files)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T13:54:42.33767Z","iopub.execute_input":"2025-01-12T13:54:42.338045Z","iopub.status.idle":"2025-01-12T13:54:42.950162Z","shell.execute_reply.started":"2025-01-12T13:54:42.338019Z","shell.execute_reply":"2025-01-12T13:54:42.948981Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def calibrate_training_data(files, CHUNKS_SIZE, path_folder, path_out, train_adc_info, axis_info, DO_MASK=True, DO_THE_NL_CORR=True, DO_DARK=True, DO_FLAT=True, TIME_BINNING=True, cut_inf=38, cut_sup=321):\n\n    cut_inf = 38\n    cut_sup = 321\n    \n    # Load index and necessary files\n    # glob.glob gives the filenames of all the paths that match the location given\n    files = glob.glob(os.path.join(path_folder + 'train/', '*/*'))\n    index = get_index(files[:], CHUNKS_SIZE)  # 22 is used as an example, adjust as needed\n\n    # Reading ADC info and axis info\n    train_adc_info = pd.read_csv(os.path.join(path_folder, 'train_adc_info.csv')).set_index('planet_id')\n    axis_info = pd.read_parquet(os.path.join(path_folder, 'axis_info.parquet'))\n    length = len(index)\n    avg_time = 0\n    time_vals = []\n    # Processing each chunk of files\n    for n, index_chunk in enumerate(index):\n        num_bars = int(n/length * 20)\n        print(f'{n}/{length} [{\"#\" * num_bars + \" \" * (20-num_bars)}] {np.round((length - n) * avg_time, 2)}s  {np.round(avg_time, 2)}s/it')\n        start = time.time()\n        for i in range(CHUNKS_SIZE):\n            # print(n, index_chunk, i)\n            AIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, cut_sup - cut_inf))\n            FGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))            \n\n            # Read and preprocess AIRS signal\n            df = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_signal.parquet'))\n            signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 356))\n            gain = train_adc_info['AIRS-CH0_adc_gain'].loc[index_chunk[i]]\n            offset = train_adc_info['AIRS-CH0_adc_offset'].loc[index_chunk[i]]\n            signal = adc_convert(signal, gain, offset)\n            dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n            dt_airs[1::2] += 0.1\n            chopped_signal = signal[:, :, cut_inf:cut_sup]\n            del df, signal, gain, offset\n\n            # Read calibration files for AIRS\n            flat_airs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n            dark_airs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n            dead_airs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n            linear_corr = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n            \n            # Mask and correct hot/dead pixels for AIRS\n            if DO_MASK:\n                chopped_signal = mask_and_correct_bad_pixels(chopped_signal, dead_airs, dark_airs)\n                AIRS_CH0_clean[i] = chopped_signal\n                del chopped_signal\n\n            # Apply non-linearity correction for AIRS\n            if DO_THE_NL_CORR:\n                linear_corr_signal = apply_linear_corr_spline(linear_corr, AIRS_CH0_clean[i])\n                AIRS_CH0_clean[i] = linear_corr_signal\n                del linear_corr_signal\n            del linear_corr\n\n            # Apply dark correction for AIRS\n            if DO_DARK:\n                cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark_airs, dt_airs)\n                AIRS_CH0_clean[i] = cleaned_signal\n                del cleaned_signal\n            del dark_airs, dt_airs\n\n            # Process FGS1 signal\n            df_fgs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_signal.parquet'))\n            fgs_signal = df_fgs.values.astype(np.float64).reshape((df_fgs.shape[0], 32, 32))\n            FGS1_gain = train_adc_info['FGS1_adc_gain'].loc[index_chunk[i]]\n            FGS1_offset = train_adc_info['FGS1_adc_offset'].loc[index_chunk[i]]\n            fgs_signal = adc_convert(fgs_signal, FGS1_gain, FGS1_offset)\n            dt_fgs1 = np.ones(len(fgs_signal)) * 0.1\n            dt_fgs1[1::2] += 0.1\n            chopped_FGS1 = fgs_signal\n            del df_fgs, fgs_signal, FGS1_gain, FGS1_offset\n\n            # Read calibration files for FGS1\n            flat_fgs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n            dark_fgs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 32))\n            dead_fgs1 = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 32))\n            linear_corr_fgs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 32))\n\n            # Mask and correct hot/dead pixels for FGS1\n            if DO_MASK:\n                chopped_FGS1 = mask_and_correct_bad_pixels(chopped_FGS1, dead_fgs1, dark_fgs)\n            FGS1_clean[i] = chopped_FGS1\n            del chopped_FGS1\n            # Apply non-linearity correction for FGS1\n            if DO_THE_NL_CORR:\n                linear_corr_signal = apply_linear_corr_spline(linear_corr_fgs, FGS1_clean[i])\n                FGS1_clean[i] = linear_corr_signal\n                del linear_corr_signal\n            del linear_corr_fgs\n\n            # Apply dark correction for FGS1\n            if DO_DARK:\n                cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark_fgs, dt_fgs1)\n                FGS1_clean[i] = cleaned_signal\n                del cleaned_signal\n            del dark_fgs, dt_fgs1\n\n        # Save data and perform time binning\n        \n        if AIRS_CH0_clean.shape[0] == 1:\n            AIRS_cds = get_cds(AIRS_CH0_clean)\n        FGS1_cds = get_cds(FGS1_clean)\n        del AIRS_CH0_clean, FGS1_clean\n                \n        if TIME_BINNING:\n            AIRS_cds_binned = bin_obs(AIRS_cds, binning=30)\n            FGS1_cds_binned = bin_obs(FGS1_cds, binning=30 * 12)\n        else:\n            AIRS_cds = AIRS_cds.transpose(0, 1, 3, 2)\n            AIRS_cds_binned = AIRS_cds\n            FGS1_cds = FGS1_cds\n\n        del AIRS_cds, FGS1_cds\n        \n        # Apply flat field correction\n        for i in range(CHUNKS_SIZE):\n            flat_airs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n            flat_fgs = pd.read_parquet(os.path.join(path_folder, f'train/{index_chunk[i]}/FGS1_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n            if DO_FLAT:\n                corrected_AIRS_cds_binned = correct_flat_field(flat_airs, dead_airs, AIRS_cds_binned[i])\n                AIRS_cds_binned[i] = corrected_AIRS_cds_binned\n                corrected_FGS1_cds_binned = correct_flat_field(flat_fgs, dead_fgs1, FGS1_cds_binned[i])\n                FGS1_cds_binned[i] = corrected_FGS1_cds_binned\n                del corrected_AIRS_cds_binned, corrected_FGS1_cds_binned\n                \n            del flat_airs, flat_fgs\n\n        # Save the processed data\n        np.save(os.path.join(path_out, f'AIRS_clean_train_{index_chunk[i]}.npy'), AIRS_cds_binned)\n        np.save(os.path.join(path_out, f'FGS1_train_{index_chunk[i]}.npy'), FGS1_cds_binned)\n        del AIRS_cds_binned, FGS1_cds_binned\n        time_vals.append(time.time() - start)\n        avg_time = sum(time_vals)/len(time_vals)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-14T13:34:09.779408Z","iopub.execute_input":"2025-01-14T13:34:09.779805Z","iopub.status.idle":"2025-01-14T13:34:09.803739Z","shell.execute_reply.started":"2025-01-14T13:34:09.779767Z","shell.execute_reply":"2025-01-14T13:34:09.802659Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"files = glob.glob(os.path.join(path_folder + 'train/', '*/*'))\n\n# index = get_index(files[:22],CHUNKS_SIZE) \n\ntrain_adc_info = pd.read_csv(os.path.join(path_folder, 'train_adc_info.csv'))\ntrain_adc_info = train_adc_info.set_index('planet_id')\naxis_info = pd.read_parquet(os.path.join(path_folder,'axis_info.parquet'))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T13:55:00.29919Z","iopub.execute_input":"2025-01-12T13:55:00.300333Z","iopub.status.idle":"2025-01-12T13:55:00.710509Z","shell.execute_reply.started":"2025-01-12T13:55:00.300293Z","shell.execute_reply":"2025-01-12T13:55:00.709468Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"calibrate_training_data(files, CHUNKS_SIZE, path_folder, path_out, train_adc_info, axis_info, \n                             DO_MASK=True, DO_THE_NL_CORR=True, DO_DARK=True, DO_FLAT=True, TIME_BINNING=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T13:55:22.095582Z","iopub.execute_input":"2025-01-12T13:55:22.096019Z","iopub.status.idle":"2025-01-12T13:55:24.323748Z","shell.execute_reply.started":"2025-01-12T13:55:22.095984Z","shell.execute_reply":"2025-01-12T13:55:24.319848Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# def load_data (file, chunk_size, nb_files): \n#     data0 = np.load(file + '_0.npy')\n#     data_all = np.zeros((nb_files*chunk_size, data0.shape[1], data0.shape[2], data0.shape[3]))\n#     data_all[:chunk_size] = data0\n#     for i in range (1, nb_files): \n#         data_all[i*chunk_size:(i+1)*chunk_size] = np.load(file + '_{}.npy'.format(i))\n#     return data_all \n\n# data_train = load_data(path_out + 'AIRS_clean_train', CHUNKS_SIZE, len(index)) \n# data_train_FGS = load_data(path_out + 'FGS1_train', CHUNKS_SIZE, len(index))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:39:37.489416Z","iopub.execute_input":"2025-01-12T15:39:37.490201Z","iopub.status.idle":"2025-01-12T15:39:37.741933Z","shell.execute_reply.started":"2025-01-12T15:39:37.490168Z","shell.execute_reply":"2025-01-12T15:39:37.740343Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# np.save('/kaggle/working/' + 'data_train.npy', data_train)\n# np.save('/kaggle/working/' + 'data_train_FGS.npy', data_train_FGS)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:39:55.455238Z","iopub.execute_input":"2025-01-12T15:39:55.455629Z","iopub.status.idle":"2025-01-12T15:39:55.50091Z","shell.execute_reply.started":"2025-01-12T15:39:55.455603Z","shell.execute_reply":"2025-01-12T15:39:55.499523Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# data_train=np.load('/kaggle/working/data_train.npy')\n# data_train_FGS=np.load('/kaggle/working/data_train_FGS.npy')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T19:18:25.519787Z","iopub.execute_input":"2025-01-12T19:18:25.520157Z","iopub.status.idle":"2025-01-12T19:18:25.950054Z","shell.execute_reply.started":"2025-01-12T19:18:25.520132Z","shell.execute_reply":"2025-01-12T19:18:25.948689Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import matplotlib.pyplot as plt \n\n# print('Shape of the training datasset: \\t')\n# print('\\n For AIRS-CH0:', data_train.shape)\n# print('\\n For FGS1:', data_train_FGS.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:43:00.56038Z","iopub.execute_input":"2025-01-12T15:43:00.560781Z","iopub.status.idle":"2025-01-12T15:43:00.566916Z","shell.execute_reply.started":"2025-01-12T15:43:00.560749Z","shell.execute_reply":"2025-01-12T15:43:00.565681Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# sns.heatmap(data_train[-1,5,:,:].T)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:46:08.24661Z","iopub.execute_input":"2025-01-12T15:46:08.247018Z","iopub.status.idle":"2025-01-12T15:46:08.795595Z","shell.execute_reply.started":"2025-01-12T15:46:08.246989Z","shell.execute_reply":"2025-01-12T15:46:08.794362Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# sns.heatmap(data_train_FGS[-1,2,:,:].T)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:52:16.243206Z","iopub.execute_input":"2025-01-12T12:52:16.2436Z","iopub.status.idle":"2025-01-12T12:52:16.692961Z","shell.execute_reply.started":"2025-01-12T12:52:16.243574Z","shell.execute_reply":"2025-01-12T12:52:16.691914Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for i in range(len(data_train_FGS)) : \n#     light_curve = data_train_FGS[i,:,:,:].sum(axis=(1,2))\n#     plt.plot(light_curve/light_curve.mean(), '-', alpha=0.3)\n# plt.xlabel('Time (frame index)')\n# plt.ylabel('Normalized flux in the frame')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T16:13:33.860825Z","iopub.execute_input":"2025-01-12T16:13:33.86122Z","iopub.status.idle":"2025-01-12T16:13:34.103834Z","shell.execute_reply.started":"2025-01-12T16:13:33.86119Z","shell.execute_reply":"2025-01-12T16:13:34.102343Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for i in range(len(data_train)) : \n#     light_curve = data_train[i,:,:,:].sum(axis=(1,2))\n#     plt.plot(light_curve/light_curve.mean(), '-', alpha=0.3)\n# plt.xlabel('Time (frame index)')\n# plt.ylabel('Normalized flux in the frame')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:45:01.952913Z","iopub.execute_input":"2025-01-12T15:45:01.95334Z","iopub.status.idle":"2025-01-12T15:45:02.236243Z","shell.execute_reply.started":"2025-01-12T15:45:01.953308Z","shell.execute_reply":"2025-01-12T15:45:02.234941Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# train_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv',\n#                            index_col='planet_id')\n# train_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',\n#                            index_col='planet_id')\n# wavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\n# axis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n","metadata":{"papermill":{"duration":0.188106,"end_time":"2024-08-03T12:42:17.37713","exception":false,"start_time":"2024-08-03T12:42:17.189024","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2025-01-12T15:46:43.224651Z","iopub.execute_input":"2025-01-12T15:46:43.225068Z","iopub.status.idle":"2025-01-12T15:46:43.396622Z","shell.execute_reply.started":"2025-01-12T15:46:43.225036Z","shell.execute_reply":"2025-01-12T15:46:43.395241Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Some facts about the data:\n- We have 673 planets for training. These planets belong to two different stars.\n- There will be roughly 800 planets for testing (but the test data is hidden).\n- The competition is a multi-output regression task with 283 targets to predict.","metadata":{"papermill":{"duration":0.008486,"end_time":"2024-08-03T12:42:17.394634","exception":false,"start_time":"2024-08-03T12:42:17.386148","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"## The FGS1 data\n\nHaving read the metadata, we'll tackle the FGS1 data (Fine Guidance System). The FGS1 measurements consist of one file per planet (673 files for 673 planets for training). For now, we ignore the calibration files.\n\nEach file contains 135,000 rows of images taken at 0.1 second time steps. Each row is a 32\\*32 image at a single wavelength.\n\nWe read a sample file:","metadata":{"papermill":{"duration":0.008446,"end_time":"2024-08-03T12:42:17.411955","exception":false,"start_time":"2024-08-03T12:42:17.403509","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"Every row of the file corresponds to an image of a star. The images come in pairs, and the second image is lighter than the first one:","metadata":{"papermill":{"duration":0.009711,"end_time":"2024-08-03T12:42:19.199727","exception":false,"start_time":"2024-08-03T12:42:19.190016","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# _, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n# sns.heatmap(f_signal.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000)\n# ax1.set_aspect('equal')\n# sns.heatmap(f_signal.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000)\n# ax2.set_aspect('equal')\n# plt.suptitle('A pair of FGS1 images')\n# plt.show()","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.559852,"end_time":"2024-08-03T12:42:19.769794","exception":false,"start_time":"2024-08-03T12:42:19.209942","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2025-01-12T16:13:53.85083Z","iopub.execute_input":"2025-01-12T16:13:53.851225Z","iopub.status.idle":"2025-01-12T16:13:54.203439Z","shell.execute_reply.started":"2025-01-12T16:13:53.851196Z","shell.execute_reply":"2025-01-12T16:13:54.201804Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"To see the time series, we first have to compute the difference between the even and the odd frames to get the net signal (67500 time steps). We then take the mean over all 1024 pixels. The net signal is very noisy, and we smoothen it by computing a moving average. The plot of the smoothened signal clearly shows that the signal intensity is reduced (i.e., the image gets darker) while the planet passes in front of the star (between time steps 23500 and 44000).\n\nThe left diagram shows a planet with a strong reduction of the signal intensity, the right diagram shows a planet with a weak reduction:","metadata":{"papermill":{"duration":0.011386,"end_time":"2024-08-03T12:42:19.791847","exception":false,"start_time":"2024-08-03T12:42:19.780461","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# planet_id = 499191466\n# f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/{planet_id}/AIRS-CH0_signal.parquet')\n\n# # net_signal = mean_signal[1::2] - mean_signal[0::2]\n# # window=800\n# # smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n# # plt.plot(net_signal, label='raw signal')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:53:24.849139Z","iopub.execute_input":"2025-01-12T15:53:24.85069Z","iopub.status.idle":"2025-01-12T15:53:28.180648Z","shell.execute_reply.started":"2025-01-12T15:53:24.850649Z","shell.execute_reply":"2025-01-12T15:53:28.179559Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# mean_signal = f_signal.mean(axis=1)\n# net_signal = mean_signal[1::2] - mean_signal[0::2]\n# cum_signal = net_signal.cumsum()\n# window=800\n# smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T15:50:52.952913Z","iopub.execute_input":"2025-01-12T15:50:52.953327Z","iopub.status.idle":"2025-01-12T15:50:53.116676Z","shell.execute_reply.started":"2025-01-12T15:50:52.953293Z","shell.execute_reply":"2025-01-12T15:50:53.115547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# _, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=(12, 4))\n\n# # planet_id = 14485303\n# # f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\n# mean_signal = f_signal.values.mean(axis=1)\n# net_signal = mean_signal[1::2] - mean_signal[0::2]\n# cum_signal = net_signal.cumsum()\n# window=800\n# smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n# ax1.set_title('FGS1: time series of planet with strong signal')\n# ax1.plot(net_signal, label='raw signal')\n# ax1.legend()\n# ax3.plot(smooth_signal, color='c', label='smoothened signal')\n# ax3.legend()\n# ax3.set_xlabel('time step')\n# for time_step in [20500, 23500, 44000, 47000]:\n#     ax3.axvline(time_step, color='gray')\n\n# # planet_id = 4249337798\n# # f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\n# mean_signal = f_signal..mean(axis=1)\n# net_signal = mean_signal[1::2] - mean_signal[0::2]\n# cum_signal = net_signal.cumsum()\n# window=800\n# smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n# ax2.set_title('FGS1: time series of planet with weak signal')\n# ax2.plot(net_signal, label='raw signal')\n# ax2.legend()\n# ax4.plot(smooth_signal, color='c', label='smoothened signal')\n# ax4.legend()\n# ax4.set_xlabel('time step')\n# for time_step in [20500, 23500, 44000, 47000]:\n#     ax4.axvline(time_step, color='gray')\n\n# # plt.suptitle('FGS1 time series', y=0.96)\n# plt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-01-12T15:54:35.602681Z","iopub.execute_input":"2025-01-12T15:54:35.603142Z","iopub.status.idle":"2025-01-12T15:54:37.003179Z","shell.execute_reply.started":"2025-01-12T15:54:35.60311Z","shell.execute_reply":"2025-01-12T15:54:37.002017Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import pandas as pd\n# import numpy as np\n# import matplotlib.pyplot as plt\n# from scipy.ndimage import gaussian_filter1d\n\n# # Example planet IDs\n# planet_ids = [14485303, 4249337798]\n\n# def load_calibrated_signal(planet_id):\n#     f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n#     # dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_calibration/dark.parquet')\n#     # dead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_calibration/dead.parquet')\n#     # flat = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_calibration/flat.parquet')\n#     # linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_calibration/linear_corr.parquet')\n\n#     # print(f\"Planet ID: {planet_id}\")\n#     # print(f\"f_signal shape: {f_signal.shape}, columns: {f_signal.columns}\")\n#     # print(f\"dark shape: {dark.shape}\")\n#     # print(f\"dead shape: {dead.shape}, columns: {dead.columns}\")\n#     # print(f\"flat shape: {flat.shape}\")\n#     # print(f\"linear_corr columns: {linear_corr.columns}\")\n\n#     # # Dark correction: Subtract dark signal from raw signal\n#     # print(\"Applying dark correction...\")\n#     # f_signal -= dark.reindex(f_signal.index, method='nearest')\n\n#     # # Dead pixel correction: Mark dead pixels as NaN\n#     # print(\"Applying dead pixel correction...\")\n#     # dead = dead.reindex(f_signal.index)\n#     # for col in dead.columns:\n#     #     dead_pixels = dead[col].to_numpy()\n#     #     if dead_pixels.shape[0] == f_signal.shape[0]:\n#     #         f_signal.loc[dead_pixels == 1, col] = np.nan\n#     # print(f\"NaN count after dead pixel correction: {f_signal.isna().sum().sum()}\")\n\n#     # # Interpolation: Fill NaN values through interpolation\n#     # f_signal = f_signal.interpolate(axis=0).ffill().bfill()\n#     # print(f\"NaN count after interpolation: {f_signal.isna().sum().sum()}\")\n\n#     # # Flat field correction: Divide by the flat field (and handle zero flat field values)\n#     # print(\"Applying flat field correction...\")\n#     # flat[flat == 0] = np.nan  # Replace zero values with NaN\n#     # f_signal = f_signal.divide(flat).fillna(0)  # Use fillna(0) to avoid issues with NaN values in flat field\n\n#     # # Linearity correction: Apply polynomial coefficients for linearity correction\n#     # print(\"Applying linearity correction...\")\n#     # for col in f_signal.columns:\n#     #     if col in linear_corr.columns:\n#     #         coefficients = linear_corr[col].values[::-1]  # Reverse coefficients for polynomial evaluation\n#     #         coefficients = np.clip(coefficients, -1e6, 1e6)  # Clip the coefficients to avoid overflow\n#     #         f_signal[col] = np.polyval(coefficients, f_signal[col].values)\n\n#     # # Ensure no NaN or Inf values after calibration\n#     # f_signal = f_signal.replace([np.inf, -np.inf], np.nan).fillna(0)\n\n#     # print(\"Calibration complete.\")\n#     # print(f\"Calibrated signal preview:\\n{f_signal.mean(axis=1).head()}\")\n#     return f_signal\n\n\n# def process_signal(f_signal):\n#     print(\"Processing signal...\")\n#     mean_signal = f_signal.mean(axis=1)\n#     print(f\"Mean signal sample: {mean_signal[:5]}\")\n\n#     # Compute the net signal as the difference between consecutive frames\n#     net_signal = mean_signal[1::2] - mean_signal[0::2]\n#     print(f\"Net signal sample: {net_signal[:5]}\")\n\n#     # Smooth the signal using a Gaussian filter\n#     smoothed_signal = gaussian_filter1d(net_signal, sigma=20)\n#     print(f\"Smoothed signal sample: {smoothed_signal[:5]}\")\n#     return net_signal, smoothed_signal\n\n# def extract_features(smoothed_signal):\n#     print(\"Extracting features...\")\n#     transit_start = np.argmin(smoothed_signal[20000:30000]) + 20000\n#     transit_end = np.argmin(smoothed_signal[40000:50000]) + 40000\n#     transit_depth = smoothed_signal[transit_start:transit_end].min()\n#     transit_duration = transit_end - transit_start\n\n#     print(f\"Transit start: {transit_start}, Transit end: {transit_end}\")\n#     print(f\"Transit depth: {transit_depth}, Transit duration: {transit_duration}\")\n#     return transit_depth, transit_duration, transit_start, transit_end\n\n\n\n# for planet_id in planet_ids:\n#     f_signal = load_calibrated_signal(planet_id)\n#     net_signal, smoothed_signal = process_signal(f_signal)\n\n\n#     transit_depth, transit_duration, transit_start, transit_end = extract_features(smoothed_signal)\n    \n\n#     ax_raw, ax_smooth = axes[i]\n    \n#     ax_raw.plot(net_signal, label=f'Planet {planet_id} Raw Signal')\n#     ax_raw.set_title(f'Raw Signal for Planet {planet_id}')\n#     ax_raw.legend()\n\n#     ax_smooth.plot(smoothed_signal, color='c', label='Smoothed Signal')\n#     ax_smooth.axvline(transit_start, color='orange', linestyle='--', label='Transit Start')\n#     ax_smooth.axvline(transit_end, color='red', linestyle='--', label='Transit End')\n#     ax_smooth.set_title(f'Smoothed Signal for Planet {planet_id}')\n#     ax_smooth.legend()\n\n#     # Annotating features\n#     ax_smooth.annotate(f'Transit Depth: {transit_depth:.2f}', xy=(transit_start, transit_depth), \n#                        xytext=(transit_start + 5000, transit_depth + 100),\n#                        arrowprops=dict(arrowstyle='->', color='black'))\n\n# plt.tight_layout()\n# plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T10:58:18.175867Z","iopub.execute_input":"2024-12-01T10:58:18.176338Z","iopub.status.idle":"2024-12-01T10:59:10.365771Z","shell.execute_reply.started":"2024-12-01T10:58:18.176302Z","shell.execute_reply":"2024-12-01T10:59:10.364346Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## The AIRS data\n\nAIRS is the other sensor of the satellite. It produces one file per planet as well. Each file contains 11,250 rows of images captured at constant time steps. Each 32 x 356 image has been flattened into 11392 columns.","metadata":{"papermill":{"duration":0.077595,"end_time":"2024-08-03T13:00:12.362301","exception":false,"start_time":"2024-08-03T13:00:12.284706","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# planet_id = 4249337798\n# a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n# a_signal","metadata":{"papermill":{"duration":2.313418,"end_time":"2024-08-03T13:00:14.755198","exception":false,"start_time":"2024-08-03T13:00:12.44178","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:20.139045Z","iopub.execute_input":"2024-08-04T11:55:20.140139Z","iopub.status.idle":"2024-08-04T11:55:22.915138Z","shell.execute_reply.started":"2024-08-04T11:55:20.140098Z","shell.execute_reply":"2024-08-04T11:55:22.913964Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# a_signal = a_signal.values.reshape(11250, 32, 356)\n\n# plt.figure(figsize=(10, 3))\n# sns.heatmap(a_signal[1])\n# plt.ylabel('spatial dimension')\n# plt.xlabel('wavelength dimension')\n# plt.show()","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.807236,"end_time":"2024-08-03T13:00:15.817499","exception":false,"start_time":"2024-08-03T13:00:15.010263","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:22.917083Z","iopub.execute_input":"2024-08-04T11:55:22.917443Z","iopub.status.idle":"2024-08-04T11:55:23.523336Z","shell.execute_reply.started":"2024-08-04T11:55:22.917413Z","shell.execute_reply":"2024-08-04T11:55:23.522136Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The data again is a time series, and we can see how the star is obscured while the planet is passing in front of it.","metadata":{"papermill":{"duration":0.08019,"end_time":"2024-08-03T13:00:15.979865","exception":false,"start_time":"2024-08-03T13:00:15.899675","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# mean_signal = a_signal.mean(axis=2).mean(axis=1)\n# net_signal = mean_signal[1::2] - mean_signal[0::2]\n# cum_signal = net_signal.cumsum()\n# window=80\n# smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n# _, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\n# ax1.plot(net_signal, label='raw net signal')\n# ax1.legend()\n# ax2.plot(smooth_signal, color='c', label='smoothened net signal')\n# ax2.legend()\n# ax2.set_xlabel('time')\n# for time_step in [20500, 23500, 44000, 47000]:\n#     ax2.axvline(time_step * 11250 // 135000, color='gray')\n# plt.suptitle('AIRS-CH0 time series', y=0.96)\n# plt.show()\n","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.727242,"end_time":"2024-08-03T13:00:16.788841","exception":false,"start_time":"2024-08-03T13:00:16.061599","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:23.524891Z","iopub.execute_input":"2024-08-04T11:55:23.525274Z","iopub.status.idle":"2024-08-04T11:55:24.131163Z","shell.execute_reply.started":"2024-08-04T11:55:23.525234Z","shell.execute_reply":"2024-08-04T11:55:24.130046Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Reading the data\n\nWe now read the FGS1 data and the AIRS-CH0 data for all 673 training planets. As the dataset doesn't fit into RAM completely, we keep only two one-dimensional time series for every planet. At the end we'll have\n1. A time series with 67500 steps per planet taken from the FGS1 data, and\n2. A time series with 5625 steps per planet taken from the AIRS-CH0 data.\n\nWe use the Jupyter `%%writefile` cell magic to save the function code to a file. This ensures that the inference notebook will process the test data in exactly the same way as this notebook processes the training data.","metadata":{"papermill":{"duration":0.011189,"end_time":"2024-08-03T12:42:20.928036","exception":false,"start_time":"2024-08-03T12:42:20.916847","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# %%writefile f_read_and_preprocess.py\n\n# def f_read_and_preprocess(dataset, adc_info, planet_ids):\n#     \"\"\"Read the FGS1 files for all planet_ids and extract the time series.\n    \n#     Parameters\n#     dataset: 'train' or 'test'\n#     adc_info: metadata dataframe, either train_adc_info or test_adc_info\n#     planet_ids: list of planet ids\n    \n#     Returns\n#     dataframe with one row per planet_id and 67500 values per row\n    \n#     \"\"\"\n#     f_raw_train = np.full((len(planet_ids), 67500), np.nan, dtype=np.float32)\n#     for i, planet_id in tqdm(list(enumerate(planet_ids))):\n#         f_signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_signal.parquet')\n#         mean_signal = f_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / 1024 # mean over the 32*32 pixels\n#         net_signal = mean_signal[1::2] - mean_signal[0::2]\n#         f_raw_train[i] = net_signal\n#     return f_raw_train\n    ","metadata":{"papermill":{"duration":1067.258124,"end_time":"2024-08-03T13:00:08.197814","exception":false,"start_time":"2024-08-03T12:42:20.93969","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T12:27:23.283354Z","iopub.execute_input":"2024-08-04T12:27:23.283949Z","iopub.status.idle":"2024-08-04T12:27:23.294819Z","shell.execute_reply.started":"2024-08-04T12:27:23.283906Z","shell.execute_reply":"2024-08-04T12:27:23.29351Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# %%time\n# exec(open('f_read_and_preprocess.py', 'r').read())\n# f_raw_train = f_read_and_preprocess('train', train_adc_info, train_labels.index)\n# with open('f_raw_train.pickle', 'wb') as f:\n#     pickle.dump(f_raw_train, f)\n","metadata":{"papermill":{"duration":1067.258124,"end_time":"2024-08-03T13:00:08.197814","exception":false,"start_time":"2024-08-03T12:42:20.93969","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T10:16:27.461497Z","iopub.execute_input":"2024-08-04T10:16:27.461961Z","iopub.status.idle":"2024-08-04T10:27:09.214198Z","shell.execute_reply.started":"2024-08-04T10:16:27.461928Z","shell.execute_reply":"2024-08-04T10:27:09.212147Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# %%writefile a_read_and_preprocess.py\n# def a_read_and_preprocess(dataset, adc_info, planet_ids):\n#     \"\"\"Read the AIRS-CH0 files for all planet_ids and extract the time series.\n    \n#     Parameters\n#     dataset: 'train' or 'test'\n#     adc_info: metadata dataframe, either train_adc_info or test_adc_info\n#     planet_ids: list of planet ids\n    \n#     Returns\n#     dataframe with one row per planet_id and 5625 values per row\n    \n#     \"\"\"\n#     a_raw_train = np.full((len(planet_ids), 5625), np.nan, dtype=np.float32)\n#     for i, planet_id in tqdm(list(enumerate(planet_ids))):\n#         signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet')\n#         mean_signal = signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / (32*356) # mean over the 32*356 pixels\n#         net_signal = mean_signal[1::2] - mean_signal[0::2]\n#         a_raw_train[i] = net_signal\n#     return a_raw_train\n    ","metadata":{"execution":{"iopub.status.busy":"2024-08-04T12:27:27.791719Z","iopub.execute_input":"2024-08-04T12:27:27.792857Z","iopub.status.idle":"2024-08-04T12:27:27.800511Z","shell.execute_reply.started":"2024-08-04T12:27:27.792811Z","shell.execute_reply":"2024-08-04T12:27:27.799003Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# %%time\n# exec(open('a_read_and_preprocess.py', 'r').read())\n# a_raw_train = a_read_and_preprocess('train', train_adc_info, train_labels.index)\n# with open('a_raw_train.pickle', 'wb') as f:\n#     pickle.dump(a_raw_train, f)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T10:49:43.17432Z","iopub.execute_input":"2024-08-04T10:49:43.1748Z","iopub.status.idle":"2024-08-04T11:10:07.5598Z","shell.execute_reply.started":"2024-08-04T10:49:43.174748Z","shell.execute_reply":"2024-08-04T11:10:07.557469Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"As a plausibility check, we plot the means of all time series:","metadata":{}},{"cell_type":"code","source":"# plt.figure(figsize=(6, 2))\n# plt.plot(f_raw_train.mean(axis=0))\n# for time_step in [20500, 23500, 44000, 47000]:\n#     plt.axvline(time_step, color='gray')\n# plt.xlabel('time step')\n# plt.title('FGS1: Overall mean')\n# plt.show()\n\n# plt.figure(figsize=(6, 2))\n# plt.plot(a_raw_train.mean(axis=0))\n# for time_step in [20500, 23500, 44000, 47000]:\n#     plt.axvline(time_step * 11250 // 135000, color='gray')\n# plt.xlabel('time step')\n# plt.title('AIRS-CH0: Overall mean')\n# plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:21:05.261032Z","iopub.execute_input":"2024-08-04T11:21:05.262209Z","iopub.status.idle":"2024-08-04T11:21:05.753651Z","shell.execute_reply.started":"2024-08-04T11:21:05.262167Z","shell.execute_reply":"2024-08-04T11:21:05.752301Z"},"_kg_hide-input":true,"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature engineering\n\nWe want to know how much darker the images get when the planet obscures the star. The time series diagrams above show that the planets reduce the brightness of the stars (on average) by 0.2 % (from 228.2 to 227.6 or from 1371 to 1368).","metadata":{}},{"cell_type":"code","source":"# %%writefile feature_engineering.py\n\n# def feature_engineering(f_raw, a_raw):\n#     \"\"\"Create a dataframe with two features from the raw data.\n    \n#     Parameters:\n#     f_raw: ndarray of shape (n_planets, 67500)\n#     a_raw: ndarray of shape (n_planets, 5625)\n    \n#     Return value:\n#     df: DataFrame of shape (n_planets, 2)\n#     \"\"\"\n#     obscured = f_raw[:, 23500:44000].mean(axis=1)\n#     unobscured = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\n#     f_relative_reduction = (unobscured - obscured) / unobscured\n#     obscured = a_raw[:, 1958:3666].mean(axis=1)\n#     unobscured = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\n#     a_relative_reduction = (unobscured - obscured) / unobscured\n\n#     df = pd.DataFrame({'a_relative_reduction': a_relative_reduction,\n#                        'f_relative_reduction': f_relative_reduction})\n    \n#     return df\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:49:32.362734Z","iopub.execute_input":"2024-08-04T11:49:32.363268Z","iopub.status.idle":"2024-08-04T11:49:32.371284Z","shell.execute_reply.started":"2024-08-04T11:49:32.363232Z","shell.execute_reply":"2024-08-04T11:49:32.370079Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# exec(open('feature_engineering.py', 'r').read())\n\n# train = feature_engineering(f_raw_train, a_raw_train)","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:33:15.805608Z","iopub.execute_input":"2024-08-04T11:33:15.806053Z","iopub.status.idle":"2024-08-04T11:33:15.833728Z","shell.execute_reply.started":"2024-08-04T11:33:15.806018Z","shell.execute_reply":"2024-08-04T11:33:15.832575Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The following scatterplot shows a strong correlation between the signal reduction when the planet is in front of the star and one of the targets we want to predict. The other targets have a similarly high correlation.","metadata":{"papermill":{"duration":0.073746,"end_time":"2024-08-03T13:00:08.349433","exception":false,"start_time":"2024-08-03T13:00:08.275687","status":"completed"},"tags":[]}},{"cell_type":"code","source":"\n# color_array = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])\n# plt.scatter(train.a_relative_reduction, train_labels.wl_1, s=15, alpha=0.5,\n#             c=color_array[train_adc_info.star])\n# plt.xlabel('relative signal reduction when planet is in front')\n# plt.ylabel('target')\n# plt.title('Correlation between relative signal reduction and target')\n# plt.gca().set_aspect('equal')\n# points = [plt.Line2D([0], [0], label=f'star {i}', marker='o', markersize=3,\n#          markeredgecolor=color_array[i], markerfacecolor=color_array[i], linestyle='') for i in range(2)]\n\n# plt.legend(handles=points)\n# plt.show()","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.478644,"end_time":"2024-08-03T13:00:08.902093","exception":false,"start_time":"2024-08-03T13:00:08.423449","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:18.162145Z","iopub.execute_input":"2024-08-04T11:33:18.162554Z","iopub.status.idle":"2024-08-04T11:33:18.521622Z","shell.execute_reply.started":"2024-08-04T11:33:18.162523Z","shell.execute_reply":"2024-08-04T11:33:18.520319Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# The model and the cross-validation\n\nTo keep things simple, we predict the targets with ridge regression.\n\nWe are interested in three cross-validation metrics:\n1. The R2 score is above 0.9, which confirms the correlation we've seen in the scatterplot.\n2. The root mean squared error will be the predicted uncertainty.\n3. The competition metric gives an indication of the leaderboard score. Unfortunately the competition metric depends on the value of `sigma_true`, which I don't know.","metadata":{"papermill":{"duration":0.076791,"end_time":"2024-08-03T13:00:09.056736","exception":false,"start_time":"2024-08-03T13:00:08.979945","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# model = Ridge(alpha=1e-12)\n\n# oof_pred = cross_val_predict(model, train, train_labels)\n\n# print(f\"# R2 score: {r2_score(train_labels, oof_pred):.3f}\")\n# sigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\n# print(f\"# Root mean squared error: {sigma_pred:.6f}\")\n# # R2 score: 0.971\n# # Root mean squared error: 0.000293","metadata":{"papermill":{"duration":0.700113,"end_time":"2024-08-03T13:00:09.837053","exception":false,"start_time":"2024-08-03T13:00:09.13694","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:20.382678Z","iopub.execute_input":"2024-08-04T11:33:20.383118Z","iopub.status.idle":"2024-08-04T11:33:20.508765Z","shell.execute_reply.started":"2024-08-04T11:33:20.383083Z","shell.execute_reply":"2024-08-04T11:33:20.507443Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# col = 1\n# plt.scatter(oof_pred[:,col], train_labels.iloc[:,col], s=15, c='lightgreen')\n# plt.gca().set_aspect('equal')\n# plt.xlabel('y_pred')\n# plt.ylabel('y_true')\n# plt.title('Comparing y_true and y_pred')\n# plt.show()","metadata":{"papermill":{"duration":0.700113,"end_time":"2024-08-03T13:00:09.837053","exception":false,"start_time":"2024-08-03T13:00:09.13694","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-04T11:33:24.335234Z","iopub.execute_input":"2024-08-04T11:33:24.33565Z","iopub.status.idle":"2024-08-04T11:33:24.628829Z","shell.execute_reply.started":"2024-08-04T11:33:24.335622Z","shell.execute_reply":"2024-08-04T11:33:24.627451Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# %%writefile competition_score.py\n# # Adapted from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\n# class ParticipantVisibleError(Exception):\n#     pass\n\n# def competition_score(\n#         solution: pd.DataFrame,\n#         submission: pd.DataFrame,\n#         naive_mean: float,\n#         naive_sigma: float,\n#         sigma_true: float,\n#         row_id_column_name='planet_id',\n#     ) -> float:\n#     '''\n#     This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n#     we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n#     x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n#     the GLL value for one spectrum is the sum of all of them.\n\n#     Inputs:\n#         - solution: Ground Truth spectra (from test set)\n#             - shape: (nsamples, n_wavelengths)\n#         - submission: Predicted spectra and errors (from participants)\n#             - shape: (nsamples, n_wavelengths*2)\n#         naive_mean: (float) mean from the train set.\n#         naive_sigma: (float) standard deviation from the train set.\n#         sigma_true: (float) essentially sets the scale of the outputs.\n#     '''\n\n#     del solution[row_id_column_name]\n#     del submission[row_id_column_name]\n\n#     if submission.min().min() < 0:\n#         raise ParticipantVisibleError('Negative values in the submission')\n#     for col in submission.columns:\n#         if not pd.api.types.is_numeric_dtype(submission[col]):\n#             raise ParticipantVisibleError(f'Submission column {col} must be a number')\n\n#     n_wavelengths = len(solution.columns)\n#     if len(submission.columns) != n_wavelengths*2:\n#         raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n#     y_pred = submission.iloc[:, :n_wavelengths].values\n#     # Set a non-zero minimum sigma pred to prevent division by zero errors.\n#     sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n#     y_true = solution.values\n\n#     GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n#     GLL_true = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true)))\n#     GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true)))\n\n#     submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n#     return float(np.clip(submit_score, 0.0, 1.0))","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.029243,"end_time":"2024-08-03T12:42:17.162838","exception":false,"start_time":"2024-08-03T12:42:17.133595","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:25.777887Z","iopub.execute_input":"2024-08-04T11:33:25.778327Z","iopub.status.idle":"2024-08-04T11:33:25.787298Z","shell.execute_reply.started":"2024-08-04T11:33:25.778293Z","shell.execute_reply":"2024-08-04T11:33:25.78591Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# %%writefile postprocessing.py\n\n# def postprocessing(pred_array, index, sigma_pred):\n#     \"\"\"Create a submission dataframe from its components\n    \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: float\n    \n#     Return value:\n#     df: DataFrame of shape (n_samples, 566) with planet_id as index\n#     \"\"\"\n#     return pd.concat([pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns),\n#                       pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)])],\n#                      axis=1)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:33:27.453077Z","iopub.execute_input":"2024-08-04T11:33:27.453955Z","iopub.status.idle":"2024-08-04T11:33:27.46076Z","shell.execute_reply.started":"2024-08-04T11:33:27.453914Z","shell.execute_reply":"2024-08-04T11:33:27.459579Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# exec(open('competition_score.py', 'r').read())\n# exec(open('postprocessing.py', 'r').read())\n\n# oof_df = postprocessing(oof_pred, train_adc_info.index, sigma_pred)\n# display(oof_df)\n\n# gll_score = competition_score(train_labels.copy().reset_index(),\n#                               oof_df.copy().reset_index(),\n#                               naive_mean=train_labels.values.mean(),\n#                               naive_sigma=train_labels.values.std(),\n#                               sigma_true=0.000003)\n# print(f\"# Estimated competition score: {gll_score:.3f}\")\n# # Estimated competition score: 0.259","metadata":{"papermill":{"duration":0.231627,"end_time":"2024-08-03T13:00:10.149224","exception":false,"start_time":"2024-08-03T13:00:09.917597","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:28.166304Z","iopub.execute_input":"2024-08-04T11:33:28.166706Z","iopub.status.idle":"2024-08-04T11:33:28.266475Z","shell.execute_reply.started":"2024-08-04T11:33:28.166678Z","shell.execute_reply":"2024-08-04T11:33:28.265186Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Refitting and saving the model","metadata":{}},{"cell_type":"code","source":"# # Refit the model to the full dataset\n# model.fit(train, train_labels)\n# with open('model.pickle', 'wb') as f:\n#     pickle.dump(model, f)\n# with open('sigma_pred.pickle', 'wb') as f:\n#     pickle.dump(sigma_pred, f)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T13:15:17.228909Z","iopub.execute_input":"2024-08-04T13:15:17.229421Z","iopub.status.idle":"2024-08-04T13:15:17.253336Z","shell.execute_reply.started":"2024-08-04T13:15:17.229387Z","shell.execute_reply":"2024-08-04T13:15:17.251897Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submission","metadata":{"papermill":{"duration":0.077479,"end_time":"2024-08-03T13:00:10.305108","exception":false,"start_time":"2024-08-03T13:00:10.227629","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# # Load the data\n# test_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n#                            index_col='planet_id')\n# sample_submission = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv',\n#                                 index_col='planet_id')\n# f_raw_test = f_read_and_preprocess('test', test_adc_info, sample_submission.index)\n# a_raw_test = a_read_and_preprocess('test', test_adc_info, sample_submission.index)\n# test = feature_engineering(f_raw_test, a_raw_test)\n\n# # Load the model\n# with open('model.pickle', 'rb') as f:\n#     model = pickle.load(f)\n# with open('sigma_pred.pickle', 'rb') as f:\n#     sigma_pred = pickle.load(f)\n\n# # Predict\n# test_pred = model.predict(test)\n\n# # Package into submission file\n# sub_df = postprocessing(test_pred, sample_submission.index, sigma_pred)\n# display(sub_df)\n# sub_df.to_csv('submission.csv')\n# #!head submission.csv","metadata":{"papermill":{"duration":1.822768,"end_time":"2024-08-03T13:00:12.204995","exception":false,"start_time":"2024-08-03T13:00:10.382227","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T13:15:48.046438Z","iopub.execute_input":"2024-08-04T13:15:48.046918Z","iopub.status.idle":"2024-08-04T13:15:51.037458Z","shell.execute_reply.started":"2024-08-04T13:15:48.046886Z","shell.execute_reply":"2024-08-04T13:15:51.036278Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{},"outputs":[],"execution_count":null}]}