{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":9110664,"sourceType":"datasetVersion","datasetId":5498833},{"sourceId":9177563,"sourceType":"datasetVersion","datasetId":5546655}],"dockerImageVersionId":30747,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# ADC 2024 starter notebook ","metadata":{}},{"cell_type":"markdown","source":"This baseline notebook is designed to offer a starting point for the competiters. Please note that the approach we've taken is not *THE* solution — it's simply ONE possible approach. Our aim is to assist participants in exploring different ways to preprocess and model the data. Please feel free to fork the notebook and save the model/data for your own exploration.\n\n","metadata":{}},{"cell_type":"markdown","source":"This notebook was prepared by Virginie Batista and Angèle Syty from the Institut d'Astrophysique de Paris, and Orphée Faucoz from Centre National d’Etudes Spatiales (CNES), with support from Gordon Yip and Tara Tahseen from University College London.","metadata":{}},{"cell_type":"markdown","source":"# READ THIS BEFORE YOU PROCEED\nThis training procedure uses the light dataset produced from this [notebook (Version 5)](https://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data). We applied all the calibration steps EXCEPT Linearity Correction with Chunksize = 1. The binned dataset is available to download [here](https://www.kaggle.com/datasets/gordonyip/binned-dataset-v3/data). *If you want to carry out all the correction, you will have to do so yourself.*\n\n\n**This notebook will only provide the model checkpoints, you are welcomed to use these checkpoints with your own script and submit to the leaderboard.** ","metadata":{}},{"cell_type":"markdown","source":"## Task overview","metadata":{}},{"cell_type":"markdown","source":"The challenge's primary objective is to process these exposures to produce a single, clean spectrum for each exoplanet, summarizing the rp/rs values across all wavelengths.\n\nThe exposure are subject to noises and the images or spectrum are not perfect. The Jitter noise has a complex signature that the ML model should recognize to produce a better spectra.\n\nDifferent techniques are possible and are up to the participant imagination to produce a novel (and hopefully better) solution to this task.\n\nHere outline our baseline approach :\n\nWe first fit  a 1D CNN to fit the mean value of the transmission spectra, taking as input the transit white curve (total flux of each image taken as a function of time).\n\nFor the second part of the baseline, to retrieve the atmopsheric features, we make the data lighter by summing up the fluxes along the y-axis, for each wavelength, resulting in 2D images of dimension (N_times, N_wavelengths). We also cut the signal to remove the out of transit in order to enhance transit depth variations between wavelengths. For the same reason, we substract the mean flux, corresponding to the average transit depth, to keep only wavelength variations around this mean. We use a 2D CNN to fit the atmospheric features. \n","metadata":{}},{"cell_type":"code","source":"import matplotlib.image as mpimg\nimport matplotlib.pyplot as plt\n\nimg = mpimg.imread('/kaggle/input/baseline-img/2nd_baseline.png')\nplt.figure(figsize=(10, 15))\nplt.imshow(img)\nplt.axis('off') \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:01.0771Z","iopub.execute_input":"2024-08-28T10:17:01.078153Z","iopub.status.idle":"2024-08-28T10:17:01.637171Z","shell.execute_reply.started":"2024-08-28T10:17:01.078118Z","shell.execute_reply":"2024-08-28T10:17:01.636235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Import library\n<a id=\"import\"></a>","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nfrom tensorflow.keras.models import load_model\nimport tensorflow as tf \nimport random \nimport os\nfrom tensorflow.keras.losses import MeanAbsoluteError\nfrom matplotlib.ticker import ScalarFormatter\nimport pandas as pd\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:01.63908Z","iopub.execute_input":"2024-08-28T10:17:01.639779Z","iopub.status.idle":"2024-08-28T10:17:06.735786Z","shell.execute_reply.started":"2024-08-28T10:17:01.639737Z","shell.execute_reply":"2024-08-28T10:17:06.73464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Setup Paths and Read Data","metadata":{}},{"cell_type":"code","source":"data_folder = '/kaggle/input/binned-dataset-v3/' # path to the folder containing the data\nauxiliary_folder = '/kaggle/input/ariel-data-challenge-2024/' # path to the folder containing the train targets and wavelengths informations","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:06.737044Z","iopub.execute_input":"2024-08-28T10:17:06.737642Z","iopub.status.idle":"2024-08-28T10:17:06.74274Z","shell.execute_reply.started":"2024-08-28T10:17:06.737609Z","shell.execute_reply":"2024-08-28T10:17:06.741434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_train = np.load(f'{data_folder}/data_train.npy')\ndata_train_FGS = np.load(f'{data_folder}/data_train_FGS.npy')\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:06.746005Z","iopub.execute_input":"2024-08-28T10:17:06.746468Z","iopub.status.idle":"2024-08-28T10:17:54.016348Z","shell.execute_reply.started":"2024-08-28T10:17:06.746432Z","shell.execute_reply":"2024-08-28T10:17:54.014985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We create a directory to save the outputs of this notebook, and define the hyperparameters of the model","metadata":{}},{"cell_type":"code","source":"output_dir = './output'\n\nSEED = 42\n\ndo_the_mcdropout_wc = True\ndo_the_mcdropout = True\n\nif not os.path.exists(output_dir):\n    os.makedirs(output_dir)\n    print(f\"Directory {output_dir} created.\")\nelse:\n    print(f\"Directory {output_dir} already exists.\")","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:54.017791Z","iopub.execute_input":"2024-08-28T10:17:54.0182Z","iopub.status.idle":"2024-08-28T10:17:54.025061Z","shell.execute_reply.started":"2024-08-28T10:17:54.018159Z","shell.execute_reply":"2024-08-28T10:17:54.024063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1D-CNN for mean transit depth","metadata":{}},{"cell_type":"markdown","source":"## Preprocessing for 1D CNN","metadata":{}},{"cell_type":"code","source":"train_solution = np.loadtxt(f'{auxiliary_folder}/train_labels.csv', delimiter = ',', skiprows = 1)\n\ntargets = train_solution[:,1:]\ntargets_mean = targets[:,1:].mean(axis = 1) # used for the 1D-CNN to extract the mean value, only AIRS wavelengths as the FGS point is not used in the white curve\nN = targets.shape[0]\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:54.026564Z","iopub.execute_input":"2024-08-28T10:17:54.027011Z","iopub.status.idle":"2024-08-28T10:17:54.081333Z","shell.execute_reply.started":"2024-08-28T10:17:54.026977Z","shell.execute_reply":"2024-08-28T10:17:54.080197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We create the dataset by adding the FGS frame, crushed in one column, at the end of the AIRS data cube.  \nThe images are normalized using the star spectrum extracted from the images themselves.","metadata":{}},{"cell_type":"code","source":"signal_AIRS_diff_transposed_binned, signal_FGS_diff_transposed_binned  = data_train, data_train_FGS\nFGS_column = signal_FGS_diff_transposed_binned.sum(axis = 2)\ndataset = np.concatenate([signal_AIRS_diff_transposed_binned, FGS_column[:,:, np.newaxis,:]], axis = 2)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:54.082729Z","iopub.execute_input":"2024-08-28T10:17:54.083058Z","iopub.status.idle":"2024-08-28T10:17:59.166878Z","shell.execute_reply.started":"2024-08-28T10:17:54.083031Z","shell.execute_reply":"2024-08-28T10:17:59.166052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"we sum up the pixels on the y-axis to transform the data into 2D images","metadata":{}},{"cell_type":"code","source":"dataset = dataset.sum(axis=3)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:17:59.16807Z","iopub.execute_input":"2024-08-28T10:17:59.168382Z","iopub.status.idle":"2024-08-28T10:18:01.065755Z","shell.execute_reply.started":"2024-08-28T10:17:59.168356Z","shell.execute_reply":"2024-08-28T10:18:01.06485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We divide the images by the star flux assuming the first and last 50 instants belong to the out of transit. ","metadata":{}},{"cell_type":"code","source":"def create_dataset_norm(dataset1, dataset2) :\n    dataset_norm1 = np.zeros(dataset1.shape)\n    dataset_norm2 = np.zeros(dataset1.shape)\n    dataset_min = dataset1.min()\n    dataset_max = dataset1.max()\n    dataset_norm1 = (dataset1 - dataset_min) / (dataset_max - dataset_min)\n    dataset_norm2 = (dataset2 - dataset_min) / (dataset_max - dataset_min)\n    return dataset_norm1, dataset_norm2\n\n\ndef norm_star_spectrum (signal) : \n    img_star = signal[:,:50].mean(axis = 1) + signal[:,-50:].mean(axis = 1)\n    return signal/img_star[:,np.newaxis,:]\n\ndataset_norm = norm_star_spectrum(dataset)\ndataset_norm = np.transpose(dataset_norm,(0,2,1))\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:01.066858Z","iopub.execute_input":"2024-08-28T10:18:01.067134Z","iopub.status.idle":"2024-08-28T10:18:01.188333Z","shell.execute_reply.started":"2024-08-28T10:18:01.067111Z","shell.execute_reply":"2024-08-28T10:18:01.187533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Split the targets and observations between valid and train\n\nWe start by computing a \"white curve\", that is actually the sum of the signal over the all image, as a function of time. We split the data and normalize the train/valid/test data.","metadata":{}},{"cell_type":"code","source":"cut_inf, cut_sup = 39, 321 # we have previously cut the data along the wavelengths to remove the edges, this is to match with the targets range in the make data file\nl = cut_sup - cut_inf + 1 \nwls = np.arange(l)\n\n\ndef split (data, N) : \n    list_planets = random.sample(range(0, data.shape[0]), N_train)\n    list_index_1 = np.zeros(data.shape[0], dtype = bool)\n    for planet in list_planets : \n        list_index_1[planet] = True\n    data_1 = data[list_index_1]\n    data_2 = data[~list_index_1]\n    return data_1, data_2, list_index_1\n\nN_train = 8*N//10\n\n# Validation and train data split\ntrain_obs, valid_obs, list_index_train = split(dataset_norm, N_train)\ntrain_targets, valid_targets = targets[list_index_train], targets[~list_index_train]\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:01.192252Z","iopub.execute_input":"2024-08-28T10:18:01.192606Z","iopub.status.idle":"2024-08-28T10:18:01.282142Z","shell.execute_reply.started":"2024-08-28T10:18:01.192577Z","shell.execute_reply":"2024-08-28T10:18:01.281115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal_AIRS_diff_transposed_binned = signal_AIRS_diff_transposed_binned.sum(axis=3)\nwc_mean = signal_AIRS_diff_transposed_binned.mean(axis=1).mean(axis=1)\nwhite_curve = signal_AIRS_diff_transposed_binned.sum(axis=2)/ wc_mean[:, np.newaxis]\n\ndef normalise_wlc(train, valid) :\n\n    wlc_train_min = train.min()\n    wlc_train_max = train.max()\n    train_norm = (train - wlc_train_min) / (wlc_train_max - wlc_train_min)\n    valid_norm = (valid - wlc_train_min) / (wlc_train_max - wlc_train_min)\n    \n    return train_norm, valid_norm\n\ndef normalize (train, valid) : \n    max_train = train.max()\n    min_train = train.min()\n    train_norm = (train - min_train) / (max_train - min_train)\n    valid_norm = (valid - min_train) / (max_train - min_train)\n    return train_norm, valid_norm, min_train, max_train\n\n# Split the light curves and targets \ntrain_wc, valid_wc = white_curve[list_index_train], white_curve[~list_index_train]\ntrain_targets_wc, valid_targets_wc = targets_mean[list_index_train], targets_mean[~list_index_train]\n\n# Normalize the wlc\ntrain_wc, valid_wc = normalise_wlc(train_wc, valid_wc)\n\n# Normalize the targets \ntrain_targets_wc_norm, valid_targets_wc_norm, min_train_valid_wc, max_train_valid_wc = normalize(train_targets_wc, valid_targets_wc)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:01.283252Z","iopub.execute_input":"2024-08-28T10:18:01.283578Z","iopub.status.idle":"2024-08-28T10:18:02.718479Z","shell.execute_reply.started":"2024-08-28T10:18:01.283552Z","shell.execute_reply":"2024-08-28T10:18:02.717674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nfor i in range (200) : \n    plt.plot(train_wc[-i], '-', alpha = 0.5)\nplt.title('Light-curves from the train set') \nplt.xlabel('Time')\nplt.ylabel('Normalized flux')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:02.719569Z","iopub.execute_input":"2024-08-28T10:18:02.719855Z","iopub.status.idle":"2024-08-28T10:18:03.188177Z","shell.execute_reply.started":"2024-08-28T10:18:02.71983Z","shell.execute_reply":"2024-08-28T10:18:03.187233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train 1D CNN\nThe model to estimate the mean of the target spectrum using the white light-curve is a 1D-CNN with Dropout layers to make a MC Dropout prediction.","metadata":{}},{"cell_type":"code","source":"from keras.layers import Input, Conv1D, MaxPooling1D, Flatten, Dense, Dropout, BatchNormalization, Concatenate,AveragePooling1D\nfrom keras.models import Model, load_model\nfrom tensorflow.keras.optimizers import Adam, SGD\nfrom tensorflow.keras.callbacks import LearningRateScheduler, ModelCheckpoint\n\n\ninput_wc = Input((187,1))\nx = Conv1D(32, 3, activation='relu')(input_wc)\nx = MaxPooling1D()(x)\nx = BatchNormalization() (x)\nx = Conv1D(64, 3, activation='relu')(x)\nx = MaxPooling1D()(x)\nx = Conv1D(128, 3, activation='relu')(x)\nx = MaxPooling1D()(x)\nx = Conv1D(256, 3, activation='relu')(x)\nx = MaxPooling1D()(x)\nx = Flatten()(x)\n\nx = Dense(500, activation='relu')(x)\nx = Dropout(0.2)(x, training = True)\nx = Dense(100, activation='relu')(x)\nx = Dropout(0.1)(x, training = True)\noutput_wc = Dense(1, activation='linear')(x)\n\nmodel_wc = Model(inputs=input_wc, outputs=output_wc)\nmodel_wc.summary()\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:03.189385Z","iopub.execute_input":"2024-08-28T10:18:03.189684Z","iopub.status.idle":"2024-08-28T10:18:04.030675Z","shell.execute_reply.started":"2024-08-28T10:18:03.189658Z","shell.execute_reply":"2024-08-28T10:18:04.029697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def scheduler(epoch, lr):\n    decay_rate = 0.2\n    decay_step = 200  \n    if epoch % decay_step == 0 and epoch:\n        return lr * decay_rate\n    return lr\n\noptimizer = SGD(0.001)\nmodel_wc.compile(optimizer=optimizer, loss='mse', metrics=[MeanAbsoluteError()])\ncallback = LearningRateScheduler(scheduler)\ncheckpoint_filepath = 'output/model_1dcnn.keras'\nmodel_ckt = ModelCheckpoint(\n    checkpoint_filepath,\n    monitor=\"val_loss\",\n    verbose=0,\n    save_best_only=True,\n    save_weights_only=False,\n    mode=\"min\",\n    save_freq=\"epoch\",\n)\n\nprint('Running ...')\nhistory = model_wc.fit(\n    x = train_wc,  \n    y = train_targets_wc_norm,\n    validation_data = (valid_wc, valid_targets_wc_norm),  \n    batch_size=16,\n    epochs= 1200,\n    shuffle=True,\n    verbose=0, \n    callbacks=[model_ckt]\n    )\nprint('Done.')\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:18:04.032106Z","iopub.execute_input":"2024-08-28T10:18:04.032945Z","iopub.status.idle":"2024-08-28T10:20:58.861904Z","shell.execute_reply.started":"2024-08-28T10:18:04.032905Z","shell.execute_reply":"2024-08-28T10:20:58.86087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1D CNN Inference","metadata":{}},{"cell_type":"code","source":"model_wc = load_model(checkpoint_filepath)\n\nplt.plot(history.history['loss'], label='Training Loss')\nplt.plot(history.history['val_loss'], label='Validation Loss')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:20:58.863523Z","iopub.execute_input":"2024-08-28T10:20:58.864384Z","iopub.status.idle":"2024-08-28T10:20:59.243754Z","shell.execute_reply.started":"2024-08-28T10:20:58.864344Z","shell.execute_reply":"2024-08-28T10:20:59.242817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then, we perform the MC Dropout to obtain the mean prediction and the uncertainty associated. We choose to compute 1000 instances.","metadata":{}},{"cell_type":"code","source":"nb_dropout_wc = 1000\n\ndef unstandardizing (data, min_train_valid, max_train_valid) : \n    return data * (max_train_valid - min_train_valid) + min_train_valid\n\ndef MC_dropout_WC (model, data, nb_dropout) : \n    predictions = np.zeros((nb_dropout, data.shape[0]))\n    for i in range(nb_dropout) : \n        predictions[i,:] = model.predict(data, verbose = 0).flatten()\n    return predictions\n\nif do_the_mcdropout_wc :\n    print('Running ...')\n    prediction_valid_wc = MC_dropout_WC(model_wc, valid_wc, nb_dropout_wc)\n    spectre_valid_wc_all = unstandardizing(prediction_valid_wc, min_train_valid_wc, max_train_valid_wc)\n    spectre_valid_wc, spectre_valid_std_wc = spectre_valid_wc_all.mean(axis = 0), spectre_valid_wc_all.std(axis = 0)\n    print('Done.')\n\nelse : \n    spectre_valid_wc = model_wc.predict(valid_wc).flatten()\n    spectre_valid_wc = unstandardizing(spectre_valid_wc, min_train_valid_wc, max_train_valid_wc)\n    spectre_valid_std_wc = 0.1*np.abs(spectre_valid_wc)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:20:59.245086Z","iopub.execute_input":"2024-08-28T10:20:59.245395Z","iopub.status.idle":"2024-08-28T10:22:06.348718Z","shell.execute_reply.started":"2024-08-28T10:20:59.245368Z","shell.execute_reply":"2024-08-28T10:22:06.347692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"residuals = spectre_valid_wc - valid_targets_wc\nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True,\n                               gridspec_kw={'height_ratios': [3, 1]})\n\nax1.errorbar(x = np.arange(len(spectre_valid_wc)), y = spectre_valid_wc, yerr =spectre_valid_std_wc, fmt = '.', color = 'k', ecolor = 'gray', label='Prediction', alpha=0.8)\nax1.fill_between(np.arange(len(spectre_valid_wc)), spectre_valid_wc - spectre_valid_std_wc, spectre_valid_wc + spectre_valid_std_wc, color = 'grey', alpha = 0.5)\nax1.vlines(np.arange(len(spectre_valid_wc)),ymin=0, ymax=spectre_valid_wc, colors='r', linestyle='dashed',alpha = 0.1)\nax1.plot(valid_targets_wc, 'r.', label='Target', alpha=0.8)\nax1.set_xlabel('Concatenated targets')\nax1.set_ylabel('$(R_p/R_s)^2$')\nax1.set_title('Prediction vs target, mean value of the spectrum, on validation dataset')\nax1.legend()\n\nax2.plot(residuals, 'b.', label='Residuals', alpha=0.8)\nax2.set_xlabel('Concatenated targets')\nax2.set_ylabel('Residuals')\nax2.axhline(0, color='black', linestyle='--', linewidth=1)  \nax2.legend()\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:06.349899Z","iopub.execute_input":"2024-08-28T10:22:06.350187Z","iopub.status.idle":"2024-08-28T10:22:07.161288Z","shell.execute_reply.started":"2024-08-28T10:22:06.350161Z","shell.execute_reply":"2024-08-28T10:22:07.160204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"residuals = valid_targets_wc - spectre_valid_wc\nprint('MSE : ', np.sqrt((residuals**2).mean())*1e6, 'ppm')","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.162458Z","iopub.execute_input":"2024-08-28T10:22:07.16275Z","iopub.status.idle":"2024-08-28T10:22:07.168677Z","shell.execute_reply.started":"2024-08-28T10:22:07.162727Z","shell.execute_reply":"2024-08-28T10:22:07.167535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# np.save(f'{output_dir}/pred_valid_wc.npy', spectre_valid_wc)\n# np.save(f'{output_dir}/targ_valid_wc.npy', valid_targets_wc)\n# np.save(f'{output_dir}/std_valid_wc.npy', spectre_valid_std_wc)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.169736Z","iopub.execute_input":"2024-08-28T10:22:07.17009Z","iopub.status.idle":"2024-08-28T10:22:07.180378Z","shell.execute_reply.started":"2024-08-28T10:22:07.170064Z","shell.execute_reply":"2024-08-28T10:22:07.17937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2D CNN for atmospheric features\n\n<a id=\"fluctu\"></a>\nWe now remove the mean value (transit depth) of the spectra to keep the atmospheric features only","metadata":{}},{"cell_type":"markdown","source":"## Preprocessing for 2D CNN","metadata":{}},{"cell_type":"code","source":"def suppress_mean(targets, mean) : \n    res = targets - np.repeat(mean.reshape((mean.shape[0], 1)), repeats = targets.shape[1], axis = 1)\n    return res\ntrain_targets, valid_targets = targets[list_index_train], targets[~list_index_train]\n\ntrain_targets_shift = suppress_mean(train_targets,  targets_mean[list_index_train])\nvalid_targets_shift = suppress_mean(valid_targets,  targets_mean[~list_index_train])","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.181654Z","iopub.execute_input":"2024-08-28T10:22:07.182216Z","iopub.status.idle":"2024-08-28T10:22:07.192868Z","shell.execute_reply.started":"2024-08-28T10:22:07.182188Z","shell.execute_reply":"2024-08-28T10:22:07.192006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We normalize the targets so that they range between -1 and 1, centered on zero","metadata":{}},{"cell_type":"code","source":"##### normalization of the targets ###\ndef targets_normalization (data1, data2) : \n    data_min = data1.min()\n    data_max = data1.max()\n    data_abs_max = np.max([data_min, data_max])  \n    data1 = data1/data_abs_max\n    data2 = data2/data_abs_max\n    return data1, data2, data_abs_max\n\ndef targets_norm_back (data, data_abs_max) : \n    return data * data_abs_max\n\ntrain_targets_norm, valid_targets_norm, targets_abs_max = targets_normalization(train_targets_shift, valid_targets_shift)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.194238Z","iopub.execute_input":"2024-08-28T10:22:07.195007Z","iopub.status.idle":"2024-08-28T10:22:07.20338Z","shell.execute_reply.started":"2024-08-28T10:22:07.194979Z","shell.execute_reply":"2024-08-28T10:22:07.202373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,5))\n\nfor i in range (240) :\n    plt.plot(wls, train_targets_norm[i], 'g-', alpha = 0.5)\nplt.plot([], [], 'g-', alpha=0.5, label='Train targets')\nfor i in range (60) : \n    plt.plot(wls, valid_targets_norm[i], 'r-', alpha = 0.7)\nplt.plot([], [], 'r-', alpha=0.5, label='Valid targets (true mean)')\n\nplt.legend()\nplt.ylabel(f'$(R_p/R_s)^2$')\nplt.title('All targets after substracting the mean value and normalization')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.204723Z","iopub.execute_input":"2024-08-28T10:22:07.205369Z","iopub.status.idle":"2024-08-28T10:22:07.990582Z","shell.execute_reply.started":"2024-08-28T10:22:07.205336Z","shell.execute_reply":"2024-08-28T10:22:07.989574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"###### Transpose #####\ntrain_obs = train_obs.transpose(0, 2, 1)\nvalid_obs = valid_obs.transpose(0, 2, 1)\nprint(train_obs.shape)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.992149Z","iopub.execute_input":"2024-08-28T10:22:07.992971Z","iopub.status.idle":"2024-08-28T10:22:07.998457Z","shell.execute_reply.started":"2024-08-28T10:22:07.992933Z","shell.execute_reply":"2024-08-28T10:22:07.99745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We cut the transit to keep the in-transit. We assume an arbitrary transit duration of 40 instants with a transit occuring between 75 and 115.","metadata":{}},{"cell_type":"code","source":"##### Substracting the out transit signal #####\ndef suppress_out_transit (data, ingress, egress) : \n    data_in = data[:, ingress:egress,:]\n    return data_in\n\ningress, egress = 75,115\ntrain_obs_in = suppress_out_transit(train_obs, ingress, egress)\nvalid_obs_in = suppress_out_transit(valid_obs, ingress, egress)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:07.999784Z","iopub.execute_input":"2024-08-28T10:22:08.000657Z","iopub.status.idle":"2024-08-28T10:22:08.008246Z","shell.execute_reply.started":"2024-08-28T10:22:08.000628Z","shell.execute_reply":"2024-08-28T10:22:08.007517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We remove the mean value of the in-transit to get relative data like the targets","metadata":{}},{"cell_type":"code","source":"###### Substract the mean #####\ndef substract_data_mean(data):\n    data_mean = np.zeros(data.shape)\n    for i in range(data.shape[0]):\n        data_mean[i] = data[i] - data[i].mean()\n    return data_mean\n\ntrain_obs_2d_mean = substract_data_mean(train_obs_in)\nvalid_obs_2d_mean = substract_data_mean(valid_obs_in)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:08.009383Z","iopub.execute_input":"2024-08-28T10:22:08.00976Z","iopub.status.idle":"2024-08-28T10:22:08.080069Z","shell.execute_reply.started":"2024-08-28T10:22:08.009736Z","shell.execute_reply":"2024-08-28T10:22:08.07925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use the same normalization as for the targets, i.e. between -1 and 1 centered on zero","metadata":{}},{"cell_type":"code","source":"##### Normalization dataset #####\ndef data_norm(data1, data2):\n    data_min = data1.min()\n    data_max = data1.max()\n    data_abs_max = np.max([data_min, data_max])  \n    data1 = data1/data_abs_max\n    data2 = data2/data_abs_max\n    return data1, data2, data_abs_max\n\n\ndef data_normback(data, data_abs_max) : \n    return data * data_abs_max\n\ntrain_obs_norm, valid_obs_norm, data_abs_max = data_norm(train_obs_2d_mean, valid_obs_2d_mean)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:08.08118Z","iopub.execute_input":"2024-08-28T10:22:08.081484Z","iopub.status.idle":"2024-08-28T10:22:08.114991Z","shell.execute_reply.started":"2024-08-28T10:22:08.081458Z","shell.execute_reply":"2024-08-28T10:22:08.11394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,5))\nfor i in range (train_obs.shape[0]) :\n    plt.plot(wls, train_obs_norm[i,10], 'g-', alpha = 0.5)\nplt.plot([], [], 'g-', alpha=0.5, label='Train targets')\nfor i in range (valid_obs.shape[0]) : \n    plt.plot(wls, valid_obs_norm[i,10], 'r-', alpha = 0.7)\nplt.plot([], [], 'r-', alpha=0.5, label='Valid targets (true mean)')\n\nplt.legend()\nplt.ylabel(f'$(R_p/R_s)^2$')\nplt.title('Train and Valid data after substracting the mean value and normalization')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:08.116263Z","iopub.execute_input":"2024-08-28T10:22:08.116589Z","iopub.status.idle":"2024-08-28T10:22:10.439434Z","shell.execute_reply.started":"2024-08-28T10:22:08.116563Z","shell.execute_reply":"2024-08-28T10:22:10.438411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train 2D CNN","metadata":{}},{"cell_type":"code","source":"from tensorflow import keras\nfrom keras.layers import Input, Conv2D, MaxPooling2D, Flatten, Dense, Concatenate, Reshape, Dropout, BatchNormalization, AveragePooling2D\nfrom keras.models import Model\nimport tensorflow as tf\nimport numpy as np\n\n## CNN 2 global normalization data\ninput_obs = Input((40,283,1))\nx = Conv2D(32, (3, 1), activation='relu', padding='same')(input_obs)\nx = MaxPooling2D((2, 1))(x)\nx = BatchNormalization() (x)\nx = Conv2D(64, (3, 1), activation='relu', padding='same')(x)\nx = MaxPooling2D((2, 1))(x)\nx = Conv2D(128, (3, 1), activation='relu', padding='same')(x)\nx = MaxPooling2D((2, 1))(x)\nx = Conv2D(256, (3, 1), activation='relu', padding='same')(x)\nx = Conv2D(32, (1, 3), activation='relu', padding='same')(x)\nx = MaxPooling2D((1, 2))(x)\nx = BatchNormalization() (x)\nx = Conv2D(64, (1, 3), activation='relu', padding='same')(x)\nx = MaxPooling2D((1, 2))(x)\nx = Conv2D(128, (1, 3), activation='relu', padding='same')(x)\nx = MaxPooling2D((1, 2))(x)\nx = Conv2D(256, (1, 3), activation='relu', padding='same')(x)\nx = MaxPooling2D((1, 2))(x)\nx = Flatten()(x)\n# DNN\nx = Dense(700, activation='relu')(x)\nx = Dropout(0.2)(x, training = True)\noutput = Dense(283, activation='linear')(x)\n\nmodel = Model(inputs=[input_obs], outputs=output)\n\ncheckpoint_filepath = 'output/model_2dcnn.keras'\nmodel_ckt2 = ModelCheckpoint(\n    checkpoint_filepath,\n    monitor=\"val_loss\",\n    verbose=0,\n    save_best_only=True,\n    save_weights_only=False,\n    mode=\"min\",\n    save_freq=\"epoch\",\n)\nmodel.compile(optimizer=Adam(0.001), loss='mse', metrics=[MeanAbsoluteError()])\nmodel.summary()","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:10.444319Z","iopub.execute_input":"2024-08-28T10:22:10.444641Z","iopub.status.idle":"2024-08-28T10:22:10.602608Z","shell.execute_reply.started":"2024-08-28T10:22:10.444614Z","shell.execute_reply":"2024-08-28T10:22:10.601631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = model.fit(\n     x = train_obs_norm,  \n     y = train_targets_norm,\n     validation_data = (valid_obs_norm, valid_targets_norm),  \n     batch_size=32,\n     epochs= 200,\n     shuffle=True,\n     verbose=0,\n    callbacks=[model_ckt2]\n )\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:22:10.603843Z","iopub.execute_input":"2024-08-28T10:22:10.604205Z","iopub.status.idle":"2024-08-28T10:24:56.643548Z","shell.execute_reply.started":"2024-08-28T10:22:10.604167Z","shell.execute_reply":"2024-08-28T10:24:56.6424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Postprocessing and visualisation","metadata":{}},{"cell_type":"markdown","source":"We obtain uncertainties on the predictions by computing a MCDropout.","metadata":{}},{"cell_type":"code","source":"nb_dropout = 5\n\ndef NN_uncertainity(model, x_test, targets_abs_max, T=5):\n    predictions = []\n    for _ in range(T):\n        pred_norm = model.predict([x_test],verbose=0)\n        pred = targets_norm_back(pred_norm, targets_abs_max)\n        predictions += [pred]  \n    mean, std = np.mean(np.array(predictions), axis=0), np.std(np.array(predictions), axis=0)\n    return mean, std\n\n\nif do_the_mcdropout :\n    spectre_valid_shift, spectre_valid_shift_std = NN_uncertainity(model, [valid_obs_norm], targets_abs_max, T = nb_dropout)\n    \nelse :\n\n    pred_valid_norm = model.predict([valid_obs_norm])\n    pred_valid = targets_norm_back(pred_valid_norm, targets_abs_max)\n    spectre_valid_shift = pred_valid\n    spectre_valid_shift_std = spectre_valid_shift*0.1\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:56.64575Z","iopub.execute_input":"2024-08-28T10:24:56.646064Z","iopub.status.idle":"2024-08-28T10:24:58.094042Z","shell.execute_reply.started":"2024-08-28T10:24:56.646036Z","shell.execute_reply":"2024-08-28T10:24:58.092907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"residuals = valid_targets_shift - spectre_valid_shift\nprint('MSE : ', np.sqrt((residuals**2).mean())*1e6, 'ppm')","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:58.095232Z","iopub.execute_input":"2024-08-28T10:24:58.095552Z","iopub.status.idle":"2024-08-28T10:24:58.101838Z","shell.execute_reply.started":"2024-08-28T10:24:58.095524Z","shell.execute_reply":"2024-08-28T10:24:58.100803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# np.save(f'{output_dir}/pred_valid_shift.npy', spectre_valid_shift)\n# np.save(f'{output_dir}/targ_valid_shift.npy', valid_targets_shift)\n# np.save(f'{output_dir}/std_valid_shift.npy', spectre_valid_shift_std)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:58.103053Z","iopub.execute_input":"2024-08-28T10:24:58.103358Z","iopub.status.idle":"2024-08-28T10:24:58.112147Z","shell.execute_reply.started":"2024-08-28T10:24:58.103327Z","shell.execute_reply":"2024-08-28T10:24:58.111179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nfor i in range (50) : \n    plt.plot(spectre_valid_shift[-i]+0.0001*i, '-', alpha = 0.5)\nplt.title('Spectra predictions for the validation set') \nplt.xlabel('Time')\nplt.ylabel('Arbitrary flux')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:58.113246Z","iopub.execute_input":"2024-08-28T10:24:58.113854Z","iopub.status.idle":"2024-08-28T10:24:58.460699Z","shell.execute_reply.started":"2024-08-28T10:24:58.113828Z","shell.execute_reply":"2024-08-28T10:24:58.459616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_valid_planets = [0, 12, 35, 60, 70]\nwavelength = np.loadtxt('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv', skiprows=1, delimiter = ',')\nuncertainty = spectre_valid_shift_std\nfor i in (list_valid_planets):\n    plt.figure()\n    plt.title('Result for the sample {} of the validation set'.format(i))\n    plt.plot(wavelength, spectre_valid_shift[i], '.k', label = 'Prediction')\n    plt.plot(wavelength, valid_targets_shift[i], color = 'tomato', label = 'Target')\n    plt.fill_between(wavelength, spectre_valid_shift[i] - spectre_valid_shift_std[i], spectre_valid_shift[i] + spectre_valid_shift_std[i], color='silver', alpha = 0.8, label = 'Uncertainty')\n    plt.legend()\n    plt.ylabel(f'$(R_p/R_s)^2$')\n    plt.xlabel(f'Wavelength ($\\mu$m)')\n    plt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:58.462279Z","iopub.execute_input":"2024-08-28T10:24:58.462615Z","iopub.status.idle":"2024-08-28T10:24:59.960337Z","shell.execute_reply.started":"2024-08-28T10:24:58.462586Z","shell.execute_reply":"2024-08-28T10:24:59.959345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Combine 1D and 2D CNN output for FINAL SPECTRA ","metadata":{}},{"cell_type":"code","source":"######## ADD THE FLUCTUATIONS TO THE MEAN ########\ndef add_the_mean (shift, mean) : \n    return shift + mean[:,np.newaxis]\n\npredictions_valid = add_the_mean(spectre_valid_shift,spectre_valid_wc)\n\npredictions_std_valid = np.sqrt(spectre_valid_std_wc[:,np.newaxis]**2 + spectre_valid_shift_std**2)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:59.961581Z","iopub.execute_input":"2024-08-28T10:24:59.961854Z","iopub.status.idle":"2024-08-28T10:24:59.96746Z","shell.execute_reply.started":"2024-08-28T10:24:59.96183Z","shell.execute_reply":"2024-08-28T10:24:59.966488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"uncertainty = predictions_std_valid\n\ndef plot_one_sample_valid(ax, p):\n    ax.set_title(f'Result for sample {p} ')\n    line1, = ax.plot(wavelength, predictions_valid[p], '.k', label='Prediction')\n    line2, = ax.plot(wavelength, valid_targets[p], color='tomato', label='Target')\n    ax.fill_between(wavelength, predictions_valid[p, :] - uncertainty[p], predictions_valid[p, :] + uncertainty[p], color='silver', alpha=0.8, label='Uncertainty')\n    ax.set_ylabel(f'$(R_p/R_s)^2$')\n    ax.set_xlabel(f'Wavelength ($\\mu$m)')\n    return line1, line2\n\n\nnum_samples = 16\nrows, cols = 4, 4\n\nfig, axs = plt.subplots(rows, cols, figsize=(15, 10))\nsamples = [1, 2, 7, 15, 20, 25, 30, 35, 40, 45, 50, 55, 6, 5, 8, 9]\nlines = []\n\nfor i, ax in enumerate(axs.flat):\n    lines.extend(plot_one_sample_valid(ax, samples[i]))\n\nfig.legend(lines[:2], ['Prediction', 'Target'], loc='upper center', ncol=3, bbox_to_anchor=(0.5, -0.05))\nfig.suptitle('Validation dataset')\nplt.tight_layout()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:24:59.96877Z","iopub.execute_input":"2024-08-28T10:24:59.969454Z","iopub.status.idle":"2024-08-28T10:25:03.140594Z","shell.execute_reply.started":"2024-08-28T10:24:59.969419Z","shell.execute_reply":"2024-08-28T10:25:03.139626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"######## PLOTS THE RESULT ########\npredictions = predictions_valid\ntargets_plot = valid_targets\nstd = predictions_std_valid\n\npredictions_concatenated_plot = np.concatenate(predictions, axis=0)\nwls_concatenated = np.arange(predictions_concatenated_plot.shape[0])\ntargets_concatenated_plot = np.concatenate(targets_plot, axis=0)\nspectre_valid_std_concatenated = np.concatenate(std, axis=0)\nresiduals = targets_concatenated_plot - predictions_concatenated_plot\nuncertainty = spectre_valid_std_concatenated\n\nfig, axs = plt.subplots(2, 1, figsize=(9, 8), gridspec_kw={'height_ratios': [3, 1]})\n\n\naxs[0].plot(wls_concatenated, predictions_concatenated_plot, '-', color='k', label=\"Prediction\")\naxs[0].plot(wls_concatenated, targets_concatenated_plot, '-', color='tomato', label=\"Target\")\naxs[0].fill_between(np.arange(len(wls_concatenated)), \n                    predictions_concatenated_plot - uncertainty, \n                    predictions_concatenated_plot + uncertainty, \n                    color='silver', alpha=1, label='Uncertainty')\naxs[0].set_xlabel('Concatenated wavelengths for all planets')\naxs[0].set_ylabel(f'$(R_p/R_s)^2$')\naxs[0].set_title('Prediction vs target, validation dataset')\naxs[0].legend()\n\naxs[1].plot(wls_concatenated, residuals, '-', color='cornflowerblue', label=\"Residual\")\naxs[1].fill_between(np.arange(len(wls_concatenated)), \n                    residuals - uncertainty, \n                    residuals + uncertainty, \n                    color='lightblue', alpha=0.9, label='Uncertainty')\naxs[1].set_xlabel('Concatenated wavelengths for all planets')\naxs[1].set_ylabel('Residual')\naxs[1].set_title('Residuals with Uncertainty')\naxs[1].legend()\n\nplt.tight_layout()\nplt.show()\n\nprint('MSE : ',np.sqrt((residuals**2).mean())*1e6, 'ppm')","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:25:03.14191Z","iopub.execute_input":"2024-08-28T10:25:03.142266Z","iopub.status.idle":"2024-08-28T10:25:04.306814Z","shell.execute_reply.started":"2024-08-28T10:25:03.142235Z","shell.execute_reply":"2024-08-28T10:25:04.305898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# np.save(f'{output_dir}/pred_valid.npy', predictions_valid)\n# np.save(f'{output_dir}/std_valid.npy', predictions_std_valid)","metadata":{"execution":{"iopub.status.busy":"2024-08-28T10:25:04.307962Z","iopub.execute_input":"2024-08-28T10:25:04.308233Z","iopub.status.idle":"2024-08-28T10:25:04.312201Z","shell.execute_reply.started":"2024-08-28T10:25:04.308209Z","shell.execute_reply":"2024-08-28T10:25:04.311269Z"},"trusted":true},"execution_count":null,"outputs":[]}]}