{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","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"}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd \nimport numpy as np \nimport matplotlib.pyplot as plt \nimport torch \nimport itertools \nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-10T17:08:36.815974Z","iopub.execute_input":"2024-10-10T17:08:36.816415Z","iopub.status.idle":"2024-10-10T17:08:36.821635Z","shell.execute_reply.started":"2024-10-10T17:08:36.816372Z","shell.execute_reply":"2024-10-10T17:08:36.820502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Testing Raw Image ","metadata":{}},{"cell_type":"markdown","source":"### FGS1 Image","metadata":{}},{"cell_type":"code","source":"planet_id = 499191466\npath_to_dataset = '/kaggle/input/ariel-data-challenge-2024'\nfgs1_signal_path = os.path.join(path_to_dataset, 'test', str(planet_id), 'FGS1_signal.parquet')\nfgs1_signal = pd.read_parquet(fgs1_signal_path)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:36.823929Z","iopub.execute_input":"2024-10-10T17:08:36.824997Z","iopub.status.idle":"2024-10-10T17:08:37.345577Z","shell.execute_reply.started":"2024-10-10T17:08:36.824939Z","shell.execute_reply":"2024-10-10T17:08:37.344135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualizing the first 32x32 image (first row reshaped)\nimage_data = fgs1_signal.iloc[0].values.reshape(32, 32)\n\n# Ploting the image using matplotlib\nplt.imshow(image_data, cmap='hot')\nplt.colorbar()\nplt.title('FGS1 Signal Image (First Row)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:37.34827Z","iopub.execute_input":"2024-10-10T17:08:37.348782Z","iopub.status.idle":"2024-10-10T17:08:37.713916Z","shell.execute_reply.started":"2024-10-10T17:08:37.348726Z","shell.execute_reply":"2024-10-10T17:08:37.712348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### AIRS-CHO Image","metadata":{}},{"cell_type":"code","source":"AIRS_CH0_signal_path = os.path.join(path_to_dataset, 'train', str(100468857), 'AIRS-CH0_signal.parquet')\nAIRS_CH0_signal = pd.read_parquet(AIRS_CH0_signal_path)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:37.715529Z","iopub.execute_input":"2024-10-10T17:08:37.715901Z","iopub.status.idle":"2024-10-10T17:08:38.971416Z","shell.execute_reply.started":"2024-10-10T17:08:37.715865Z","shell.execute_reply":"2024-10-10T17:08:38.970181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_data = AIRS_CH0_signal.iloc[0].values.reshape(32, 356)\n\n# Using MatPlotlib \nplt.imshow(image_data, cmap='viridis', aspect='auto')\nplt.colorbar(label='Signal Intensity')\nplt.title(f'Combine AIRS-CH0 Signal Data for Planet {planet_id}')\nplt.xlabel('Time Step * Pixel Y')\nplt.ylabel('Pixel X')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:38.974682Z","iopub.execute_input":"2024-10-10T17:08:38.975182Z","iopub.status.idle":"2024-10-10T17:08:39.480778Z","shell.execute_reply.started":"2024-10-10T17:08:38.975126Z","shell.execute_reply":"2024-10-10T17:08:39.479405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Using MatPlotlib \nplt.imshow(image_data, cmap='viridis')\nplt.colorbar(label='Signal Intensity')\nplt.title(f'Combine AIRS-CH0 Signal Data for Planet {planet_id}')\nplt.xlabel('Time Step * Pixel Y')\nplt.ylabel('Pixel X')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.482339Z","iopub.execute_input":"2024-10-10T17:08:39.482834Z","iopub.status.idle":"2024-10-10T17:08:39.883849Z","shell.execute_reply.started":"2024-10-10T17:08:39.482791Z","shell.execute_reply":"2024-10-10T17:08:39.882592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calibration \n\n1. Analog-to-Digital Conversion (ADC) Correction: Gain and Offset \n2. Masking Hot/Dead Pixels \n3. Dark Current Subtraction \n4. Linearity Correction \n5. Correlated Double Sampling (CDS) \n6. Flat Field Correction \n7. Time Binning (Optional)","metadata":{}},{"cell_type":"markdown","source":"# Aanalog-to-Digital Conversion (ADC): Applying Gain and Offset\n\n- **Purpose**: ADC conversion transforms raw pixel voltage readings into meaningful physical units (e.g., electrons or photons). \n\n- **Reasoning**: Applying gain and offset first ensures that all subsequent calibration steps operate on data in the correct physical scale. This is crucial because calibration parameters(like dark current or flat fields) are typically defined in these physical units. ","metadata":{}},{"cell_type":"code","source":"def ADC_convert(signal, gain, offset): \n    ''' \n    Parameters: \n        signal (np.ndarrray): Raw signal data. \n        gain (Float): Gain values for ADC conversion \n        Offset (Float): Offset value for ADC conversion.\n        \n    Returns: \n        np.ndarray: Calibrated signal in physical units.\n        \n    '''\n    \n    signal = signal.astype(np.float64)\n    signal = signal * gain + offset \n    return signal ","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.885523Z","iopub.execute_input":"2024-10-10T17:08:39.886007Z","iopub.status.idle":"2024-10-10T17:08:39.893682Z","shell.execute_reply.started":"2024-10-10T17:08:39.885948Z","shell.execute_reply":"2024-10-10T17:08:39.892003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Masking Hot/Dead Pixels ","metadata":{}},{"cell_type":"code","source":"from astropy.stats import sigma_clip \n\ndef mask_hot_dead(signal, dead, dark):\n    \"\"\"\n    Masks dead and hot pixels in the signal data.\n    \n    Parameters:\n        signal (np.ndarray): Signal data.\n        dead (np.ndarray): Dead pixel map.\n        dark (np.ndarray): Dark frame for identifying hot pixels.\n        \n    Returns:\n        np.ma.MaskedArray: Masked signal data.\n    \"\"\"\n    \n    # Identifying hot pixels using sigma clipping on the dark frame \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    \n    # Mask dead and hot pixels \n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal) \n    return signal \n    ","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.895603Z","iopub.execute_input":"2024-10-10T17:08:39.896265Z","iopub.status.idle":"2024-10-10T17:08:39.909126Z","shell.execute_reply.started":"2024-10-10T17:08:39.896111Z","shell.execute_reply":"2024-10-10T17:08:39.907638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dark Current Subtraction","metadata":{}},{"cell_type":"code","source":"def clean_dark(signal, dead, dark, dt): \n    \"\"\"\n    Subtracts dark current from the signal.\n    \n    Parameters:\n        signal (np.ma.MaskedArray): Signal data with masked bad pixels.\n        dead (np.ndarray): Dead pixel map.\n        dark (np.ndarray): Dark current map.\n        dt (np.ndarray): Integration time for each frame.\n        \n    Returns:\n        np.ma.MaskedArray: Dark-subtracted signal.\n    \"\"\"\n    dark = np.ma.masked_where(dead, dark)\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n    signal -= dark * dt[:, np.newaxis, np.newaxis]\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.910989Z","iopub.execute_input":"2024-10-10T17:08:39.911627Z","iopub.status.idle":"2024-10-10T17:08:39.942487Z","shell.execute_reply.started":"2024-10-10T17:08:39.911566Z","shell.execute_reply":"2024-10-10T17:08:39.940902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Linearity Correction","metadata":{}},{"cell_type":"code","source":"def apply_linear_corr(linear_corr, clean_signal): \n    \"\"\"\n    Applies linearity correction to the signal data.\n    \n    Parameters:\n        linear_corr (np.ndarray): Coefficients for the inverse polynomial.\n        clean_signal (np.ma.MaskedArray): Dark-subtracted signal.\n        \n    Returns:\n        np.ma.MaskedArray: Linearity-corrected signal.\n    \"\"\"\n    \"\"\"\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(range(clean_signal.shape[1]), range(clean_signal.shape[2])): \n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal\n    \"\"\"\n    \n    # Ensure linear_corr is flipped along axis 0 to get proper polynomial coefficients\n    linear_corr = np.flip(linear_corr, axis=0)  # Shape: (degree, x, y)\n\n    # Get the shape of the signal\n    time, x, y = clean_signal.shape\n\n    # Initialize the corrected signal array\n    corrected_signal = np.zeros_like(clean_signal, dtype=np.float64)\n\n    # Apply polynomial correction for each pixel (x, y) across all time steps\n    for i in range(x):\n        for j in range(y):\n            # Extract the polynomial coefficients for pixel (i, j)\n            coeffs = linear_corr[:, i, j]\n            # Apply polynomial correction for all time steps for pixel (i, j)\n            corrected_signal[:, i, j] = np.polyval(coeffs, clean_signal[:, i, j])\n\n    # Maintain the mask from the input signal\n    corrected_signal = np.ma.MaskedArray(corrected_signal, mask=clean_signal.mask)\n\n    return corrected_signal","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.944141Z","iopub.execute_input":"2024-10-10T17:08:39.94474Z","iopub.status.idle":"2024-10-10T17:08:39.956946Z","shell.execute_reply.started":"2024-10-10T17:08:39.944683Z","shell.execute_reply":"2024-10-10T17:08:39.955598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlated Double Sampling (CDS)","metadata":{}},{"cell_type":"code","source":"def get_cds(signal):\n    \"\"\"\n    Applies Correlated Double Sampling (CDS) to the signal data.\n    \n    Parameters:\n        signal (np.ndarray): Linearity-corrected signal (3D: time, x, y).\n        \n    Returns:\n        np.ndarray: CDS-processed signal.\n    \"\"\"\n    # Perform CDS by subtracting alternating frames\n    cds = signal[1::2, :, :] - signal[::2, :, :]\n    return cds","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.96182Z","iopub.execute_input":"2024-10-10T17:08:39.962735Z","iopub.status.idle":"2024-10-10T17:08:39.976926Z","shell.execute_reply.started":"2024-10-10T17:08:39.962681Z","shell.execute_reply":"2024-10-10T17:08:39.975538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Flat Field Correction","metadata":{}},{"cell_type":"code","source":"def correct_flat_field(flat, dead, signal): \n    \"\"\"\n    Applies flat field correction to the signal data.\n    \n    Parameters:\n        flat (np.ndarray): Flat field map.\n        dead (np.ndarray): Dead pixel map.\n        signal (np.ndarray): CDS-processed signal.\n        \n    Returns:\n        np.ndarray: Flat field-corrected signal.\n    \"\"\"\n    # Mask out dead pixels in the flat field\n    flat = np.ma.masked_where(dead, flat)\n    \n    # Tile the flat field to match the number of time steps in the signal\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    \n    # Apply flat field correction by dividing the signal by the flat field\n    signal = signal / flat \n    \n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.978908Z","iopub.execute_input":"2024-10-10T17:08:39.979835Z","iopub.status.idle":"2024-10-10T17:08:39.991965Z","shell.execute_reply.started":"2024-10-10T17:08:39.979774Z","shell.execute_reply":"2024-10-10T17:08:39.990501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Testing Calibration Workflow","metadata":{}},{"cell_type":"code","source":"# File paths\npath_to_dataset = '/kaggle/input/ariel-data-challenge-2024'\nplanet_id = 100468857\n\nfgs1_signal_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'FGS1_signal.parquet')\nairs_signal_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_signal.parquet')","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:39.993543Z","iopub.execute_input":"2024-10-10T17:08:39.995123Z","iopub.status.idle":"2024-10-10T17:08:40.015883Z","shell.execute_reply.started":"2024-10-10T17:08:39.995035Z","shell.execute_reply":"2024-10-10T17:08:40.013689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the raw signals\nfgs1_signal = pd.read_parquet(fgs1_signal_path)\nairs_signal = pd.read_parquet(airs_signal_path)\n\n# Load the train_adc_info.csv file\ntrain_adc_info = pd.read_csv(os.path.join(path_to_dataset, 'train_adc_info.csv'))\n\n# Set 'planet_id' as the index of the DataFrame\ntrain_adc_info = train_adc_info.set_index('planet_id')\n\n# Access the gain and offset values for AIRS-CH0 using the planet ID\ngain_airs = train_adc_info.loc[planet_id, 'AIRS-CH0_adc_gain']\noffset_airs = train_adc_info.loc[planet_id, 'AIRS-CH0_adc_offset']","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:40.017579Z","iopub.execute_input":"2024-10-10T17:08:40.018034Z","iopub.status.idle":"2024-10-10T17:08:41.556558Z","shell.execute_reply.started":"2024-10-10T17:08:40.017975Z","shell.execute_reply":"2024-10-10T17:08:41.555105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply ADC conversion\nairs_signal_calibrated = ADC_convert(airs_signal.values, gain_airs, offset_airs)\nprint(\"ADC conversion done. Shape of AIRS signal:\", airs_signal_calibrated.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:41.557921Z","iopub.execute_input":"2024-10-10T17:08:41.558325Z","iopub.status.idle":"2024-10-10T17:08:42.414934Z","shell.execute_reply.started":"2024-10-10T17:08:41.558284Z","shell.execute_reply":"2024-10-10T17:08:42.413533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Masking Hot/Dead Pixels for AIRS-CH0 \n# Load dead and dark frames\ndead_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dead.parquet')).values.reshape((32, 356))[:, 39:322]\ndark_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dark.parquet')).values.reshape((32, 356))[:, 39:322]","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:42.41632Z","iopub.execute_input":"2024-10-10T17:08:42.4168Z","iopub.status.idle":"2024-10-10T17:08:42.469576Z","shell.execute_reply.started":"2024-10-10T17:08:42.416748Z","shell.execute_reply":"2024-10-10T17:08:42.468057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Reshape the raw signal to its original shape (11250, 32, 356) before wavelength cut\nairs_signal_reshaped_full = airs_signal_calibrated.reshape(11250, 32, 356)\n\n# Modify the wavelength cut to select 283 pixels (keeping one more pixel)\nairs_signal_reshaped = airs_signal_reshaped_full[:, :, 39:322]  # From 39 to 322 to get 283 pixels\n\n# Now apply the masking to the reshaped signal\nairs_signal_masked = mask_hot_dead(airs_signal_reshaped, dead_airs, dark_airs)\n\n# Check the result\nprint(\"Masking applied. Shape of AIRS signal after masking:\", airs_signal_masked.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:42.471054Z","iopub.execute_input":"2024-10-10T17:08:42.471574Z","iopub.status.idle":"2024-10-10T17:08:45.31799Z","shell.execute_reply.started":"2024-10-10T17:08:42.471516Z","shell.execute_reply":"2024-10-10T17:08:45.316749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load integration time from axis_info.parquet\naxis_info = pd.read_parquet(os.path.join(path_to_dataset, 'axis_info.parquet'))\n\n# Extract integration time for AIRS-CH0 (assuming it's available under the column 'AIRS-CH0-integration_time')\ndt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n\n# Adjust integration time according to the 29.08 update, adding 0.1s to every second value\ndt_airs[1::2] += 0.1\n\n# Apply dark current subtraction\nairs_signal_dark_subtracted = clean_dark(airs_signal_masked, dead_airs, dark_airs, dt_airs)\n\n# Check the result\nprint(\"Dark current subtraction applied. Shape of AIRS signal after subtraction:\", airs_signal_dark_subtracted.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:45.319441Z","iopub.execute_input":"2024-10-10T17:08:45.319808Z","iopub.status.idle":"2024-10-10T17:08:46.874691Z","shell.execute_reply.started":"2024-10-10T17:08:45.31977Z","shell.execute_reply":"2024-10-10T17:08:46.87352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the linearity correction coefficients from the calibration file\nlinear_corr_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/linear_corr.parquet')).values.reshape((6, 32, 356))[:, :, 39:322]\n\n# Apply linearity correction to the dark-subtracted signal\nairs_signal_linear_corrected = apply_linear_corr(linear_corr_airs, airs_signal_dark_subtracted)\n\n# Check the result\nprint(\"Linearity correction applied. Shape of AIRS signal after correction:\", airs_signal_linear_corrected.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:46.876003Z","iopub.execute_input":"2024-10-10T17:08:46.876423Z","iopub.status.idle":"2024-10-10T17:08:58.169533Z","shell.execute_reply.started":"2024-10-10T17:08:46.876339Z","shell.execute_reply":"2024-10-10T17:08:58.16795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply CDS to the linearity-corrected signal\nairs_signal_cds = get_cds(airs_signal_linear_corrected)\n\n# Check the result\nprint(\"Correlated double sampling applied. Shape of AIRS signal after CDS:\", airs_signal_cds.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:58.171635Z","iopub.execute_input":"2024-10-10T17:08:58.172199Z","iopub.status.idle":"2024-10-10T17:08:58.443881Z","shell.execute_reply.started":"2024-10-10T17:08:58.172127Z","shell.execute_reply":"2024-10-10T17:08:58.442498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flat_field_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/flat.parquet')\nflat = pd.read_parquet(flat_field_path).values.reshape((32, 356))[:, 39:322]\ndead_pixel_map_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dead.parquet')\ndead = pd.read_parquet(dead_pixel_map_path).values.reshape((32, 356))[:, 39:322]\nairs_signal_flat_corrected = correct_flat_field(flat, dead, airs_signal_cds)\n\n# Check the result\nprint(\"Flat field correction applied. Shape of AIRS signal after flat field correction:\", airs_signal_flat_corrected.shape)","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:08:58.445424Z","iopub.execute_input":"2024-10-10T17:08:58.445818Z","iopub.status.idle":"2024-10-10T17:09:00.285007Z","shell.execute_reply.started":"2024-10-10T17:08:58.445777Z","shell.execute_reply":"2024-10-10T17:09:00.281206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nplt.imshow(airs_signal_flat_corrected[0, :, :], cmap='viridis', aspect='auto')\nplt.colorbar()\nplt.title('Calibrated Signal (After All Corrections)')\nplt.xlabel('Wavelength Pixel')\nplt.ylabel('Spatial Pixel')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:09:00.286539Z","iopub.execute_input":"2024-10-10T17:09:00.28694Z","iopub.status.idle":"2024-10-10T17:09:00.697956Z","shell.execute_reply.started":"2024-10-10T17:09:00.286898Z","shell.execute_reply":"2024-10-10T17:09:00.696779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check the mask in the calibrated signal\nmask = np.ma.getmask(airs_signal_flat_corrected)\n\n# Print the number of masked (hot/dead) pixels\nprint(f\"Number of masked pixels: {mask.sum()}\")\n\n# Visualize the mask (where 1 indicates a masked pixel)\nplt.imshow(mask[0, :, :], cmap='gray', aspect='auto')\nplt.colorbar()\nplt.title('Masked Pixels in the Calibrated Signal')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:09:00.699431Z","iopub.execute_input":"2024-10-10T17:09:00.69979Z","iopub.status.idle":"2024-10-10T17:09:01.079378Z","shell.execute_reply.started":"2024-10-10T17:09:00.699754Z","shell.execute_reply":"2024-10-10T17:09:01.077998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.ndimage import gaussian_filter\n\n# Fill masked values using Gaussian filtering for interpolation\nsignal_filled = airs_signal_flat_corrected.filled(gaussian_filter(airs_signal_flat_corrected, sigma=1))\n\n# Plot the filled signal\nplt.imshow(signal_filled[0, :, :], cmap='viridis', aspect='auto')\nplt.colorbar()\nplt.title('Calibrated Signal with Interpolated Masked Pixels')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:09:01.080822Z","iopub.execute_input":"2024-10-10T17:09:01.081215Z","iopub.status.idle":"2024-10-10T17:09:06.663678Z","shell.execute_reply.started":"2024-10-10T17:09:01.081172Z","shell.execute_reply":"2024-10-10T17:09:06.662097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Workflow","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\nfrom scipy.fft import fft\n\npath_to_dataset = '/kaggle/input/ariel-data-challenge-2024'\ntrain_dir = os.path.join(path_to_dataset, 'train')\n\nmax_planets_to_process = 5 \n\n# Placeholder for the complete feature set \nall_planet_features = []\n\n# Iterate through all planet IDs in the dataset \nplanet_ids = [planet for planet in os.listdir(train_dir) if os.path.isdir(os.path.join(train_dir, planet))]\nplanet_ids = planet_ids[:max_planets_to_process]\n\n# Iterate over each planet ID to perform calibration and feature extraction \nfor planet_id in planet_ids: \n    try: \n        # File paths for the current planet \n        fgs1_signal_path = os.path.join(train_dir, planet_id, 'FGS1_signal.parquet')\n        airs_signal_path = os.path.join(train_dir, planet_id, 'AIRS-CH0_signal.parquet')\n        \n        # Load the signal data for the current planet\n        fgs1_signal = pd.read_parquet(fgs1_signal_path)\n        airs_signal = pd.read_parquet(airs_signal_path)\n        \n        # Apply ADC conversion\n        airs_signal_calibrated = ADC_convert(airs_signal.values, gain_airs, offset_airs)\n        print(\"ADC conversion done. Shape of AIRS signal:\", airs_signal_calibrated.shape)\n        \n        # Masking Hot/Dead Pixels for AIRS-CH0 \n        # Load dead and dark frames\n        dead_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dead.parquet')).values.reshape((32, 356))[:, 39:322]\n        dark_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dark.parquet')).values.reshape((32, 356))[:, 39:322]\n        \n        # Reshape the raw signal to its original shape (11250, 32, 356) before wavelength cut\n        airs_signal_reshaped_full = airs_signal_calibrated.reshape(11250, 32, 356)\n\n        # Modify the wavelength cut to select 283 pixels (keeping one more pixel)\n        airs_signal_reshaped = airs_signal_reshaped_full[:, :, 39:322]  # From 39 to 322 to get 283 pixels\n\n        # Now apply the masking to the reshaped signal\n        airs_signal_masked = mask_hot_dead(airs_signal_reshaped, dead_airs, dark_airs)\n\n        # Check the result\n        print(\"Masking applied. Shape of AIRS signal after masking:\", airs_signal_masked.shape)\n        \n        # Load integration time from axis_info.parquet\n        axis_info = pd.read_parquet(os.path.join(path_to_dataset, 'axis_info.parquet'))\n\n        # Extract integration time for AIRS-CH0 (assuming it's available under the column 'AIRS-CH0-integration_time')\n        dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n\n        # Adjust integration time according to the 29.08 update, adding 0.1s to every second value\n        dt_airs[1::2] += 0.1\n\n        # Apply dark current subtraction\n        airs_signal_dark_subtracted = clean_dark(airs_signal_masked, dead_airs, dark_airs, dt_airs)\n\n        # Check the result\n        print(\"Dark current subtraction applied. Shape of AIRS signal after subtraction:\", airs_signal_dark_subtracted.shape)\n        \n        # Load the linearity correction coefficients from the calibration file\n        linear_corr_airs = pd.read_parquet(os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/linear_corr.parquet')).values.reshape((6, 32, 356))[:, :, 39:322]\n\n        # Apply linearity correction to the dark-subtracted signal\n        airs_signal_linear_corrected = apply_linear_corr(linear_corr_airs, airs_signal_dark_subtracted)\n\n        # Check the result\n        print(\"Linearity correction applied. Shape of AIRS signal after correction:\", airs_signal_linear_corrected.shape)\n        \n        # Apply CDS to the linearity-corrected signal\n        airs_signal_cds = get_cds(airs_signal_linear_corrected)\n\n        # Check the result\n        print(\"Correlated double sampling applied. Shape of AIRS signal after CDS:\", airs_signal_cds.shape)\n        \n        flat_field_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/flat.parquet')\n        flat = pd.read_parquet(flat_field_path).values.reshape((32, 356))[:, 39:322]\n        dead_pixel_map_path = os.path.join(path_to_dataset, 'train', str(planet_id), 'AIRS-CH0_calibration/dead.parquet')\n        dead = pd.read_parquet(dead_pixel_map_path).values.reshape((32, 356))[:, 39:322]\n        airs_signal_flat_corrected = correct_flat_field(flat, dead, airs_signal_cds)\n\n        # Check the result\n        print(\"Flat field correction applied. Shape of AIRS signal after flat field correction:\", airs_signal_flat_corrected.shape)\n        \n        # Flatten the 3D array into 2D: (n_samples, n_features)\n        n_samples, n_x, n_y = airs_signal_flat_corrected.shape\n        flattened_data = airs_signal_flat_corrected.reshape(n_samples, -1)\n        \n        # Step 1: Normalization and Standardization \n        scaler = StandardScaler()\n        normalized_data = scaler.fit_transform(flattened_data)\n        \n        # Step 2: Dimensionality Reduction using PCA \n        # Apply PCA to reduce dimensionality while retaining important vairance \n        pca = PCA(n_components=50) # retain 50 principal components \n        principal_components = pca.fit_transform(normalized_data)\n        \n        # Step 3: Generate Time-Series Features \n        # Calculate mean, variance, and Fourier Transform components \n        time_series_features = []\n        \n        for i in range(n_samples): \n            signal_sample = airs_signal_flat_corrected[i]\n            \n            # Mean and Variance \n            mean_value = np.mean(signal_sample)\n            variance_value = np.var(signal_sample)\n            \n            # Apply FFT to detect periodic components (E.g., Jitter)\n            fft_values = np.abs(fft(signal_sample.flatten()))[:100]\n            \n            # Combine all extracted features \n            features = [mean_value, variance_value] + list(fft_values)\n            time_series_features.append(features)\n            \n        # Append features for current planet to the master list \n        all_planet_features.append({\n            'planet_id': planet_id, \n            'features': time_series_features, \n            'principal_components': principal_components\n        })\n        \n    except Exception as e: \n        print(f\"Error processing planet {planet_id}: {e}\")\n              \n# Convert the collected features to a dataframe for easier manipulation \nfinal_features = []\n\nfor planet in all_planet_features:\n              planet_id = planet['planet_id']\n              features_df = pd.DataFrame(planet['features'])\n              features_df['planet_id'] = planet_id \n              final_features.append(features_df)\n# Concatenate all features Dataframes \nall_features_df = pd.concat(final_features, ignore_index=True)\n\n# Save the extracted features to a CSV file for further analysis\nall_features_df.to_csv('extracted_features_for_all_planets.csv', index=False)\n\nprint(\"Features have been saved to 'extracted_features_for_all_planets.csv'.\")\n# Display the first few rows of the extracted features DataFrame for a quick check\nprint(all_features_df.head())","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:50:07.041154Z","iopub.execute_input":"2024-10-10T17:50:07.041722Z","iopub.status.idle":"2024-10-10T17:52:51.40949Z","shell.execute_reply.started":"2024-10-10T17:50:07.041677Z","shell.execute_reply":"2024-10-10T17:52:51.407996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_features_df.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-10T17:53:14.230938Z","iopub.execute_input":"2024-10-10T17:53:14.231409Z","iopub.status.idle":"2024-10-10T17:53:14.264527Z","shell.execute_reply.started":"2024-10-10T17:53:14.231362Z","shell.execute_reply":"2024-10-10T17:53:14.263383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}