{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":205066650,"sourceType":"kernelVersion"},{"sourceId":205193436,"sourceType":"kernelVersion"}],"dockerImageVersionId":30786,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport pickle\nimport multiprocessing\nimport math\nfrom scipy.optimize import curve_fit\nimport scipy\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-29T19:05:31.144308Z","iopub.execute_input":"2024-10-29T19:05:31.144715Z","iopub.status.idle":"2024-10-29T19:05:32.944533Z","shell.execute_reply.started":"2024-10-29T19:05:31.144676Z","shell.execute_reply":"2024-10-29T19:05:32.9433Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_labels_orig = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',\n                        index_col='planet_id')\ntrain_labels = train_labels_orig.to_numpy()\nstar = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv').star.to_numpy()\npoints_list = pickle.load(open('/kaggle/input/adc24-preprocessing/points_list.p', 'br'))\nadvanced_fitting_results = pickle.load(open('/kaggle/input/adc24-preprocessing/advanced_fitting_results.p', 'br'))\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv',\n                           index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:32.946333Z","iopub.execute_input":"2024-10-29T19:05:32.946827Z","iopub.status.idle":"2024-10-29T19:05:34.726693Z","shell.execute_reply.started":"2024-10-29T19:05:32.94679Z","shell.execute_reply":"2024-10-29T19:05:34.725514Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"new_train_data_1 = pickle.load(open('/kaggle/input/adc24-model-2/new_train_data_1.p', 'br'))\nnew_train_data_2 = pickle.load(open('/kaggle/input/adc24-model-2/new_train_data_2.p', 'br'))","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:34.728586Z","iopub.execute_input":"2024-10-29T19:05:34.728977Z","iopub.status.idle":"2024-10-29T19:05:34.761583Z","shell.execute_reply.started":"2024-10-29T19:05:34.72894Z","shell.execute_reply":"2024-10-29T19:05:34.760428Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"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\n\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\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\n    train = (train_1+train_2)/2\n    return train_1, train_2, train","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:34.763795Z","iopub.execute_input":"2024-10-29T19:05:34.764702Z","iopub.status.idle":"2024-10-29T19:05:34.772791Z","shell.execute_reply.started":"2024-10-29T19:05:34.764661Z","shell.execute_reply":"2024-10-29T19:05:34.771579Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def 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)\n\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)\n\n\nnaive_mean = np.mean(train_labels)\nnaive_sigma = np.std(train_labels)","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:34.774701Z","iopub.execute_input":"2024-10-29T19:05:34.775206Z","iopub.status.idle":"2024-10-29T19:05:34.796674Z","shell.execute_reply.started":"2024-10-29T19:05:34.775133Z","shell.execute_reply":"2024-10-29T19:05:34.795508Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile model_2_postprocessing.py\ndef model_2_postprocessing(relred1, relred2, 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, 8)\n    train = (2.5*train_1+train_2)/3.5\n\n    diff = np.std(train[:, :185], axis = 1)\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 25)\n    train = (2.5*train_1+train_2)/3.5\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    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    pred_0 = pred.copy()\n\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 5)\n    train = (2.5*train_1+train_2)/3.5\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\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 10)\n    train = (2.5*train_1+train_2)/3.5\n    pred = train.copy()\n    pred = np.where(np.expand_dims(diff, axis = 1)+train*0>0.00007, pred, pred_0)\n    pred[:, 220:] = pred_0[:, 220:]\n    pred_0 = pred.copy()\n\n\n    train_1, train_2, train = get_freq_averaged(relred1, relred2, 1)\n    train = (2.5*train_1+train_2)/3.5\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\n    pred[:, 220:] = pred_0[:, 220:]\n\n    if test:\n        fit_pol_1 = pickle.load(open(f'{data_path}/fit_pol_1.p', 'br'))\n        popt = pickle.load(open(f'{data_path}/popt.p', 'br'))\n\n    if not test:\n        fit_pol_1 = np.polynomial.polynomial.Polynomial.fit(np.mean(pred, axis = 1), np.mean(targets-pred, axis = 1), 1)\n        pickle.dump(fit_pol_1, open('fit_pol_1.p', 'bw'))\n    pred = pred+np.expand_dims(fit_pol_1(np.mean(pred, axis = 1)), 1)\n\n    train_sigma = train[:, :185].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.p', 'bw'))\n\n    # Perform curve fitting\n    sigmas_fit_pred = func((x, y), *popt)\n    sigma_pred_final = np.tile(sigmas_fit_pred[:, None], (1, 283))+0.0000025\n\n    return pred, sigma_pred_final","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:34.798789Z","iopub.execute_input":"2024-10-29T19:05:34.799276Z","iopub.status.idle":"2024-10-29T19:05:34.811839Z","shell.execute_reply.started":"2024-10-29T19:05:34.799227Z","shell.execute_reply":"2024-10-29T19:05:34.810514Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"relred1 = new_train_data_1.copy()\nrelred2 = new_train_data_2.copy()","metadata":{"execution":{"iopub.status.busy":"2024-10-29T19:05:34.815404Z","iopub.execute_input":"2024-10-29T19:05:34.815761Z","iopub.status.idle":"2024-10-29T19:05:34.825246Z","shell.execute_reply.started":"2024-10-29T19:05:34.815726Z","shell.execute_reply":"2024-10-29T19:05:34.823972Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('model_2_postprocessing.py', 'r').read())\npred, sigma_pred_final = model_2_postprocessing(relred1, relred2, targets = train_labels_orig.to_numpy().copy(), test = False, data_path = None)\n\nsub_df = postprocessing(pred,\n                            train_adc_info.index,\n                                sigma_pred=sigma_pred_final)\n\ngll_score = competition_score(train_labels_orig.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_orig.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_orig.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-10-29T19:05:34.827053Z","iopub.execute_input":"2024-10-29T19:05:34.827529Z","iopub.status.idle":"2024-10-29T19:05:35.866924Z","shell.execute_reply.started":"2024-10-29T19:05:34.82748Z","shell.execute_reply":"2024-10-29T19:05:35.865642Z"},"trusted":true},"outputs":[],"execution_count":null}]}