{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.13"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":204906070,"sourceType":"kernelVersion"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":1087.424982,"end_time":"2024-08-03T13:00:18.164826","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-08-03T12:42:10.739844","version":"2.5.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\nimport pickle\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\n\nimport os\nimport itertools\nfrom scipy.interpolate import splrep, BSpline\nimport multiprocessing\nimport plotly.express as px\nfrom matplotlib  import cm\n\nfrom scipy.optimize import curve_fit \nfrom mpl_toolkits.mplot3d import Axes3D ","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-11-03T14:45:09.685576Z","iopub.execute_input":"2024-11-03T14:45:09.686088Z","iopub.status.idle":"2024-11-03T14:45:13.841303Z","shell.execute_reply.started":"2024-11-03T14:45:09.686041Z","shell.execute_reply":"2024-11-03T14:45:13.840378Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv',\n                           index_col='planet_id')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',\n                           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')\nstar = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv').star.to_numpy()","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-11-03T14:45:13.843335Z","iopub.execute_input":"2024-11-03T14:45:13.843867Z","iopub.status.idle":"2024-11-03T14:45:14.178307Z","shell.execute_reply.started":"2024-11-03T14:45:13.843838Z","shell.execute_reply":"2024-11-03T14:45:14.177382Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"advanced_fitting_results = pickle.load(open('/kaggle/input/adc24-preprocessing/advanced_fitting_results.p', 'br'))\nresults_fit_on_sections = pickle.load(open('/kaggle/input/adc24-preprocessing/results_fit_on_sections.p', 'br'))\ndata_list_pols = pickle.load(open('/kaggle/input/adc24-preprocessing/data_list_pols.p', 'br'))\ndata_list_pols = np.asarray(data_list_pols)","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:14.179752Z","iopub.execute_input":"2024-11-03T14:45:14.18007Z","iopub.status.idle":"2024-11-03T14:45:27.740092Z","shell.execute_reply.started":"2024-11-03T14:45:14.180042Z","shell.execute_reply":"2024-11-03T14:45:27.739005Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"baseline_list = [x[0] for x in advanced_fitting_results]\ndelta_list = [x[1] for x in advanced_fitting_results]\nerr_list = [x[2] for x in advanced_fitting_results]\nbaseline_list_2 = [x[3] for x in advanced_fitting_results]\ndelta_list_2 = [x[4] for x in advanced_fitting_results]\nerr_list_2 = [x[5] for x in advanced_fitting_results]\n\nrelred1 = np.asarray(delta_list)/np.asarray(baseline_list)\nrelred2 = np.asarray(delta_list_2)/np.asarray(baseline_list_2)\nrelred1 = relred1[:, ::-1]\nrelred2 = relred2[:, ::-1]","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.743016Z","iopub.execute_input":"2024-11-03T14:45:27.743654Z","iopub.status.idle":"2024-11-03T14:45:27.819826Z","shell.execute_reply.started":"2024-11-03T14:45:27.743618Z","shell.execute_reply":"2024-11-03T14:45:27.818992Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"targets = train_labels.to_numpy()\ntargets_orig = targets.copy()","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.820891Z","iopub.execute_input":"2024-11-03T14:45:27.821249Z","iopub.status.idle":"2024-11-03T14:45:27.826185Z","shell.execute_reply.started":"2024-11-03T14:45:27.821223Z","shell.execute_reply":"2024-11-03T14:45:27.825164Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# SCORES/validation","metadata":{}},{"cell_type":"code","source":"%%writefile postprocessing.py\n\ndef postprocessing(pred_array, index, sigma_pred):\n    \"\"\"Create a submission dataframe from its components\n    \n    Parameters:\n    pred_array: ndarray of shape (n_samples, 283)\n    index: pandas.Index of length n_samples with name 'planet_id'\n    sigma_pred: float\n    \n    Return value:\n    df: DataFrame of shape (n_samples, 566) with planet_id as index\n    \"\"\"\n    return pd.concat([pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns),\n                      pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)])],\n                     axis=1)","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.82739Z","iopub.execute_input":"2024-11-03T14:45:27.827651Z","iopub.status.idle":"2024-11-03T14:45:27.84548Z","shell.execute_reply.started":"2024-11-03T14:45:27.827627Z","shell.execute_reply":"2024-11-03T14:45:27.844434Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile competition_score.py\n# Adapted from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef competition_score(\n        solution: pd.DataFrame,\n        submission: pd.DataFrame,\n        naive_mean: float,\n        naive_sigma: float,\n        sigma_true: float,\n        row_id_column_name='planet_id',\n    ) -> float:\n    '''\n    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        sigma_true: (float) essentially sets the scale of the outputs.\n    '''\n\n    del solution[row_id_column_name]\n    del submission[row_id_column_name]\n\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n    for col in submission.columns:\n        if not pd.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be a number')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != n_wavelengths*2:\n        raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n    y_true = solution.values\n\n    GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true)))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true)))\n\n    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    \n    GLL_pred_indv = scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    return float(np.clip(submit_score, 0.0, 1.0)), (GLL_pred_indv - GLL_mean/(len(GLL_pred_indv)*len(GLL_pred_indv[0])))/(GLL_true - GLL_mean)","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.846809Z","iopub.execute_input":"2024-11-03T14:45:27.847131Z","iopub.status.idle":"2024-11-03T14:45:27.859108Z","shell.execute_reply.started":"2024-11-03T14:45:27.847103Z","shell.execute_reply":"2024-11-03T14:45:27.857984Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('postprocessing.py', 'r').read())\nexec(open('competition_score.py', 'r').read())\nadc_info =  train_adc_info\nnaive_mean = train_labels.values.mean()\nnaive_sigma = train_labels.values.std()","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.860443Z","iopub.execute_input":"2024-11-03T14:45:27.860747Z","iopub.status.idle":"2024-11-03T14:45:27.875236Z","shell.execute_reply.started":"2024-11-03T14:45:27.86072Z","shell.execute_reply":"2024-11-03T14:45:27.874371Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile model_1_postprocessing.py\ndef model_1_postprocessing(trelred1, relred2, results_fit_on_sections, targets = None, test = True, data_path = None):\n    def get_freq_averaged(relred1, relred2, freq_window):\n        laa = relred1.copy()\n        signaL_mean_freq = []\n        for freq_idx in range(282):\n            min_freq = max(0, freq_idx-freq_window)\n            max_freq = min(freq_idx+freq_window+1, 282)\n            signaL_mean_freq.append(np.mean(laa[:, min_freq:max_freq], axis = 1))\n        train_1 = np.transpose(np.asarray(signaL_mean_freq))\n        laa = relred2.copy()\n        signaL_mean_freq = []\n        for freq_idx in range(282):\n            min_freq = max(0, freq_idx-freq_window)\n            max_freq = min(freq_idx+freq_window+1, 282)\n            signaL_mean_freq.append(np.mean(laa[:, min_freq:max_freq], axis = 1))\n        train_2 = np.transpose(np.asarray(signaL_mean_freq))\n        train_1 = np.concatenate([train_1[:, :1], train_1], axis = 1)\n        train_2 = np.concatenate([train_2[:, :1], train_2], axis = 1)\n        train = (train_1+train_2)/2\n        return train_1, train_2, train\n\n    def get_results_fit_on_sections_data(results_fit_on_sections):\n        data_list = []\n        for idx in range(len(results_fit_on_sections)):\n            [rms_list_1, fit_pol_list_1, rms_list_2, fit_pol_list_2,\n                        rms_list_3, fit_pol_list_3, rms_list_4, fit_pol_list_4,\n                        rms_list_5, fit_pol_list_5, rms_list_6, fit_pol_list_6] = results_fit_on_sections[idx]\n            fit_pol_list_2 = fit_pol_list_2[::-1]\n            fit_pol_list_4 = fit_pol_list_4[::-1]\n            fit_pol_list_6 = fit_pol_list_6[::-1]\n            data = [fit_pol_list_1[-1], fit_pol_list_2[-1], fit_pol_list_3[-1], fit_pol_list_4[-1], fit_pol_list_5[-1], fit_pol_list_6[-1]]\n            data_list.append(data)\n        data_list = np.asarray(data_list)\n        return data_list\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 15)\n    data_list =  get_results_fit_on_sections_data(results_fit_on_sections)\n\n    diff_4 = np.abs(data_list[:,2,0]-data_list[:,1,0])\n    diff_5 = np.abs(data_list[:,4,0]-data_list[:,3,0])\n    diff_6 = np.abs(data_list[:,2,1]-data_list[:,1,1])\n    diff_7 = np.abs(data_list[:,4,1]-data_list[:,3,1])\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 15)\n    diff = np.std(train[:, :190], axis = 1)\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 40)\n    factor = 0.25\n    factor_3 = 0.25\n    factor_2 =1.6\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = np.mean(train[:, :220], axis = 1, keepdims = True)+train*0\n    pred = pred+(train-np.mean(train, axis = 1, keepdims = True))*0.14\n\n    factor_x = 0.9\n    factor_x_2 = ((np.asarray(range(283))[::-1]/283)[None, :])**5\n    pred = (pred+train*factor_x*factor_x_2)/(1+factor_x*factor_x_2)\n\n    pred_0 = pred.copy()\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 25)\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = np.mean(train[:, :220], axis = 1, keepdims = True)+train*0\n    pred = pred+(train-np.mean(train, axis = 1, keepdims = True))*0.17\n    factor_x = 0.9\n    factor_x_2 = ((np.asarray(range(283))[::-1]/283)[None, :])**2.5\n    pred = (pred+train*factor_x*factor_x_2)/(1+factor_x*factor_x_2)\n    indices = (np.expand_dims(diff, axis = 1)+train*0>0.00005)\n    pred = np.where(indices, pred, pred_0)\n    pred_x = pred_0.copy()\n    pred[:, 200:] = pred_0[:, 200:]\n    pred_0 = pred.copy()\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 30)\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = train.copy()\n    pred = np.where(np.expand_dims(diff, axis = 1)+train*0>0.000071, pred, pred_0)\n    pred[:, 220:] = pred_0[:, 220:]\n    pred_0 = pred.copy()\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 25)\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = train.copy()\n    pred = np.where(np.expand_dims(diff, axis = 1)+train*0>0.00008, pred, pred_0)\n    pred[:, 220:] = pred_0[:, 220:]\n    pred_0 = pred.copy()\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 20)\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = train.copy()\n    indices = (np.expand_dims(diff, axis = 1)+train*0>0.000071)\n    pred = np.where(np.expand_dims(diff, axis = 1)+train*0>0.0001, pred, pred_0)\n    pred[:, 220:] = pred_0[:, 220:]\n    pred_0 = pred.copy()\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 15)\n    train_X = (factor_2*train_1/diff_4[:, None]**factor+train_2/diff_5[:, None]**factor)/(factor_2/diff_4[:, None]**factor+1/diff_5[:, None]**factor)\n    train_Y = (factor_2*train_1/diff_6[:, None]**factor_3+train_2/diff_7[:, None]**factor_3)/(factor_2/diff_6[:, None]**factor_3+1/diff_7[:, None]**factor_3)\n    train = (train_X+train_Y)/2\n    pred = train.copy()\n    pred = np.where(np.expand_dims(diff, axis = 1)+train*0>0.00015, pred, pred_0)\n    pred[:, 220:] = pred_0[:, 220:]\n    pred_0 = pred.copy()\n\n    if test:\n        popt = pickle.load(open(f'{data_path}/popt_model_1.p', 'br'))\n        fit_pol_1 = pickle.load(open(f'{data_path}/fit_pol_1_model_1.p', 'br'))\n\n    if not test:\n        fit_pol_1 = np.polynomial.polynomial.Polynomial.fit(np.mean(pred[:, :220], axis = 1), np.mean(targets-pred, axis = 1), 1)\n        pickle.dump(fit_pol_1, open('fit_pol_1_model_1.p', 'bw'))\n    pred = pred+np.expand_dims(fit_pol_1(np.mean(pred, axis = 1)), 1)\n\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 18)\n    train_sigma = train[:, :180].copy()\n    diff = np.std(train_sigma, axis = 1)\n    train_sigma_2 = np.abs(train_1[:, :185].copy()-train_2[:, :185].copy())\n    diff_2 = np.mean(train_sigma_2**2, axis = 1)**0.5\n\n    x = diff\n    y = diff_2\n\n    def func(xy, a, b, c): \n        x, y = xy \n        return a + b*x + c*y\n\n    if not test:\n        z = z=np.mean(((pred-targets)**2), axis = 1)**0.5\n        data = np.array([x, y, z]).T \n        popt, pcov = curve_fit(func, (x, y), z)\n        pickle.dump(popt, open('popt_model_1.p', 'bw'))\n\n    sigmas_fit_pred = func((x, y), *popt)\n    sigma_pred_final = np.tile(sigmas_fit_pred[:, None], (1, 283))*1.1\n\n    return pred, sigma_pred_final","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.876848Z","iopub.execute_input":"2024-11-03T14:45:27.877315Z","iopub.status.idle":"2024-11-03T14:45:27.889854Z","shell.execute_reply.started":"2024-11-03T14:45:27.877286Z","shell.execute_reply":"2024-11-03T14:45:27.888724Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('model_1_postprocessing.py', 'r').read())\npred, sigma_pred_final = model_1_postprocessing(relred1, relred2, results_fit_on_sections, targets = targets_orig.copy(), test = False, data_path = None)\nsub_df = postprocessing(pred,\n                        adc_info.index,\n                            sigma_pred=sigma_pred_final)\ngll_score = competition_score(train_labels.copy().reset_index(),\n                              sub_df.copy().reset_index(),\n                              naive_mean=naive_mean,\n                              naive_sigma=naive_sigma,\n                              sigma_true=0.00001)\nprint(gll_score[0])\ngll_score = competition_score(train_labels.copy().reset_index()[star == 0],\n      sub_df.copy().reset_index()[star == 0],\n                          naive_mean=naive_mean,\n                          naive_sigma=naive_sigma,\n                          sigma_true=0.00001)\nprint(gll_score[0])\n\ngll_score = competition_score(train_labels.copy().reset_index()[star == 1],\n                          sub_df.copy().reset_index()[star == 1],\n                          naive_mean=naive_mean,\n                          naive_sigma=naive_sigma,\n                          sigma_true=0.00001)\nprint(gll_score[0])","metadata":{"execution":{"iopub.status.busy":"2024-11-03T14:45:27.892815Z","iopub.execute_input":"2024-11-03T14:45:27.893255Z","iopub.status.idle":"2024-11-03T14:45:28.49432Z","shell.execute_reply.started":"2024-11-03T14:45:27.893219Z","shell.execute_reply":"2024-11-03T14:45:28.493117Z"},"trusted":true},"outputs":[],"execution_count":null}]}