{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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":"nvidiaTeslaT4","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":9198335,"sourceType":"datasetVersion","datasetId":5561110}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\nimport os\nimport torch\nimport torch.nn as nn\nfrom torch.utils.data import Dataset, DataLoader\nimport torch.optim as optim\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import RidgeCV\nfrom sklearn.metrics import r2_score, mean_squared_error\nimport random\nfrom astropy.stats import sigma_clip\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:07.625083Z","iopub.execute_input":"2024-08-21T19:38:07.625505Z","iopub.status.idle":"2024-08-21T19:38:07.63279Z","shell.execute_reply.started":"2024-08-21T19:38:07.625473Z","shell.execute_reply":"2024-08-21T19:38:07.631754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data\n\nWe start by reading the metadata:","metadata":{}},{"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')\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\nsample_submission = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv',\n                           index_col='planet_id')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n# wavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\n# axis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:07.634644Z","iopub.execute_input":"2024-08-21T19:38:07.635014Z","iopub.status.idle":"2024-08-21T19:38:07.744692Z","shell.execute_reply.started":"2024-08-21T19:38:07.634986Z","shell.execute_reply":"2024-08-21T19:38:07.743724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Some facts about the data:\n- We have 673 planets for training.\n- There are roughly 800 planets for testing (but the test data is hidden).\n- The competition is a multi-output regression task with 283 wavelengths to predict. The first one is from FGS1, the other 282 are from AIRS.","metadata":{}},{"cell_type":"code","source":"assert (train_adc_info.index == train_labels.index).all()\nassert (test_adc_info.index == sample_submission.index).all()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:07.745793Z","iopub.execute_input":"2024-08-21T19:38:07.746107Z","iopub.status.idle":"2024-08-21T19:38:07.751084Z","shell.execute_reply.started":"2024-08-21T19:38:07.746083Z","shell.execute_reply":"2024-08-21T19:38:07.750068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The regression target\n\n`train_labels` contains the training targets for the regression task. There is one row per planet, and the values are the sizes of the planets. You may ask two questions now:\n1. Question: What is the unit of these values? Answer: The value is the square of the radius of the planet divided by the radius of the star \\\\((\\frac{r}{R})^2\\\\). As planets always are much smaller than their host stars, the values are positive but much smaller than 1.\n2. Question: Why are there 283 values per row? Does the planet have 283 different sizes? Answer: We can measure the size of the planet by looking at it at different wavelengths. If the planet's atmosphere is intransparent for the specific wavelength, the planet looks slightly larger.\n","metadata":{}},{"cell_type":"code","source":"train_labels","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:07.753256Z","iopub.execute_input":"2024-08-21T19:38:07.753529Z","iopub.status.idle":"2024-08-21T19:38:07.793477Z","shell.execute_reply.started":"2024-08-21T19:38:07.753505Z","shell.execute_reply":"2024-08-21T19:38:07.792679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(6, 2))\nplt.title(\"Histogram of the planets' sizes (regression targets)\", fontsize=18)\nplt.hist(train_labels.values.ravel(), bins=20, density=True, color='olive')\nplt.xlabel(r\"Planet's size $(\\frac{r}{R})^2$\", fontsize=14)\nplt.ylabel('Density', fontsize=14)\nplt.xlim(0, 0.008)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:07.794575Z","iopub.execute_input":"2024-08-21T19:38:07.794911Z","iopub.status.idle":"2024-08-21T19:38:08.091418Z","shell.execute_reply.started":"2024-08-21T19:38:07.79488Z","shell.execute_reply":"2024-08-21T19:38:08.090466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reading and preprocessing the FGS1 data\n\nThe FGS1 measurements consist of one file per planet (673 files for 673 planets for training). For now, we ignore the calibration files.\n\nEach file contains 135,000 rows of images taken at 0.1 second time steps. Each row is a 32\\*32 image at a single wavelength.\n\nWe read a sample file:","metadata":{}},{"cell_type":"code","source":"planet_id = 14485303\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\nf_signal","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:08.092647Z","iopub.execute_input":"2024-08-21T19:38:08.092974Z","iopub.status.idle":"2024-08-21T19:38:08.556115Z","shell.execute_reply.started":"2024-08-21T19:38:08.092948Z","shell.execute_reply":"2024-08-21T19:38:08.554747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Every row of the file corresponds to an image of a star:","metadata":{}},{"cell_type":"code","source":"sns.heatmap(f_signal.iloc[1].values.reshape(32, 32))\nplt.gca().set_aspect('equal')\nplt.title('FGS1 image')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:08.557565Z","iopub.execute_input":"2024-08-21T19:38:08.55788Z","iopub.status.idle":"2024-08-21T19:38:08.948864Z","shell.execute_reply.started":"2024-08-21T19:38:08.557854Z","shell.execute_reply":"2024-08-21T19:38:08.947878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To see the so-called *light curve*, a time series, we first have to compute the difference between the even and the odd frames to get the net signal. The net signal is very noisy, and we smoothen it by computing a moving average. The plot of the smoothened signal clearly shows that the image gets darker while the planet passes in front of the star (between time steps 22000 and 45000).","metadata":{}},{"cell_type":"code","source":"mean_signal = f_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=800\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n_, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\nax1.plot(net_signal, label='raw net signal')\nax1.legend()\nax2.plot(smooth_signal, color='c', label='smoothened net signal')\nax2.legend()\nax2.set_xlabel('time')\nplt.suptitle('FGS1 light curve', y=0.96)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:08.950114Z","iopub.execute_input":"2024-08-21T19:38:08.950555Z","iopub.status.idle":"2024-08-21T19:38:09.904372Z","shell.execute_reply.started":"2024-08-21T19:38:08.950522Z","shell.execute_reply":"2024-08-21T19:38:09.903459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now read the FGS1 data for all 673 training planets. We keep only three values for every planet: The mean signal before, during and after the planet passes in front of the star.","metadata":{}},{"cell_type":"code","source":"%%time\ndef f_read_and_preprocess(dataset, adc_info):\n    \"\"\"Read the FGS1 files for all planet_ids and extract the signal.\n    \n    Parameters\n    dataset: 'train' or 'test'\n    adc_info: metadata dataframe, either train_adc_info or test_adc_info\n    \n    Returns\n    dataframe with one row per planet_id\n    \n    \"\"\"\n    planet_ids = adc_info.index\n    phases = []\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_signal.parquet')\n        mean_signal = f_signal.values.mean(axis=1) # mean over the 32*32 pixels\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        gain = adc_info.FGS1_adc_gain.values[i]\n        step = 8000\n        phases.append([net_signal[i*step:(i+1)*step].mean() * gain for i in range(len(net_signal) // step + 1)])\n        \n    return pd.DataFrame(phases,\n                        columns=[f\"phase_{i}\" for i in range(len(net_signal) // step + 1)],\n                        index=planet_ids)\n    \n#train = f_read_and_preprocess('train', train_adc_info)","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:09.905516Z","iopub.execute_input":"2024-08-21T19:38:09.905795Z","iopub.status.idle":"2024-08-21T19:38:09.914384Z","shell.execute_reply.started":"2024-08-21T19:38:09.90577Z","shell.execute_reply":"2024-08-21T19:38:09.913529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndef a_read_and_preprocess(dataset, train_adc_info):\n    planet_ids = train_adc_info.index\n    phases = []\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet')\n        mean_signal = a_signal.values.reshape(11250, 32, 356).mean(axis=1).mean(axis=0)\n        phases.append(mean_signal)\n    return phases\n#train_a = a_read_and_preprocess('train', train_adc_info)","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:09.918051Z","iopub.execute_input":"2024-08-21T19:38:09.918366Z","iopub.status.idle":"2024-08-21T19:38:09.930388Z","shell.execute_reply.started":"2024-08-21T19:38:09.918341Z","shell.execute_reply":"2024-08-21T19:38:09.929573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train_aa = pd.DataFrame(train_a, index = train_adc_info.index)","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:09.931457Z","iopub.execute_input":"2024-08-21T19:38:09.932048Z","iopub.status.idle":"2024-08-21T19:38:09.939669Z","shell.execute_reply.started":"2024-08-21T19:38:09.932017Z","shell.execute_reply":"2024-08-21T19:38:09.938848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_final = pd.read_csv(\"/kaggle/input/final-df/final_df.csv\").set_index('planet_id')#\n#train_final = pd.merge(train, train_aa, left_index=True, right_index=True)\ntrain_final","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:09.940745Z","iopub.execute_input":"2024-08-21T19:38:09.941094Z","iopub.status.idle":"2024-08-21T19:38:10.04278Z","shell.execute_reply.started":"2024-08-21T19:38:09.941064Z","shell.execute_reply":"2024-08-21T19:38:10.041895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train_final.to_csv(\"final_df.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:10.044118Z","iopub.execute_input":"2024-08-21T19:38:10.0448Z","iopub.status.idle":"2024-08-21T19:38:10.049001Z","shell.execute_reply.started":"2024-08-21T19:38:10.044767Z","shell.execute_reply":"2024-08-21T19:38:10.047809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following hidden cell contains the function which evaluates the competition metric.","metadata":{}},{"cell_type":"markdown","source":"# Image processing","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{random.choice(train_adc_info.index)}/AIRS-CH0_signal.parquet')\na_signal = a_signal.values.reshape(11250, 32, 356)\nplt.figure(figsize=(10, 3))\nplt.title('AIRS-CH0 image')\nsns.heatmap(a_signal[5])\nplt.ylabel('spatial dimension')\nplt.xlabel('wavelength dimension')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:10.050272Z","iopub.execute_input":"2024-08-21T19:38:10.050864Z","iopub.status.idle":"2024-08-21T19:38:12.99273Z","shell.execute_reply.started":"2024-08-21T19:38:10.050817Z","shell.execute_reply":"2024-08-21T19:38:12.99167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{random.choice(train_adc_info.index)}/FGS1_signal.parquet')\na_signal = a_signal.values.reshape(135000, 32, 32)\nplt.figure(figsize=(10, 3))\nplt.title('AIRS-CH0 image')\nsns.heatmap(a_signal[5])\nplt.ylabel('spatial dimension')\nplt.xlabel('wavelength dimension')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:12.994169Z","iopub.execute_input":"2024-08-21T19:38:12.994572Z","iopub.status.idle":"2024-08-21T19:38:15.037455Z","shell.execute_reply.started":"2024-08-21T19:38:12.994533Z","shell.execute_reply":"2024-08-21T19:38:15.036421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_signal = a_signal.mean(axis=2).mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=80\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n_, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\nax1.plot(net_signal, label='raw signal')\nax1.legend()\nax2.plot(smooth_signal, color='c', label='smoothened signal')\nax2.legend()\nax2.set_xlabel('time')\nplt.suptitle('AIRS-CH0 light curve', y=0.96)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:15.038958Z","iopub.execute_input":"2024-08-21T19:38:15.039358Z","iopub.status.idle":"2024-08-21T19:38:16.063988Z","shell.execute_reply.started":"2024-08-21T19:38:15.039321Z","shell.execute_reply":"2024-08-21T19:38:16.062972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = torch.tensor(train_final.values, dtype=torch.float32)\ny = torch.tensor(train_labels.values, dtype=torch.float32)\n\n# Dataset and DataLoader\nclass PlanetaryDataset(Dataset):\n    def __init__(self, X, y):\n        self.X = X\n        self.y = y\n\n    def __len__(self):\n        return len(self.X)\n\n    def __getitem__(self, idx):\n        return self.X[idx], self.y[idx]\n\n# Creating dataset and dataloaders\ndataset = PlanetaryDataset(X, y)\ntrain_loader = DataLoader(dataset, batch_size=32, shuffle=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:16.065429Z","iopub.execute_input":"2024-08-21T19:38:16.0658Z","iopub.status.idle":"2024-08-21T19:38:16.075204Z","shell.execute_reply.started":"2024-08-21T19:38:16.065766Z","shell.execute_reply":"2024-08-21T19:38:16.074112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#https://github.com/ehsanasgari/pytorch-seq2seq/blob/master/lstm_attention.py\nclass EncoderMLP(nn.Module):\n    def __init__(self, input_size, hidden_size):\n        super(EncoderMLP, self).__init__()\n        self.fc1 = nn.Linear(input_size, hidden_size)\n        self.fc2 = nn.Linear(hidden_size, hidden_size)\n        self.relu = nn.ReLU()\n\n    def forward(self, x):\n        x = self.relu(self.fc1(x))\n        x = self.relu(self.fc2(x))\n        return x\n\nclass RegressionDecoder(nn.Module):\n    def __init__(self, encoder_hidden_dim, output_dim):\n        super(RegressionDecoder, self).__init__()\n        self.fc1 = nn.Linear(encoder_hidden_dim, encoder_hidden_dim)\n        self.fc2 = nn.Linear(encoder_hidden_dim, output_dim)\n        self.relu = nn.ReLU()\n\n    def forward(self, encoded):\n        x = self.relu(self.fc1(encoded))\n        output = self.fc2(x)\n        return output\n\nclass RegressionModel(nn.Module):\n    def __init__(self, input_size, encoder_hidden_size, output_dim):\n        super(RegressionModel, self).__init__()\n        self.encoder = EncoderMLP(input_size, encoder_hidden_size)\n        self.decoder = RegressionDecoder(encoder_hidden_size, output_dim)\n\n    def forward(self, x):\n        encoded = self.encoder(x)\n        output = self.decoder(encoded)\n        return output","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:16.076825Z","iopub.execute_input":"2024-08-21T19:38:16.077158Z","iopub.status.idle":"2024-08-21T19:38:16.088697Z","shell.execute_reply.started":"2024-08-21T19:38:16.077132Z","shell.execute_reply":"2024-08-21T19:38:16.087788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels.shape[1]","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:16.08983Z","iopub.execute_input":"2024-08-21T19:38:16.090122Z","iopub.status.idle":"2024-08-21T19:38:16.105396Z","shell.execute_reply.started":"2024-08-21T19:38:16.090097Z","shell.execute_reply":"2024-08-21T19:38:16.104464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Hyperparameters\ninput_size = train_final.shape[1]  # Number of input features (e.g., 365)\nencoder_hidden_size = 128  # Hidden size in the encoder\noutput_dim = train_labels.shape[1]  # Number of output wavelengths (e.g., 283)\nlearning_rate = 0.001\nnum_epochs = 100\n\n# Instantiate the model, loss function, and optimizer\nmodel = RegressionModel(input_size, encoder_hidden_size, output_dim)\ncriterion = nn.MSELoss()  # or nn.L1Loss()\noptimizer = optim.Adam(model.parameters(), lr=learning_rate)\n\n# Training loop\nfor epoch in range(num_epochs):\n    model.train()\n    running_loss = 0.0\n    for inputs, targets in train_loader:\n        optimizer.zero_grad()\n        outputs = model(inputs)\n        loss = criterion(outputs, targets)\n        loss.backward()\n        optimizer.step()\n        running_loss += loss.item()\n    if (epoch + 1) % 20 == 0:\n        print(f'Epoch [{epoch+1}/{num_epochs}], Loss: {loss.item():.4f}')","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:16.106553Z","iopub.execute_input":"2024-08-21T19:38:16.106852Z","iopub.status.idle":"2024-08-21T19:38:23.656061Z","shell.execute_reply.started":"2024-08-21T19:38:16.106807Z","shell.execute_reply":"2024-08-21T19:38:23.655084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"torch.save(model.state_dict(), 'regression_model.pth')\n# Instantiate the model again with the same architecture\nmodel = RegressionModel(input_size, encoder_hidden_size, output_dim)\n\n# Load the saved model weights\nmodel.load_state_dict(torch.load('regression_model.pth'))","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:23.657768Z","iopub.execute_input":"2024-08-21T19:38:23.658168Z","iopub.status.idle":"2024-08-21T19:38:23.670483Z","shell.execute_reply.started":"2024-08-21T19:38:23.658133Z","shell.execute_reply":"2024-08-21T19:38:23.669408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_a = a_read_and_preprocess('test', test_adc_info)\ntest_a = pd.DataFrame(test_a).reset_index(drop = True)\ntest = f_read_and_preprocess('test', test_adc_info).reset_index(drop = True)\ntest_final = pd.merge(test, test_a, left_index=True, right_index=True)\ntest_final","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:23.67159Z","iopub.execute_input":"2024-08-21T19:38:23.671894Z","iopub.status.idle":"2024-08-21T19:38:25.341966Z","shell.execute_reply.started":"2024-08-21T19:38:23.671864Z","shell.execute_reply":"2024-08-21T19:38:25.34107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"to_preds = torch.tensor(test_final.values, dtype=torch.float32)\nto_preds.shape","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.34343Z","iopub.execute_input":"2024-08-21T19:38:25.343848Z","iopub.status.idle":"2024-08-21T19:38:25.351686Z","shell.execute_reply.started":"2024-08-21T19:38:25.343795Z","shell.execute_reply":"2024-08-21T19:38:25.35067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Calculate Sigma\nmodel.eval()\nwith torch.no_grad():\n    predictions_all = model(torch.tensor(train_final.values, dtype=torch.float32)).numpy()\nsigma_p = mean_squared_error(train_labels, predictions_all, squared=False)\nsigma_p","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.352859Z","iopub.execute_input":"2024-08-21T19:38:25.353136Z","iopub.status.idle":"2024-08-21T19:38:25.375756Z","shell.execute_reply.started":"2024-08-21T19:38:25.353112Z","shell.execute_reply":"2024-08-21T19:38:25.374883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate predictions\nmodel.eval()\nwith torch.no_grad():\n    predictions = model(to_preds).numpy()\npredictions = predictions.clip(train_labels.values.min(), train_labels.values.max())","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.376982Z","iopub.execute_input":"2024-08-21T19:38:25.377348Z","iopub.status.idle":"2024-08-21T19:38:25.386495Z","shell.execute_reply.started":"2024-08-21T19:38:25.377311Z","shell.execute_reply":"2024-08-21T19:38:25.38561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df = pd.DataFrame(predictions)\npred_df.columns = train_labels.columns\n\nsigma_df = pd.DataFrame({x.replace('wl_', 'sigma_'): sigma_p for x in pred_df.columns}, index=pred_df.index)\npred_df = pd.concat([pred_df, sigma_df], axis=1)\ntrain_labelss = pd.concat([train_labels, sigma_df], axis=1)#.rename_axis('planet_id').reset_index()\npred_df.insert(0, 'planet_id', 499191466)\npred_df","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.387957Z","iopub.execute_input":"2024-08-21T19:38:25.388292Z","iopub.status.idle":"2024-08-21T19:38:25.418651Z","shell.execute_reply.started":"2024-08-21T19:38:25.388261Z","shell.execute_reply":"2024-08-21T19:38:25.417883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Adapted from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef competition_score(solution: pd.DataFrame, submission: pd.DataFrame,\n                      naive_mean: float, naive_sigma: float,sigma_true: float,\n                      sigma_p:float\n    ) -> float:\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    print(n_wavelengths)\n    print(len(submission.columns))\n    if len(submission.columns) != n_wavelengths:\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    n, m = 673,566\n    sigma_pred = np.full((n, m), sigma_p)\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    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    return float(np.clip(submit_score, 0.0, 1.0))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-21T19:38:25.419709Z","iopub.execute_input":"2024-08-21T19:38:25.419998Z","iopub.status.idle":"2024-08-21T19:38:25.431131Z","shell.execute_reply.started":"2024-08-21T19:38:25.419976Z","shell.execute_reply":"2024-08-21T19:38:25.430283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score = competition_score(solution=train_labelss[:673].fillna(sigma_p), \n                          submission=pred_df.drop('planet_id',axis=1),#pd.DataFrame(predictions), \n                          naive_mean=train_labels.values.mean(),\n                          naive_sigma=train_labels.values.std(),\n                          sigma_true=0.00001,\n                          sigma_p=sigma_p\n                          )\nprint(f'Competition Score: {score:.4f}')","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.43636Z","iopub.execute_input":"2024-08-21T19:38:25.43664Z","iopub.status.idle":"2024-08-21T19:38:25.50545Z","shell.execute_reply.started":"2024-08-21T19:38:25.436617Z","shell.execute_reply":"2024-08-21T19:38:25.504517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df.shape","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.506567Z","iopub.execute_input":"2024-08-21T19:38:25.506887Z","iopub.status.idle":"2024-08-21T19:38:25.512975Z","shell.execute_reply.started":"2024-08-21T19:38:25.506854Z","shell.execute_reply":"2024-08-21T19:38:25.512037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-08-21T19:38:25.514036Z","iopub.execute_input":"2024-08-21T19:38:25.514325Z","iopub.status.idle":"2024-08-21T19:38:25.526378Z","shell.execute_reply.started":"2024-08-21T19:38:25.514301Z","shell.execute_reply":"2024-08-21T19:38:25.5256Z"},"trusted":true},"execution_count":null,"outputs":[]}]}