{"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.14"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":205005104,"sourceType":"kernelVersion"}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":22.138909,"end_time":"2024-10-31T14:47:46.242774","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-10-31T14:47:24.103865","version":"2.6.0"},"widgets":{"application/vnd.jupyter.widget-state+json":{"state":{"00b8162669e7414e90a014cfca8f0112":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"ProgressStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"ProgressStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","bar_color":null,"description_width":""}},"0238f920abea48309804c8de0b48eb1a":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"041a24d2e34241f68c5c8ef8a32b0243":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HTMLModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HTMLModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HTMLView","description":"","description_tooltip":null,"layout":"IPY_MODEL_b9340b8f100d4aba9591a8c2b989404e","placeholder":"​","style":"IPY_MODEL_7356325c4b1345998b8624f16dd4076a","value":"100%"}},"053831b861454c6899ad5e69465a6e1f":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"FloatProgressModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"FloatProgressModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"ProgressView","bar_style":"success","description":"","description_tooltip":null,"layout":"IPY_MODEL_787becb0bd6348d8a04e2ef198c4e410","max":1,"min":0,"orientation":"horizontal","style":"IPY_MODEL_9bd136706b75460ea87a95dd4f1e3e5b","value":1}},"16b39aad2d854563988feec0dd4f1769":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"DescriptionStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"DescriptionStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","description_width":""}},"41ffde370d194a588a727633cb379a33":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"DescriptionStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"DescriptionStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","description_width":""}},"6a05c68e1971435b841f83e544593819":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"DescriptionStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"DescriptionStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","description_width":""}},"6aa7eed6a9d842d38cf0fa84b045d436":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HTMLModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HTMLModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HTMLView","description":"","description_tooltip":null,"layout":"IPY_MODEL_0238f920abea48309804c8de0b48eb1a","placeholder":"​","style":"IPY_MODEL_6a05c68e1971435b841f83e544593819","value":" 1/1 [00:02&lt;00:00,  2.80s/it]"}},"7356325c4b1345998b8624f16dd4076a":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"DescriptionStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"DescriptionStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","description_width":""}},"759f3759ec9a411e830b1ffc0e05fcc1":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HTMLModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HTMLModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HTMLView","description":"","description_tooltip":null,"layout":"IPY_MODEL_d3f0745206f642d4af2fd6347052ce94","placeholder":"​","style":"IPY_MODEL_16b39aad2d854563988feec0dd4f1769","value":" 1/1 [00:09&lt;00:00,  9.15s/it]"}},"787becb0bd6348d8a04e2ef198c4e410":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"832eff621c0f4e98a6327e111c567b95":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"FloatProgressModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"FloatProgressModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"ProgressView","bar_style":"success","description":"","description_tooltip":null,"layout":"IPY_MODEL_87a9baf408784d6fa7303e586ef1e531","max":1,"min":0,"orientation":"horizontal","style":"IPY_MODEL_00b8162669e7414e90a014cfca8f0112","value":1}},"87a9baf408784d6fa7303e586ef1e531":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"9bd136706b75460ea87a95dd4f1e3e5b":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"ProgressStyleModel","state":{"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"ProgressStyleModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"StyleView","bar_color":null,"description_width":""}},"a7fd359050f7481e93b84f9753f2d428":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"ac9bdd6d2ca44fc59b82441540e834d4":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HTMLModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HTMLModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HTMLView","description":"","description_tooltip":null,"layout":"IPY_MODEL_f7cff3a03c6d44389a8f81adc41df5c0","placeholder":"​","style":"IPY_MODEL_41ffde370d194a588a727633cb379a33","value":"100%"}},"b39d35cc66364c529b8bf64930fc20dc":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"b3ac0fb979484bf5acdd000c0cc46bad":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HBoxModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HBoxModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HBoxView","box_style":"","children":["IPY_MODEL_041a24d2e34241f68c5c8ef8a32b0243","IPY_MODEL_053831b861454c6899ad5e69465a6e1f","IPY_MODEL_6aa7eed6a9d842d38cf0fa84b045d436"],"layout":"IPY_MODEL_b39d35cc66364c529b8bf64930fc20dc"}},"b9340b8f100d4aba9591a8c2b989404e":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"d21b5333ed264cb9be558a80bc1f9427":{"model_module":"@jupyter-widgets/controls","model_module_version":"1.5.0","model_name":"HBoxModel","state":{"_dom_classes":[],"_model_module":"@jupyter-widgets/controls","_model_module_version":"1.5.0","_model_name":"HBoxModel","_view_count":null,"_view_module":"@jupyter-widgets/controls","_view_module_version":"1.5.0","_view_name":"HBoxView","box_style":"","children":["IPY_MODEL_ac9bdd6d2ca44fc59b82441540e834d4","IPY_MODEL_832eff621c0f4e98a6327e111c567b95","IPY_MODEL_759f3759ec9a411e830b1ffc0e05fcc1"],"layout":"IPY_MODEL_a7fd359050f7481e93b84f9753f2d428"}},"d3f0745206f642d4af2fd6347052ce94":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}},"f7cff3a03c6d44389a8f81adc41df5c0":{"model_module":"@jupyter-widgets/base","model_module_version":"1.2.0","model_name":"LayoutModel","state":{"_model_module":"@jupyter-widgets/base","_model_module_version":"1.2.0","_model_name":"LayoutModel","_view_count":null,"_view_module":"@jupyter-widgets/base","_view_module_version":"1.2.0","_view_name":"LayoutView","align_content":null,"align_items":null,"align_self":null,"border":null,"bottom":null,"display":null,"flex":null,"flex_flow":null,"grid_area":null,"grid_auto_columns":null,"grid_auto_flow":null,"grid_auto_rows":null,"grid_column":null,"grid_gap":null,"grid_row":null,"grid_template_areas":null,"grid_template_columns":null,"grid_template_rows":null,"height":null,"justify_content":null,"justify_items":null,"left":null,"margin":null,"max_height":null,"max_width":null,"min_height":null,"min_width":null,"object_fit":null,"object_position":null,"order":null,"overflow":null,"overflow_x":null,"overflow_y":null,"padding":null,"right":null,"top":null,"visibility":null,"width":null}}},"version_major":2,"version_minor":0}}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**LIBRARIES**","metadata":{"papermill":{"duration":0.013097,"end_time":"2024-10-31T14:47:27.404853","exception":false,"start_time":"2024-10-31T14:47:27.391756","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport scipy\nfrom sklearn.decomposition import PCA\nfrom numpy.polynomial import Polynomial\nfrom tqdm.contrib.concurrent import process_map  ","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2024-11-05T18:53:46.086838Z","iopub.execute_input":"2024-11-05T18:53:46.087281Z","iopub.status.idle":"2024-11-05T18:53:46.093296Z","shell.execute_reply.started":"2024-11-05T18:53:46.087239Z","shell.execute_reply":"2024-11-05T18:53:46.092102Z"},"papermill":{"duration":4.166752,"end_time":"2024-10-31T14:47:31.584136","exception":false,"start_time":"2024-10-31T14:47:27.417384","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Configuration","metadata":{"papermill":{"duration":0.011679,"end_time":"2024-10-31T14:47:31.607842","exception":false,"start_time":"2024-10-31T14:47:31.596163","status":"completed"},"tags":[]}},{"cell_type":"code","source":"BINNING = 12 # Binning considered for AIRS signal (x12 for FGS1)\nSIGMA_TRANSITIONS = 60 # Sigma considered to construct Gaussian derivatives for transitions detection","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:46.108772Z","iopub.execute_input":"2024-11-05T18:53:46.109197Z","iopub.status.idle":"2024-11-05T18:53:46.114378Z","shell.execute_reply.started":"2024-11-05T18:53:46.109157Z","shell.execute_reply":"2024-11-05T18:53:46.113055Z"},"papermill":{"duration":0.025152,"end_time":"2024-10-31T14:47:31.644705","exception":false,"start_time":"2024-10-31T14:47:31.619553","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"MODE = 'test' # or 'train'\n\nINPUT_PATH = \"/kaggle/input/ariel-data-challenge-2024/\"\nRAW_DATA_PATH = os.path.join(INPUT_PATH, MODE)\n\nadc_info = pd.read_csv(os.path.join(INPUT_PATH, f'{MODE}_adc_info.csv')).set_index('planet_id')\naxis_info = pd.read_parquet(os.path.join(INPUT_PATH,'axis_info.parquet'))\nwavelengths = np.array(pd.read_csv(os.path.join(INPUT_PATH,'wavelengths.csv')))[0]\n\nWAVELENGTH_SELECTION =  np.array([0]+list(range(36, 318))) # relative to data generated by load_calibrated_data function\nNB_V = len(WAVELENGTH_SELECTION) # 283\n\nPLANET_NAMES = list(adc_info.index)\nSTARS = np.array(adc_info[\"star\"]).astype(\"i4\")\nSTAR_KEYS = np.unique(STARS)","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:46.134951Z","iopub.execute_input":"2024-11-05T18:53:46.1354Z","iopub.status.idle":"2024-11-05T18:53:46.356Z","shell.execute_reply.started":"2024-11-05T18:53:46.135358Z","shell.execute_reply":"2024-11-05T18:53:46.354509Z"},"papermill":{"duration":0.173358,"end_time":"2024-10-31T14:47:31.835185","exception":false,"start_time":"2024-10-31T14:47:31.661827","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data preprocessing","metadata":{}},{"cell_type":"code","source":"##################### Sensors characteristics #######################\n\nFGS1 = 'FGS1'\nAIRS = 'AIRS-CH0'\nSENSORS = [FGS1, AIRS]\n\nn_timesteps = { \n    FGS1: axis_info[FGS1 + '-axis0-h'].dropna().values.shape[0], # 135000\n    AIRS: axis_info[AIRS + '-axis0-h'].dropna().values.shape[0]  # 11250\n}\nsensor_tws_sizes = {\n    # [timesteps, wavelengths, spatial_size]\n    FGS1: [n_timesteps[FGS1], 1, 1024],\n    AIRS: [n_timesteps[AIRS], 356, 32]\n}\n\n# dts: total integration time for each non destructive read (for dark correction)\n# |--ground--|--NDR0--|--wait--|--NDR1--|--reset--|\n# 0.0       0.1       0.2      4.7      4.8       4.9    \ndts = { \n    FGS1: np.ones(n_timesteps[FGS1]) * 0.1, # 0.1 0.1 0.1 0.1 ...\n    AIRS: axis_info[AIRS + '-integration_time'].dropna().values # 0.1 4.5 0.1 4.5 ...\n}\nfor sensor in SENSORS: dts[sensor][1::2] += (0.1 + dts[sensor][0::2]) # TBC if needed: this adds (ground + NDR0) to wait time\nfor sensor in SENSORS: dts[sensor] += 0.1 # TBC if needed: this adds readout times\n\n##################### Raw data loading #######################\n\ndef read_data(planet_id, sensor):\n    \"\"\" Read the data for a given planet and sensor \"\"\"\n    n_t, n_l, n_s = sensor_tws_sizes[sensor]\n    path = f'{RAW_DATA_PATH}/{planet_id}'\n    signal = pd.read_parquet(f'{path}/{sensor}_signal.parquet').values.astype(np.float64).reshape(n_t, n_s, -1).transpose(2,0,1)\n    flat = pd.read_parquet(f'{path}/{sensor}_calibration/flat.parquet').values.astype(np.float64).reshape(n_s, -1).transpose(1,0)\n    dark = pd.read_parquet(f'{path}/{sensor}_calibration/dark.parquet').values.astype(np.float64).reshape(n_s, -1).transpose(1,0)\n    read = pd.read_parquet(f'{path}/{sensor}_calibration/read.parquet').values.astype(np.float64).reshape(n_s, -1).transpose(1,0)\n    dead = pd.read_parquet(f'{path}/{sensor}_calibration/dead.parquet').values.astype(bool).reshape(n_s, -1).transpose(1,0)\n    linear_corr = pd.read_parquet(f'{path}/{sensor}_calibration/linear_corr.parquet').values.astype(np.float64).reshape((6, n_s, -1)).transpose(2,0,1)\n    gain = adc_info[f'{sensor}_adc_gain'].loc[planet_id]\n    offset = adc_info[f'{sensor}_adc_offset'].loc[planet_id]\n    \n    return {\n        'signal': signal,\n        'flat': flat,\n        'dark': dark,\n        'read': read,\n        'dead': dead,\n        'linear_corr': linear_corr,\n        'gain': gain,\n        'offset': offset,\n        'dts': dts[sensor]\n    }\n\n##################### Calibration steps #######################\n\ndef ADC_convert(signal, gain, offset):\n    return signal / gain + offset\n\ndef mask_dead(mask, dead):\n    \"\"\" return new mask \"\"\"\n    for i in np.where(dead)[0]: \n        mask = np.delete(mask, np.where(mask == i))\n    return mask\n\ndef apply_linear_corr(linear_corr, signal):\n    for x in range(signal.shape[1]):\n        signal[:, x] = Polynomial(linear_corr[:, x])(signal[:, x])\n    return signal\n    \ndef clean_dark(signal, dark, dts):\n    dark = np.tile(dark, (signal.shape[0], 1))\n    signal -= dark * dts[:, np.newaxis]\n    return signal\n\ndef clean_read(signal, read):\n    read = np.tile(read * 2, (signal.shape[0], 1))\n    signal -= read\n    return signal\n\ndef correct_flat_field(signal, flat):\n    flat = np.tile(flat, (signal.shape[0], 1))\n    signal = signal / flat\n    return signal\n\n# if we want to ignore certain pixels\ndef relevant_pixels(data, i_wl, pct = 10):\n    \"\"\" return mask of relevant pixels \"\"\"\n    signal = data['signal'][i_wl]\n        \n    N = 10 # only evaluate one every 10 samples \n    scale = signal.shape[0] / (135000)\n    ref_signal = signal[int(10000 * scale) // 2 * 2:int(15000 * scale) // 2 * 2] # only evaluate part of graph immediately before transit\n    simple_cds = ref_signal[1::2*N] - ref_signal[0::2*N]\n\n    simple_cds = simple_cds.sum(axis = 0) / data['flat'][i_wl]\n    simple_cds_idx = np.flip(np.argsort(simple_cds))\n    n = int(simple_cds.shape[0] * pct / 100.)\n    mask = simple_cds_idx[:n]\n    \n    return mask_dead(mask, data['dead'][i_wl])\n\ndef calibrate(data, i_wl, mask):\n    signal = data['signal'][i_wl]\n    m_signal = ADC_convert(signal[:, mask], data['gain'], data['offset'])\n    \n    np.clip(m_signal, 0., None, m_signal)\n    m_signal = apply_linear_corr(data['linear_corr'][i_wl][:, mask], m_signal)\n    \n    m_signal = clean_dark(m_signal, data['dark'][i_wl][mask], data['dts'])\n    m_signal = clean_read(m_signal, data['read'][i_wl][mask])\n    m_signal = correct_flat_field(m_signal, data['flat'][i_wl][mask])\n    \n    mean_signal = m_signal.mean(axis=1)\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n\n    return net_signal\n\n########################### Main calibration function #######################\n\ndef load_calibrated_data(planet_name, f_signal_pct = 100, a_signal_pct = 100):\n    \"\"\"\n    Prepare or load calibrated data\n    FGS1 signal binned x12 to be similar to AIRS signals\n    Wavelength in returned array: FGS1 then all 356 AIRS in reverse order\n    \n    planet_name: name the planet\n    f_signal_pct: [0, 100] only select a subset (in percent) of the most illuminated pixels for FGS1\n    a_signal_pct: [0, 100] only select a subset (in percent) of the most illuminated pixels for AIRS\n    \"\"\"\n    # load data for planet\n    f_data = read_data(planet_name, 'FGS1')\n    a_data = read_data(planet_name, 'AIRS-CH0')\n\n    # Calibrations \n    # * FGS1 signal\n    f_mask = relevant_pixels(f_data, 0, pct = f_signal_pct)\n    f_signal = calibrate(f_data, 0, f_mask)\n    f_data = None # free memory?\n    f_signal_12 = f_signal.reshape(-1, 12).mean(axis=1) # 135000/2 => 11250/2\n    # * AIRS signals\n    a_signals = []    \n    for i_wl in range(356):\n        a_mask = relevant_pixels(a_data, i_wl, pct = a_signal_pct)\n        a_signal = calibrate(a_data, i_wl, a_mask)\n        a_signals.append(a_signal)\n    a_data = None # free memory?\n    a_signals.append(f_signal_12)\n    a_signals = np.array(a_signals)[::-1] # /!\\ wavelengths are in reversed order!\n    \n    return a_signals\n\n########################### Retrieve dead pixels #######################\n\ndef get_dead_pixels(planet_name):\n    path = \"%s/%d/AIRS-CH0_calibration/dead.parquet\" % (RAW_DATA_PATH, planet_name)\n    dead_table = pd.read_parquet(path).values.astype(\"?\").reshape((32, 356))\n    deads = np.where(dead_table)\n    dead_pixels = []\n    for i in range(len(deads[0])):\n        v = 356-deads[1][i]\n        p = deads[0][i]    \n        if 8 <= p < 24:\n            dead_pixels.append((v, p))\n        \n    return dead_pixels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-05T18:53:46.358516Z","iopub.execute_input":"2024-11-05T18:53:46.35904Z","iopub.status.idle":"2024-11-05T18:53:46.397122Z","shell.execute_reply.started":"2024-11-05T18:53:46.358971Z","shell.execute_reply":"2024-11-05T18:53:46.395848Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Low-level timeseries processing","metadata":{"papermill":{"duration":0.011636,"end_time":"2024-10-31T14:47:32.015515","exception":false,"start_time":"2024-10-31T14:47:32.003879","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def clip_outliers(data, sigma, n = 15):\n    \"\"\" sliding in-place removal of outliers\"\"\"\n\n    n = min(n, data.size)\n    n -= n%2\n    off = n//2\n\n    cumsum = np.cumsum(data)\n    cumsum2 = np.cumsum(np.square(data))\n      \n    mean_win = np.empty(data.shape, \"f8\")\n    std_win = np.empty(data.shape, \"f8\")        \n\n    mean_win[off:-off] = (cumsum[n:]-cumsum[:-n])/n\n    std_win[off:-off] = np.sqrt((cumsum2[n:]-cumsum2[:-n])/n - np.square(mean_win[off:-off]))\n\n    mean_win[:off] = mean_win[off]\n    mean_win[-off:] = mean_win[-off-1]    \n    std_win[:off] = std_win[off]\n    std_win[-off:] = std_win[-off-1]   \n    \n    data = np.clip(data, mean_win-std_win*sigma, mean_win+std_win*sigma)\n    \n    return data\n\n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\ndef safe_savgol_filter(data, window_size, order=2):\n    \"\"\" Data smoothing via savgol_filter \"\"\"\n    if data.size < window_size:\n        return data\n    return scipy.signal.savgol_filter(data, window_size, order)\n\n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\ndef binning(timeseries, bin_size):\n    \"\"\"reduce timeseries size along first axis of a factor bin_size\"\"\"\n    new_size = timeseries.shape[1]//bin_size # incomplete last bit is dropped\n    limit = bin_size*new_size\n    bin_sum = timeseries[:, 0:limit:bin_size].copy()\n    bin_sum2 = np.square(bin_sum)\n    for i in range(1, bin_size):\n        bin_sum[:, :] += timeseries[:, i:limit:bin_size]\n        bin_sum2[:, :] += np.square(timeseries[:, i:limit:bin_size])\n    bin_mean = bin_sum/bin_size\n    bin_var = (bin_sum2-np.square(bin_mean))/bin_size    \n    return bin_mean, bin_var","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:46.398707Z","iopub.execute_input":"2024-11-05T18:53:46.400197Z","iopub.status.idle":"2024-11-05T18:53:46.4135Z","shell.execute_reply.started":"2024-11-05T18:53:46.400141Z","shell.execute_reply":"2024-11-05T18:53:46.412052Z"},"papermill":{"duration":9.472374,"end_time":"2024-10-31T14:47:41.499714","exception":false,"start_time":"2024-10-31T14:47:32.02734","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data pre-loading","metadata":{}},{"cell_type":"code","source":"PATH_CACHE_CALIB = '/tmp/cache'\nos.makedirs(PATH_CACHE_CALIB, exist_ok=True)\n\ndef get_planet_data(planet_name, bin_size=1): \n    path_planet = \"%s/%d_calibrated_signal.npy\" % (PATH_CACHE_CALIB, planet_name)\n    if os.path.isfile(path_planet):\n        data = np.load(path_planet)\n    else:\n        data = load_calibrated_data(planet_name, f_signal_pct = 50, a_signal_pct = 50)\n        np.save(path_planet, data)\n    \n    if bin_size > 1:\n        bin_mean, bin_var = binning(data, bin_size)\n        return bin_mean\n    else:\n        return data\n\n# preload data\ndef preload(p):\n    get_planet_data(p)\n\nprocess_map(preload, PLANET_NAMES, max_workers=4)","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:46.41591Z","iopub.execute_input":"2024-11-05T18:53:46.416695Z","iopub.status.idle":"2024-11-05T18:53:56.347446Z","shell.execute_reply.started":"2024-11-05T18:53:46.416641Z","shell.execute_reply":"2024-11-05T18:53:56.346011Z"},"papermill":{"duration":9.472374,"end_time":"2024-10-31T14:47:41.499714","exception":false,"start_time":"2024-10-31T14:47:32.02734","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Ingress and egress transitions","metadata":{"papermill":{"duration":0.015652,"end_time":"2024-10-31T14:47:31.862944","exception":false,"start_time":"2024-10-31T14:47:31.847292","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# gaussian derivatives\ndef dgauss(sig):\n    xs = np.arange(-3.*sig, 3.*sig+1)\n    den = 2.*sig*sig\n    ys = np.exp(-np.square(xs)/den)\n    dys = -2*xs/den*ys\n    return dys\n\ndef d2gauss(sig):\n    xs = np.arange(-3.*sig, 3.*sig+1)\n    den = 2.*sig*sig\n    ys = np.exp(-np.square(xs)/den)\n    d2ys = np.square(2/den)*ys*(xs-sig)*(xs+sig)\n    return d2ys","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:56.349957Z","iopub.execute_input":"2024-11-05T18:53:56.350451Z","iopub.status.idle":"2024-11-05T18:53:56.359972Z","shell.execute_reply.started":"2024-11-05T18:53:56.350398Z","shell.execute_reply":"2024-11-05T18:53:56.35871Z"},"papermill":{"duration":0.030978,"end_time":"2024-10-31T14:47:31.953487","exception":false,"start_time":"2024-10-31T14:47:31.922509","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def find_transit_edges(S, sigma):\n    \"\"\" Find the centers of the transitions \"\"\"\n\n    Sc = np.convolve(S, dgauss(sigma), mode=\"valid\")\n    off = int((S.size-Sc.size)/2)\n    mid = Sc.size//2\n    \n    transit_start = np.argmin(Sc[3:mid-3])+off+3\n    transit_end = np.argmax(Sc[mid+3:-3])+off+mid+3\n\n    return transit_start, transit_end\n\n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\ndef find_transit_slopes(S, transit_start, transit_end, sigma):\n    \"\"\"find the width of the transitions\"\"\"\n    \n    Sc2 = np.convolve(S, d2gauss(sigma), mode=\"valid\")\n    off = int((S.size-Sc2.size)/2)\n\n    t1 = transit_start - off\n    t2 = transit_end - off\n    \n    sz = 2*sigma\n    t1a = np.argmin(Sc2[t1-sz:t1+1])+t1-sz+off\n    t1b = np.argmax(Sc2[t1:t1+sz+1])+t1+off\n    t2a = np.argmax(Sc2[t2-sz:t2+1])+t2-sz+off\n    t2b = np.argmin(Sc2[t2:t2+sz+1])+t2+off\n\n    return t1a, t1b, t2a, t2b","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:56.361688Z","iopub.execute_input":"2024-11-05T18:53:56.362189Z","iopub.status.idle":"2024-11-05T18:53:56.376905Z","shell.execute_reply.started":"2024-11-05T18:53:56.362138Z","shell.execute_reply":"2024-11-05T18:53:56.375875Z"},"papermill":{"duration":0.025857,"end_time":"2024-10-31T14:47:31.991415","exception":false,"start_time":"2024-10-31T14:47:31.965558","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## SNR and weights","metadata":{}},{"cell_type":"code","source":"# Compute SNR and weights\ndef star_statistics():\n    \"\"\" SNR statistics by star and weights per planet \"\"\"\n    \n    all_means = np.zeros((len(PLANET_NAMES), 357), \"f8\")\n    all_vars  = np.zeros((len(PLANET_NAMES), 357), \"f8\")\n    default_transition_width = 150\n    \n    for k,p in enumerate(PLANET_NAMES):    \n        data = get_planet_data(p)\n        S = np.mean(data, axis=0)\n        S /= S.mean()\n        start, end = find_transit_edges(S, SIGMA_TRANSITIONS)\n        outs1 = data[:, :start-default_transition_width]\n        outs2 = data[:, end+default_transition_width+1:]\n        all_means[k, :] = (outs1.mean(axis=1) + outs2.mean(axis=1)) / 2.\n        all_vars [k, :] = (outs1.var (axis=1) + outs2.var (axis=1)) / 2.\n    \n    ideal_weights = all_means / all_vars # at same noise level (divide by std), weight by snr\n    snr = all_means / np.sqrt(all_vars)\n        \n    star_snrs = {}\n    for star in STAR_KEYS:\n        star_snrs[star] = snr[STARS==star].mean(axis=0)\n        \n    return ideal_weights, star_snrs\n\nPLANETS_WEIGHTS, STARS_SNRS = star_statistics()","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:53:56.379288Z","iopub.execute_input":"2024-11-05T18:53:56.379716Z","iopub.status.idle":"2024-11-05T18:54:03.761777Z","shell.execute_reply.started":"2024-11-05T18:53:56.379673Z","shell.execute_reply":"2024-11-05T18:54:03.760441Z"},"papermill":{"duration":0.058988,"end_time":"2024-10-31T14:47:41.571062","exception":false,"start_time":"2024-10-31T14:47:41.512074","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Correlation with true spectra","metadata":{"papermill":{"duration":0.012364,"end_time":"2024-10-31T14:47:41.602507","exception":false,"start_time":"2024-10-31T14:47:41.590143","status":"completed"},"tags":[]}},{"cell_type":"code","source":"NN = 160 # correlation is only evaluated on the first wavelengths where the SNR is high\n\nrandom_labels = np.load('/kaggle/input/adc-2024-space-coders-labels-generation/labels.npy')\ntrain_labels = pd.read_csv(\"/kaggle/input/adc-2024-space-coders-labels-generation/train_labels.csv\").set_index('planet_id')\nall_labels = np.vstack((train_labels, random_labels))\n\nNORMALIZED_LABELS = (all_labels.transpose()-all_labels[:,:NN].mean(axis=1)) / all_labels[:,:NN].std(axis=1)\nNORMALIZED_LABELS = NORMALIZED_LABELS.transpose()\n\ndef find_best_correllation(sp, star):\n    \"\"\" Correllation with some true labels \"\"\"\n    sp_ = (sp - sp[:NN].mean()) / sp[:NN].std()\n    corrs = [np.corrcoef(sp[:NN], NORMALIZED_LABELS[i, :NN])[0, 1] for i in range(len(NORMALIZED_LABELS))]\n\n    return max(corrs)","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:54:03.763209Z","iopub.execute_input":"2024-11-05T18:54:03.763588Z","iopub.status.idle":"2024-11-05T18:54:04.130489Z","shell.execute_reply.started":"2024-11-05T18:54:03.76355Z","shell.execute_reply":"2024-11-05T18:54:04.129229Z"},"papermill":{"duration":0.302988,"end_time":"2024-10-31T14:47:41.964679","exception":false,"start_time":"2024-10-31T14:47:41.661691","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Raw spectrum extraction","metadata":{"papermill":{"duration":0.012959,"end_time":"2024-10-31T14:47:41.994813","exception":false,"start_time":"2024-10-31T14:47:41.981854","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def compute_weights(size, t1a, t1b, t2a, t2b):\n    \"\"\" Compute fitting weights (ignore transitions) \"\"\"\n    ws = np.ones(size, \"f8\")\n    ws[t1a-1:t1b+2] = 1.e-5\n    ws[t2a-1:t2b+2] = 1.e-5\n    return ws\n    \n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nCACHE_POLY = {}\ndef estimate_transit_depth(signal0, signalk, start, end, reset_cache=False, deg0=None):\n    \"\"\" Estimate the transit depth for a wavelength \"\"\"\n\n    global CACHE_POLY # keep fit params from last wavelength to facilitate convergence\n    if reset_cache: CACHE_POLY = {}\n\n    # 1) Wide detrending\n    xdata = np.arange(len(signal0))\n    t1a, t1b, t2a, t2b = find_transit_slopes(signal0, start, end, SIGMA_TRANSITIONS//BINNING)\n    S = signal0/signal0[(xdata<t1a)|(xdata>t2b)].mean()\n    signal0 = S.copy()\n\n    SMOOTHING_W = 1400//BINNING # window used for signal smoothing before curve_fit\n    for left, right in [(0, t1a+1), (t1b, t2a+1), (t2b, len(xdata))]:\n        S[left:right] = safe_savgol_filter(clip_outliers(S[left:right], 1., 15), SMOOTHING_W)\n    \n    weights = compute_weights(S.size, t1a, t1b, t2a, t2b)      \n    transit_mask = np.zeros_like(S)\n    transit_mask[start:end+1] = 1.\n    \n    def fit_fn(x, da, *P): return Polynomial(P)(x)* (1.-transit_mask*da) \n\n    # Evaluate best degree\n    degs = [2, 4, 5]\n    if not deg0 is None:\n        degs = [deg for deg in degs if deg <= deg0]\n    best_fit = (1e10, 0., 0., [])\n    for deg in degs:\n     \n        if deg in CACHE_POLY:\n            init_input = CACHE_POLY[deg]\n        else:\n            init_input = [1e-3, 1.0]+[0.]*(deg) \n            \n        popt, pcov = scipy.optimize.curve_fit(fit_fn, xdata, S, sigma=1./weights, p0=init_input)\n        CACHE_POLY[deg] = popt # save fit params for next wavelength\n        pcorr = fit_fn(xdata, *popt)\n        pcorr0 = fit_fn(xdata, 0, *popt[1:])   \n        \n        rms = np.sqrt(np.average(np.square(pcorr-S), weights=weights))\n        rms_wgt = deg*rms\n        poly = popt[1:]\n      \n        area_ratio = 1.0 - np.mean(S[t1b+2:t2a-1]/pcorr0[t1b+2:t2a-1])\n        best_fit = min(best_fit, (rms_wgt, area_ratio, deg, poly))\n \n    (rms_wgt, best_ratio, deg, poly) = best_fit\n    rms = rms_wgt/deg\n\n    # 2) Narrow adjustment\n    Sk = None\n    if signalk is not None:\n        \n        t1a, t1b, t2a, t2b = find_transit_slopes(signalk, start, end, SIGMA_TRANSITIONS//BINNING)\n        Sk = signalk/signalk[(xdata<t1a)|(xdata>t2b)].mean()\n        signalk = Sk.copy()\n\n        for left, right in [(0, t1a+1), (t1b, t2a+1), (t2b, len(xdata))]:\n            Sk[left:right] = safe_savgol_filter(clip_outliers(Sk[left:right], 1., 15), SMOOTHING_W)\n        \n        def fit_fn(x, da, a, b): \n            return (a + b * Polynomial(poly)(x))* (1.-transit_mask*da) \n    \n        init_input = [best_ratio, 0.0, 1.]            \n        popt, pcov = scipy.optimize.curve_fit(fit_fn, xdata, Sk, sigma=1./weights, p0=init_input)\n        pcorr = fit_fn(xdata, *popt)\n        pcorr0 = fit_fn(xdata, 0, *popt[1:])   \n        \n        rms = np.sqrt(np.average(np.square(pcorr-Sk), weights=weights))\n        best_ratio = 1.0-np.mean(Sk[t1b+2:t2a-1]/pcorr0[t1b+2:t2a-1])\n        \n    return best_ratio, rms, deg\n\n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\ndef build_spectrum(p):\n    \"\"\"build the raw spectrum\"\"\"\n    \n    # get weighted calibrated data\n    star = STARS[p]\n    data = get_planet_data(PLANET_NAMES[p], BINNING)\n    data = (data.transpose()*PLANETS_WEIGHTS[p]).transpose()\n    \n    # lower weights of dead pixels\n    for (v, x) in get_dead_pixels(PLANET_NAMES[p]):\n        if x in [15, 16]:\n            data[v, :] *= 0.01\n        elif x in [14, 17]:\n            data[v, :] *= 0.1 \n        elif x in [13, 18]:\n            data[v, :] *= 0.5     \n    \n    # overall signal to determine best max degree\n    Sall = np.mean(data[WAVELENGTH_SELECTION[1:150]], axis=0)\n    start, end = find_transit_edges(Sall, SIGMA_TRANSITIONS//BINNING)\n    _, _, deg_avg = estimate_transit_depth(Sall, None, start, end, True)\n                           \n    # prepare output\n    spectrum = np.zeros((NB_V), \"f8\")\n    rmss = np.zeros((NB_V), \"f8\")\n    x = np.arange((data.shape[1]))    \n\n    # Loop on wavelengths\n    \n    # FGS1\n    S = data[WAVELENGTH_SELECTION[0]]\n    spectrum[0], rmss[0], _ = estimate_transit_depth(S, None, start, end) \n\n    # AIRS-CH0\n    for k, v in enumerate(WAVELENGTH_SELECTION[1:]): \n        # wide window\n        N_wl = 100  \n        left = max(1, v-N_wl)\n        right = min(357, v+N_wl+1)\n        S = np.mean(data[left:right, :], axis=0)\n        S /= S.mean()\n\n        # narrow window\n        if k < 200: Nk = 8\n        else: Nk = 20\n        if star <= 1 and k >= 180 and k <= 186: Nk = 5 # reduce Nk for CH4 peak\n        \n        left_k = v-Nk\n        right_k = v+Nk+1\n        Sk = np.mean(data[left_k:right_k, :], axis=0)\n        Sk /= Sk.mean()\n\n        spectrum[k+1], rmss[k+1], _ = estimate_transit_depth(S, Sk, start, end, deg0=deg_avg)\n\n    return np.concatenate([spectrum, rmss])","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:54:04.132077Z","iopub.execute_input":"2024-11-05T18:54:04.132438Z","iopub.status.idle":"2024-11-05T18:54:04.159797Z","shell.execute_reply.started":"2024-11-05T18:54:04.1324Z","shell.execute_reply":"2024-11-05T18:54:04.15854Z"},"papermill":{"duration":0.09207,"end_time":"2024-10-31T14:47:42.146229","exception":false,"start_time":"2024-10-31T14:47:42.054159","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Main loop","metadata":{"papermill":{"duration":0.012677,"end_time":"2024-10-31T14:47:42.171701","exception":false,"start_time":"2024-10-31T14:47:42.159024","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def process_planet(p):\n    \"\"\"process fully one planet\"\"\"\n\n    # Compute raw spectrum\n    raw_spectrum = build_spectrum(p)\n\n    # Add some offset\n    # we should have estimated foreground here instead :-)\n    coeff = np.minimum(1 + 0.128 * np.sqrt(np.clip(raw_spectrum[:NB_V], 2e-4, 0.9)), 1.009) \n    spectrum = np.clip(coeff*raw_spectrum[:NB_V], 2e-4, 0.9)\n\n    # Estimate correlation of the raw spectrum with some possible labels\n    star = STARS[p]\n    true_corr = find_best_correllation(spectrum, star)\n\n    # Smooth the spectrum\n    if star <= 1:\n        spectrum[:180] = safe_savgol_filter(spectrum[:180], 15)\n        # none between 180 and 186 (CH4 peak)\n        spectrum[187:] = safe_savgol_filter(spectrum[187:], 15)\n    else:\n        spectrum = safe_savgol_filter(spectrum, 15)\n    \n    # Evaluation of the spectrum dynamics, based on first wavelengths\n    spectrum_head = spectrum[:190+1]\n    begin_sp_range = np.max(spectrum_head) - np.min(spectrum_head)\n    begin_sp_mean = spectrum_head.mean() \n    begin_sp_std = spectrum_head.std() \n    low_dyn = true_corr < (0.89 - 200 * begin_sp_range)\n\n    # coefficients for computation of prediction and sigma (as barycenters between \"avg\" and \"standard\" values)\n    dyncoeff_pred = 1.\n    dyncoeff_sigma = 1.\n    if low_dyn:\n        dyncoeff_pred = 0.15 + max(true_corr, 0.) * (0.6 if star <= 2 else 0.4)\n        dyncoeff_sigma = 0.165\n    \n    # trailing ramp\n    if low_dyn or star > 0:\n        TRAILING_RAMP_IND = 200+1\n        x_ramp = np.arange(NB_V-TRAILING_RAMP_IND)\n        linepoly = np.polyfit(x_ramp, spectrum[TRAILING_RAMP_IND:], 1)\n        line = np.poly1d(linepoly)\n        spectrum[TRAILING_RAMP_IND:] = line(x_ramp)\n    \n    # Limit raw spectrum enveloppe\n    if low_dyn or star > 2:\n        spectrum = np.maximum(spectrum, begin_sp_mean - 2 * begin_sp_std)       \n        spectrum[45:-45] = np.maximum(spectrum[45:-45], begin_sp_mean - 1 * begin_sp_std)\n        if star <= 1:\n            spectrum[80:-80] = np.maximum(spectrum[80:-80], begin_sp_mean - 0.8 * begin_sp_std)\n    \n        spectrum[:84] = np.minimum(spectrum[:84], begin_sp_mean + 1. * begin_sp_std)  \n        spectrum[84:124] = np.minimum(spectrum[84:124], begin_sp_mean + 2 * begin_sp_std)\n        spectrum[124:150] = np.minimum(spectrum[124:150], begin_sp_mean + 1.5 * begin_sp_std)\n        spectrum[150:176] = np.minimum(spectrum[150:176], begin_sp_mean + 2 * begin_sp_std)\n        spectrum[192:] = np.minimum(spectrum[192:], begin_sp_mean + 1.5 * begin_sp_std)  \n        spectrum[210:] = np.minimum(spectrum[210:], begin_sp_mean + 1 * begin_sp_std)  \n\n\n    # flattening if low dynamic\n    mean = np.average(spectrum, weights=STARS_SNRS[star][WAVELENGTH_SELECTION])\n    spectrum_low =  mean*np.ones((NB_V), \"f8\")\n\n    delta_ref = np.abs(spectrum - mean)  \n        \n    # UNCERTAINTY\n    # low_dynamic uncertainty estimation\n    clip_max = 2.4e-5\n    clip_min = max(7.5e-6, 0.07*begin_sp_range)\n    factor   = 0.01 if star <= 2 else 0.015 \n    std0_low = max(min(factor*np.max(spectrum), clip_max), clip_min)\n    delta_ref_factor = 0.45 if star <= 2 else 0.5\n    stds_low = np.maximum(std0_low, delta_ref[:] * delta_ref_factor)    \n    \n    # high_dynamic uncertainty estimation                               \n    clip_max = 6.5e-5 if star <= 2 else 8.0e-5\n    clip_min = max(4.0e-5, 0.06 * begin_sp_range)\n    factor = 0.025\n    std0_high = max(min(factor*np.max(spectrum), clip_max), clip_min)                                \n    stds_high = np.maximum(std0_high, delta_ref[:] * 0.34)\n\n    \n    # RESULTS\n    results = {\"standard\" : np.concatenate([spectrum, stds_high]),\n               \"avg\" : np.concatenate([spectrum_low, stds_low])}\n             \n    # INFOS\n    infos = {\"low_dyn\" : low_dyn,\n             \"dyncoeff_pred\" : dyncoeff_pred,\n             \"dyncoeff_sigma\" : dyncoeff_sigma}  \n\n    return results, infos\n\n\n#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nalls = process_map(process_planet, range(len(PLANET_NAMES)), max_workers=4)\npack = [elem[0] for elem in alls]\ninfos = [elem[1] for elem in alls]\n\nresults = {key : np.array([result[key] for result in pack]) for key in pack[0].keys()}\nall_infos = {key : np.array([elem[key] for elem in infos]) for key in infos[0].keys()}","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:54:04.161261Z","iopub.execute_input":"2024-11-05T18:54:04.161836Z","iopub.status.idle":"2024-11-05T18:54:14.491565Z","shell.execute_reply.started":"2024-11-05T18:54:04.161786Z","shell.execute_reply":"2024-11-05T18:54:14.490309Z"},"papermill":{"duration":2.918105,"end_time":"2024-10-31T14:47:45.102212","exception":false,"start_time":"2024-10-31T14:47:42.184107","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Spectra postprocessing with PCA","metadata":{}},{"cell_type":"code","source":"def PCA_projection(spectrums, nb_pca):\n    \"\"\"projection on reduced PCA space\"\"\"\n    pca = PCA()\n    proj = pca.fit_transform(spectrums)  \n    proj[:, nb_pca:] = 0\n    return pca.inverse_transform(proj)[:, :]\n    \n# correct standard (high dynamics) spectra with PCA\nAVG_STAR_SNR = np.mean(np.array([STARS_SNRS[star][WAVELENGTH_SELECTION] for star in STAR_KEYS]), axis=0)\nraws = results[\"standard\"][:, :NB_V].copy()\nresults[\"standard\"][:, :NB_V] = PCA_projection(raws*AVG_STAR_SNR, 7)/AVG_STAR_SNR\n\n# correct\nfor star in [0, 1, 2]:\n    if not star in STAR_KEYS: continue\n    selec = STARS == star\n    coeff = STARS_SNRS[star][WAVELENGTH_SELECTION]    \n    results[\"standard\"][selec, :NB_V] = PCA_projection(raws[selec, :]*coeff, 5)/coeff","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:54:14.493322Z","iopub.execute_input":"2024-11-05T18:54:14.493717Z","iopub.status.idle":"2024-11-05T18:54:14.513382Z","shell.execute_reply.started":"2024-11-05T18:54:14.493666Z","shell.execute_reply":"2024-11-05T18:54:14.51187Z"},"papermill":{"duration":2.918105,"end_time":"2024-10-31T14:47:45.102212","exception":false,"start_time":"2024-10-31T14:47:42.184107","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Submission file finalization","metadata":{}},{"cell_type":"code","source":"# Compute true spectrum\ndef true_spectrum(standard, avg, dyncoeff_sigma, dyncoeff_pred):\n    \"\"\"compute true spectrum weighting with dyncoeff\"\"\"\n    spectrum_true = dyncoeff_pred*standard[:, :NB_V].transpose() + (1-dyncoeff_pred)*avg[:, :NB_V].transpose()\n    stds_true = np.sqrt(dyncoeff_sigma*np.square(standard[:, NB_V:].transpose()) + (1-dyncoeff_sigma)*np.square(avg[:, NB_V:].transpose()))    \n    return np.concatenate([spectrum_true.transpose(), stds_true.transpose()], axis=1)\n    \nresults[\"true\"] = true_spectrum(results[\"standard\"], results[\"avg\"], all_infos[\"dyncoeff_sigma\"], all_infos[\"dyncoeff_pred\"])\n\n# Adjust sigmas that were finally too pessimistic \nresults[\"true\"][(STARS <= 2) & (~all_infos['low_dyn']), NB_V:] *= 0.8\nresults[\"true\"][(STARS <= 2), NB_V:] *= 0.9\n\n# save submit file\nindexes = np.array([p for p in PLANET_NAMES])\ncolumns = [\"wl_%d\" % i for i in range(1, NB_V+1)] + [\"sigma_%d\" % i for i in range(1, NB_V+1)]\ndf = pd.DataFrame(np.array(results[\"true\"]), columns=columns, index=indexes)\ndf.index.name = \"planet_id\"\ndf.to_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-11-05T18:54:14.514875Z","iopub.execute_input":"2024-11-05T18:54:14.516191Z","iopub.status.idle":"2024-11-05T18:54:14.536987Z","shell.execute_reply.started":"2024-11-05T18:54:14.516138Z","shell.execute_reply":"2024-11-05T18:54:14.53584Z"},"papermill":{"duration":2.918105,"end_time":"2024-10-31T14:47:45.102212","exception":false,"start_time":"2024-10-31T14:47:42.184107","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null}]}