{"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"}],"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 matplotlib.pyplot as plt\nimport pandas as pd\nimport scipy","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:54.903069Z","iopub.execute_input":"2024-10-27T16:26:54.903496Z","iopub.status.idle":"2024-10-27T16:26:56.0664Z","shell.execute_reply.started":"2024-10-27T16:26:54.903452Z","shell.execute_reply":"2024-10-27T16:26:56.065208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Real labels\n\nFor testing","metadata":{}},{"cell_type":"code","source":"labels_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv')\nsolution_df = labels_df.reset_index().drop('index',axis=1)\nsolution_df","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.068476Z","iopub.execute_input":"2024-10-27T16:26:56.06902Z","iopub.status.idle":"2024-10-27T16:26:56.237919Z","shell.execute_reply.started":"2024-10-27T16:26:56.068978Z","shell.execute_reply":"2024-10-27T16:26:56.236632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scoring","metadata":{"execution":{"iopub.status.busy":"2024-10-21T21:08:10.728818Z","iopub.execute_input":"2024-10-21T21:08:10.729305Z","iopub.status.idle":"2024-10-21T21:08:10.734728Z","shell.execute_reply.started":"2024-10-21T21:08:10.729256Z","shell.execute_reply":"2024-10-21T21:08:10.733467Z"}}},{"cell_type":"code","source":"class ParticipantVisibleError(Exception):\n    pass\n\n\ndef score(\n        solution: pd.DataFrame,\n        submission: pd.DataFrame,\n        row_id_column_name: str,\n        naive_mean: float,\n        naive_sigma: float,\n        sigma_true: float\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        print(list(solution.columns))\n        raise ParticipantVisibleError(f'Wrong number of columns in the submission {n_wavelengths*2}vs {len(submission.columns)}')\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    return float(np.clip(submit_score, 0.0, 1.0))\n","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.239587Z","iopub.execute_input":"2024-10-27T16:26:56.240137Z","iopub.status.idle":"2024-10-27T16:26:56.258338Z","shell.execute_reply.started":"2024-10-27T16:26:56.240091Z","shell.execute_reply":"2024-10-27T16:26:56.256785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generation","metadata":{}},{"cell_type":"code","source":"def polynomial(x, *c_i):\n    pol = np.sum(p * (x)**(i) for i, p in enumerate(c_i))\n    return pol \n\ndef get_labels(dep_type:str = 'pol', n_planets=10, x_data = None):\n    result = []\n    if x_data is None:\n        x_data = np.linspace(0,1, 282)\n    if dep_type == 'pol':\n        par_limits = {'c0': (4e-3, 5.0e-3), 'c1':(-5.0e-6, 5.0e-6), 'c2':(-5.0e-6, 5.0e-6), 'c3':(-3.e-6, 3.0e-6), \n                      'c4':(-3.e-6, 3.0e-6), 'c5':(-3.e-7, 3.0e-7)}\n        pol_order = len(par_limits)\n        params = {}\n        for par in par_limits:\n            params[par] = np.random.uniform(par_limits[par][0], par_limits[par][1], n_planets)\n        print(params)\n        for i in range(n_planets):\n            c_i = [params[f'c{j}'][i] for j in range(pol_order)]\n            result += [polynomial(x_data, *c_i)]\n    vals = np.array(result)\n    val_cols = [f'wl_{i}' for i in range(vals.shape[1])]\n    return pd.DataFrame(vals, columns=val_cols).reset_index().rename({'index':'planet_id'}, axis=1)\npol_labels = get_labels()","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.262677Z","iopub.execute_input":"2024-10-27T16:26:56.263178Z","iopub.status.idle":"2024-10-27T16:26:56.286642Z","shell.execute_reply.started":"2024-10-27T16:26:56.263132Z","shell.execute_reply":"2024-10-27T16:26:56.285268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plot normalized","metadata":{}},{"cell_type":"code","source":"plt.plot(pol_labels.values[:,1:].T/pol_labels.values[:,1:].mean(axis=1));","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.288192Z","iopub.execute_input":"2024-10-27T16:26:56.288706Z","iopub.status.idle":"2024-10-27T16:26:56.604932Z","shell.execute_reply.started":"2024-10-27T16:26:56.288662Z","shell.execute_reply":"2024-10-27T16:26:56.603632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reconstruction","metadata":{}},{"cell_type":"code","source":"def get_reco_df(labels, bias_func, err_func,  bias_pars: dict=None, err_pars:dict=None):\n    if bias_pars is None:\n        bias_pars = {}\n    if err_pars is None:\n        err_pars = {}\n    reco = bias_func(labels.values[:,1:], **bias_pars)\n    assert reco.shape == labels.values[:,1:].shape\n    errors = err_func(labels.values[:,1:], **err_pars)\n    assert errors.shape == labels.values[:,1:].shape\n\n    values = np.concatenate([reco, errors], axis=1)\n    val_cols = [f'wl_{i}' for i in range(labels.values[:,1:].shape[1])]\n    err_cols = [f'sigma_{i}' for i in range(labels.values[:,1:].shape[1])]\n    df = pd.DataFrame(values, columns=val_cols+err_cols)\n    return df.reset_index().rename({'index':'planet_id'}, axis=1)\n\nget_reco_df(pol_labels, lambda x: (np.ones((x.shape[1], 1))*x.mean(axis=1)).T, lambda x:(np.ones((x.shape[1], 1))*x.std(axis=1)).T ) ","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.606508Z","iopub.execute_input":"2024-10-27T16:26:56.606978Z","iopub.status.idle":"2024-10-27T16:26:56.663656Z","shell.execute_reply.started":"2024-10-27T16:26:56.606927Z","shell.execute_reply":"2024-10-27T16:26:56.662378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mean value only, no value bias, linear scan by errors","metadata":{}},{"cell_type":"code","source":"x_vals = np.linspace(0.000001,0.001, 100)\nresult = []\nfor err_factor in x_vals:\n    pol_reco_df = get_reco_df(pol_labels, \n                               lambda x: (np.ones((x.shape[1], 1))*x.mean(axis=1)).T, \n                               lambda x:np.ones((x.shape))*err_factor) \n    s = score(pd.DataFrame(pol_labels), pd.DataFrame(pol_reco_df),\n              naive_mean=pol_labels.drop(['planet_id'],axis=1).values.mean(), \n              row_id_column_name='planet_id',\n              naive_sigma=pol_labels.drop(['planet_id'],axis=1).values.std(),                              \n              sigma_true=1e-5)\n    result += [s]\nplt.plot(x_vals, result)\nplt.xlabel('Error bias factor')","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:26:56.665328Z","iopub.execute_input":"2024-10-27T16:26:56.665808Z","iopub.status.idle":"2024-10-27T16:27:02.032487Z","shell.execute_reply.started":"2024-10-27T16:26:56.665751Z","shell.execute_reply":"2024-10-27T16:27:02.031329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mean value only, scan via a small bias, errors are constant","metadata":{}},{"cell_type":"code","source":"x_vals = np.linspace(0.99,1.01, 100)\nerr_val = 1e-5\nresult = []\nfor val_factor in x_vals:\n    pol_reco_df = get_reco_df(pol_labels, \n                               lambda x: (np.ones((x.shape[1], 1))*x.mean(axis=1)).T*val_factor, \n                               lambda x:(np.ones((x.shape))*err_val))\n    \n    s = score(pd.DataFrame(pol_labels), pd.DataFrame(pol_reco_df),\n              naive_mean=pol_labels.drop(['planet_id'],axis=1).values.mean(), \n              row_id_column_name='planet_id',\n              naive_sigma=pol_labels.drop(['planet_id'],axis=1).values.std(),                              \n              sigma_true=1e-5)\n    result += [s]\nplt.plot(x_vals, result)\nplt.xlabel('Value bias factor')\nplt.ylabel('Score');","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:02.03435Z","iopub.execute_input":"2024-10-27T16:27:02.034927Z","iopub.status.idle":"2024-10-27T16:27:06.30322Z","shell.execute_reply.started":"2024-10-27T16:27:02.034882Z","shell.execute_reply":"2024-10-27T16:27:06.30219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_vals = np.linspace(0.9,1.1, 100)\nerr_val = 1e-5\nresult = []\nfor val_factor in x_vals:\n    pol_reco_df = get_reco_df(pol_labels, \n                               lambda x: (np.ones((x.shape[1], 1))*x.mean(axis=1)).T*val_factor, \n                               lambda x:(np.ones((x.shape))*err_val))\n    \n    s = score(pd.DataFrame(pol_labels), pd.DataFrame(pol_reco_df),\n              naive_mean=pol_labels.drop(['planet_id'],axis=1).values.mean(), \n              row_id_column_name='planet_id',\n              naive_sigma=pol_labels.drop(['planet_id'],axis=1).values.std(),                              \n              sigma_true=1e-5)\n    result += [s]\nplt.plot(x_vals, result)\nplt.xlabel('Value bias factor')\nplt.ylabel('Score');","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:06.304771Z","iopub.execute_input":"2024-10-27T16:27:06.305142Z","iopub.status.idle":"2024-10-27T16:27:10.44173Z","shell.execute_reply.started":"2024-10-27T16:27:06.305102Z","shell.execute_reply":"2024-10-27T16:27:10.440667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Biased value scan, constant error","metadata":{}},{"cell_type":"code","source":"x_vals = np.linspace(-1,1, 51)\nlinestyles = ['-', '--', '-.', '-', ':']\nfor err_val, style in zip(np.linspace(0.00008, 0.00012, 5), linestyles):\n    result = []\n    for val_factor in x_vals:\n        pol_reco_df = get_reco_df(pol_labels, \n                                   lambda x: x+val_factor*1e-4, \n                                   lambda x:(np.ones((x.shape))*err_val))\n\n\n        s = score(pd.DataFrame(pol_labels), pd.DataFrame(pol_reco_df),\n                  naive_mean=pol_labels.drop(['planet_id'],axis=1).values.mean(), \n                  row_id_column_name='planet_id',\n                  naive_sigma=pol_labels.drop(['planet_id'],axis=1).values.std(),                              \n                  sigma_true=1e-5)\n        if val_factor == 1:\n            print(val_factor*1e-4, err_val, s, )\n        result += [s]\n    \n    plt.plot(x_vals, result, label=f'Error {np.round(err_val,5)}', linewidth=2, linestyle=style)\nplt.xlabel('Value bias factor * 1e-4')\nplt.legend()\nplt.grid()\nplt.ylabel('Score');","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:10.445114Z","iopub.execute_input":"2024-10-27T16:27:10.445485Z","iopub.status.idle":"2024-10-27T16:27:20.910958Z","shell.execute_reply.started":"2024-10-27T16:27:10.445445Z","shell.execute_reply":"2024-10-27T16:27:20.909722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_vals","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:20.912257Z","iopub.execute_input":"2024-10-27T16:27:20.912629Z","iopub.status.idle":"2024-10-27T16:27:20.920802Z","shell.execute_reply.started":"2024-10-27T16:27:20.912588Z","shell.execute_reply":"2024-10-27T16:27:20.919625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"5e-4 /3e-5","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:20.922425Z","iopub.execute_input":"2024-10-27T16:27:20.922904Z","iopub.status.idle":"2024-10-27T16:27:20.933715Z","shell.execute_reply.started":"2024-10-27T16:27:20.922849Z","shell.execute_reply":"2024-10-27T16:27:20.932603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Biased value and proportinal error scan","metadata":{}},{"cell_type":"code","source":"x_vals = np.linspace(-0.0002,0.0002, 100)\nresult = []\nfor val_factor in x_vals:\n    pol_reco_df = get_reco_df(pol_labels, \n                               lambda x: x+1e-4, \n                               lambda x:(np.ones((x.shape))*abs(val_factor)))\n    \n    \n    s = score(pd.DataFrame(pol_labels), pd.DataFrame(pol_reco_df),\n              naive_mean=pol_labels.drop(['planet_id'],axis=1).values.mean(), \n              row_id_column_name='planet_id',\n              naive_sigma=pol_labels.drop(['planet_id'],axis=1).values.std(),                              \n              sigma_true=1e-5)\n    result += [s]\nplt.plot(x_vals, result)\nplt.xlabel('Value bias factor')\nplt.ylabel('Score');","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:20.935179Z","iopub.execute_input":"2024-10-27T16:27:20.936041Z","iopub.status.idle":"2024-10-27T16:27:25.234737Z","shell.execute_reply.started":"2024-10-27T16:27:20.93598Z","shell.execute_reply":"2024-10-27T16:27:25.233535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pol_reco_df","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:25.236422Z","iopub.execute_input":"2024-10-27T16:27:25.236937Z","iopub.status.idle":"2024-10-27T16:27:25.272259Z","shell.execute_reply.started":"2024-10-27T16:27:25.23688Z","shell.execute_reply":"2024-10-27T16:27:25.271128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pol_labels","metadata":{"execution":{"iopub.status.busy":"2024-10-27T16:27:25.273682Z","iopub.execute_input":"2024-10-27T16:27:25.274071Z","iopub.status.idle":"2024-10-27T16:27:25.307972Z","shell.execute_reply.started":"2024-10-27T16:27:25.27403Z","shell.execute_reply":"2024-10-27T16:27:25.306872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}