{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"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"},"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":"\n## Introduction to Fast Data Cleaning Approach\nIn data science, efficient data cleaning is crucial for handling large datasets, especially in fields like astronomy where data volumes can be immense. The traditional approach using masked arrays allows for flexibility in masking invalid or missing data, but it can introduce performance bottlenecks, particularly with large-scale datasets.\n\nIn the this notebook, I adopt a more performant approach by utilizing NaN (Not a Number) values for masking. This method leverages the optimized operations available in NumPy and Pandas, which handle NaN values efficiently and allow for faster computations compared to traditional masked arrays. By replacing masked elements with NaNs, I streamline data processing and improve the speed of operations like filtering, aggregating, and statistical analysis.\n\nAdditionally, I keep the data dimensions two dimensional while processing. This makes it slightly harder to reason about dimensions during data processing, but ensures better data-locality while processing the large files of this challenge.\n\nThis code is inspired by others in this competition. Enormous thanks goes to: @AmbrosM and @GordonYip.\n\nCode licensed under [\"Rejoice of Fun License (ROFL)\"](https://lmy.medium.com/rofl-an-open-source-license-that-promotes-fun-in-coding-620388502891). Attribution appreciated :-)","metadata":{}},{"cell_type":"code","source":"import pickle\nfrom pathlib import Path\nfrom tqdm import tqdm\n\nimport numpy as np\nimport pandas as pd\nimport scipy.stats\nfrom astropy.stats import sigma_clip\n\nimport matplotlib.pyplot as plt","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":[],"execution":{"iopub.status.busy":"2024-08-17T13:57:57.249194Z","iopub.execute_input":"2024-08-17T13:57:57.249664Z","iopub.status.idle":"2024-08-17T13:57:57.256576Z","shell.execute_reply.started":"2024-08-17T13:57:57.24963Z","shell.execute_reply":"2024-08-17T13:57:57.255368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read metadata\n\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv', index_col='planet_id')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv', index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')","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":"2024-08-17T13:47:48.470611Z","iopub.execute_input":"2024-08-17T13:47:48.471197Z","iopub.status.idle":"2024-08-17T13:47:48.822537Z","shell.execute_reply.started":"2024-08-17T13:47:48.471164Z","shell.execute_reply":"2024-08-17T13:47:48.821296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile preprocess_fast.py\n\n\ndef ADC_convert(signal, gain, offset):\n    signal /= gain\n    signal += offset\n    return signal\n\ndef mask_hot_dead(signal, dead, dark):\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    hot_mask = hot.reshape((-1,))\n    dead_mask = (dead == 1.0).reshape((-1,))\n\n    signal[:, hot_mask] = np.NaN\n    signal[:, dead_mask] = np.NaN\n    return signal\n\ndef clean_dark(signal, dark, dt):\n    dark_current = (dt[:, np.newaxis] * dark)\n    signal -= dark_current\n    return signal\n\ndef clean_flat(signal, flat):\n    signal = (signal) / (flat)\n    return signal\n\ndef apply_linear_corr(linear_corr,signal):\n    for i in range(signal.shape[1]):\n        poli = np.poly1d(np.flip(linear_corr[:, i]))\n        signal[:, i] = poli(signal[:, i])\n    return signal\n\ndef bin_obs(signal ,binning):\n    signal_binned = np.zeros((signal.shape[0]//binning, signal.shape[1]))\n    for i in range(signal.shape[0]//binning):\n        signal_binned[i, :] = np.mean(signal[i*binning:(i+1)*binning, :], axis=0)\n    return signal_binned\n\ndef airs_preprocess(dataset, adc_info, axis_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    axis_info: axis info, includes gain/offset/integration\n    planet_ids: list of planet ids\n    \n    Returns\n    dataframe with one row per planet_id and 5625//binning values per row\n    \n    \"\"\"\n    binning = 60\n    \n    AIRS_CH0_gain, AIRS_CH0_offset = adc_info['AIRS-CH0_adc_gain'].values, adc_info['AIRS-CH0_adc_offset'].values\n    dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n    \n    # planet x time-binned data x frequency\n    a_raw = np.full((len(planet_ids), 5625//binning, 356), np.nan, dtype=np.float32)\n    \n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet').values.astype(np.float64)\n        flat = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_calibration/flat.parquet').values.astype(np.float64).reshape((1, -1))\n        dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_calibration/dark.parquet').values.astype(np.float64).reshape((1, -1))\n        dead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_calibration/dead.parquet').values.astype(np.float64).reshape((1, -1))\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_calibration/linear_corr.parquet').values.astype(np.float64).reshape(6,-1)\n        \n        signal = ADC_convert(signal, AIRS_CH0_gain[i], AIRS_CH0_offset[i])\n        \n        signal = mask_hot_dead(signal, dead, dark) # with NaN, avoiding use of masked arrays\n        \n        signal = apply_linear_corr(linear_corr,signal) # operating on non-masked *much* faster\n        \n        signal = clean_dark(signal, dark, dt_airs)\n        \n        signal = signal[1::2, :] - signal[0::2, :] # Correlated double sampling\n        signal = bin_obs(signal, binning)\n        \n        signal = clean_flat(signal, flat)\n        \n        # mean over pixels, preserve wavelength dimension\n        signal = signal.reshape((signal.shape[0], 32, 356))\n        signal = np.nanmean(signal, axis=1) # nanmean because we're masking with nan\n\n        a_raw[i] = signal\n    return a_raw\n\ndef fgs_preprocess(dataset, adc_info, axis_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    axis_info: axis info, includes gain/offset/integration\n    planet_ids: list of planet ids\n    \n    Returns\n    dataframe with one row per planet_id and 67500//binning values per row\n    \n    \"\"\"\n    binning = 720\n        \n    # planet x time-binned data\n    fgs_raw = np.full((len(planet_ids), 67500//binning), np.nan, dtype=np.float32)\n    \n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_signal.parquet').values.astype(np.float64)\n        flat = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_calibration/flat.parquet').values.astype(np.float64).reshape((1, -1))\n        dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_calibration/dark.parquet').values.astype(np.float64).reshape((1, -1))\n        dead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_calibration/dead.parquet').values.astype(np.float64).reshape((1, -1))\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_calibration/linear_corr.parquet').values.astype(np.float64).reshape(6,-1)\n                \n        FGS1_gain = train_adc_info['FGS1_adc_gain'].values[i]\n        FGS1_offset = train_adc_info['FGS1_adc_offset'].values[i]\n        \n        signal = ADC_convert(signal, FGS1_gain, FGS1_offset)\n            \n        signal = mask_hot_dead(signal, dead, dark) # with NaN, avoiding use of masked arrays\n        \n        signal = apply_linear_corr(linear_corr,signal) # operating on non-masked *much* faster\n        \n        dt_fgs1 = np.ones(len(signal))*0.1\n        signal = clean_dark(signal, dark, dt_fgs1)\n        \n        signal = signal[1::2, :] - signal[0::2, :] # Correlated double sampling\n        signal = bin_obs(signal, binning)\n        \n        signal = clean_flat(signal, flat)\n        \n        # mean over pixels\n        signal = np.nanmean(signal, axis=1) # nanmean because we're masking with nan\n\n        fgs_raw[i] = signal\n    return fgs_raw\n    ","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:09:02.219544Z","iopub.execute_input":"2024-08-17T14:09:02.219961Z","iopub.status.idle":"2024-08-17T14:09:02.231476Z","shell.execute_reply.started":"2024-08-17T14:09:02.219931Z","shell.execute_reply":"2024-08-17T14:09:02.230219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nexec(open('preprocess_fast.py', 'r').read())\nprint('AIRS Data:')\nairs_raw_train = airs_preprocess('train', train_adc_info, axis_info, train_labels.index[0:5])\nprint('FGS Data:')\nfgs_raw_train = fgs_preprocess('train', train_adc_info, axis_info, train_labels.index[0:5])","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:09:08.30934Z","iopub.execute_input":"2024-08-17T14:09:08.309888Z","iopub.status.idle":"2024-08-17T14:09:50.517889Z","shell.execute_reply.started":"2024-08-17T14:09:08.309853Z","shell.execute_reply":"2024-08-17T14:09:50.51657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Obligatory Light Curve Plots\n\nJust to make sure we didn't break something too badly... 🙃","metadata":{}},{"cell_type":"code","source":"for i in range(5):\n    plt.plot((fgs_raw_train[i,:])/np.mean(fgs_raw_train[i,:]))","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:05:39.210089Z","iopub.execute_input":"2024-08-17T14:05:39.210944Z","iopub.status.idle":"2024-08-17T14:05:39.498297Z","shell.execute_reply.started":"2024-08-17T14:05:39.21091Z","shell.execute_reply":"2024-08-17T14:05:39.496908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(5):\n    plt.plot((np.mean(airs_raw_train[i,:, :], axis=1)/np.mean(airs_raw_train[i,:, :])))","metadata":{"execution":{"iopub.status.busy":"2024-08-17T14:08:52.318014Z","iopub.execute_input":"2024-08-17T14:08:52.318405Z","iopub.status.idle":"2024-08-17T14:08:52.5422Z","shell.execute_reply.started":"2024-08-17T14:08:52.318378Z","shell.execute_reply":"2024-08-17T14:08:52.54103Z"},"trusted":true},"execution_count":null,"outputs":[]}]}