{"metadata":{"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":191936343,"sourceType":"kernelVersion"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true},"kernelspec":{"display_name":"Python 3 (ipykernel)","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.12.0"}},"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\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom fastprogress import master_bar, progress_bar\nfrom torch.optim import Adam\nfrom torch.optim.lr_scheduler import CosineAnnealingLR\nfrom torch.utils.data import Dataset, DataLoader\nfrom torchvision.transforms import transforms\nimport torchvision.models as models\nfrom torchvision.models.efficientnet import _efficientnet_conf, _efficientnet\nfrom functools import partial\nimport random, os\n\ndef seed_everything(seed: int):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = False\n    \nseed_everything(42)","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metric","metadata":{}},{"cell_type":"code","source":"class ParticipantVisibleError(Exception):\n    pass\ndef ariel_score(\n        solution,\n        submission,\n        naive_mean,\n        naive_sigma,\n        sigma_true\n    ):\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    if submission.min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n\n    n_wavelengths = 283\n\n    y_pred = submission[:, :n_wavelengths]\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission[:, n_wavelengths:], a_min=10**-15, a_max=None)\n    y_true = solution\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    #print(GLL_pred, GLL_true, GLL_mean)\n    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    return submit_score #float(np.clip(submit_score, 0.0, 1.0))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Raw data binned 15 time measurements\n\nUnfortunately, I have lost the exact code that I received the dataset with, I really hope that it matches the one in the inference notebook.","metadata":{}},{"cell_type":"code","source":"train_adc_info = pd.read_csv('./train_adc_info.csv',\n                           index_col='planet_id')\ntrain_labels = pd.read_csv('./train_labels.csv',\n                           index_col='planet_id')\n\npre_train = np.load('./ariel_calibrated_train.npy')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pre_train.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Transit phase detector\n\nA simple algorithm for finding points where there is a rapid drop in luminosity","metadata":{}},{"cell_type":"code","source":"def phase_detector(signal):\n    phase1, phase2 = None, None\n    best_drop = 0\n    for i in range(50,150):        \n        t1 = signal[i:i+20].max() - signal[i:i+20].min()\n        if t1 > best_drop:\n            phase1 = i+20+5\n            best_drop = t1\n    \n    best_drop = 0\n    for i in range(200,300):\n        t1 = signal[i:i+20].max() - signal[i:i+20].min()\n        if t1 > best_drop:\n            phase2 = i-5\n            best_drop = t1\n    \n    return phase1, phase2","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset\nWe will use (flux - flux_transit) / flux_star as features. If we didn't have a distortion of the star's light, then our targets would be exactly among these signs.","metadata":{}},{"cell_type":"code","source":"train = pre_train.copy()\nfor i in range(len(train_adc_info)):\n    p1,p2 = phase_detector(pre_train[i,:,1:].mean(axis=1))\n    train[i] = (train[i] - pre_train[i,p1:p2].mean(axis=0)) / pre_train[i,list(range(p1-40)) + list(range(p2+40,375))].mean(axis=0) * 1000.0","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n = 205\nplt.plot(pre_train[n,:,1:].mean(axis=1))\np1, p2 = phase_detector(pre_train[n,:,1:].mean(axis=1))\nplt.axvline(p1)\nplt.axvline(p2)\nplt.axvline(p1-40)\nplt.axvline(p2+40)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.mean(), train.std()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model \nAs a model, I use MobileNet_v2 as the least overfitting. Initially, it accepts a 3-channel image as input, so I added a layer that converts the input into a three-channel one.\n\nThe output has 3 numbers for each wavelength corresponding to 0.1587,0.5,0.8413 percentiles.","metadata":{}},{"cell_type":"code","source":"\ndef create_model_mnet2():\n    model = models.mobilenet_v3_small(dropout=0.0, norm_layer = nn.Identity)\n    model.features[0][0] = nn.Conv2d(3, 16, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False)\n    model.classifier[3] = nn.Linear(in_features=1024, out_features=283*3, bias=True)\n    return model\n\nclass ImpModel(torch.nn.Module):\n    def __init__(self):\n        super(ImpModel, self).__init__()\n\n        self.filter = nn.Sequential(\n            nn.Conv2d(1, 3, kernel_size=(3,1), stride = (2,1), bias=False),\n            nn.LeakyReLU()\n        )\n        self.model_1d = create_model_mnet2()\n        \n    def forward(self, x):\n        x = self.filter(x)\n        x = self.model_1d(x)\n        return x","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loss\nThis loss function predicts not only the value itself, but also its two percentiles. With their help, assuming that the error is distributed normally, we can get an estimate of the standard deviation.\n\nhttps://www.kaggle.com/code/vyacheslavefimov/quantile-loss-quantile-regression","metadata":{}},{"cell_type":"code","source":"def q_loss(quantiles, y_pred, target):\n    losses = []\n    for i, q in enumerate(quantiles):        \n        errors = target - y_pred[..., i]\n        losses.append(torch.max((q - 1) * errors, q * errors).unsqueeze(-1))\n    losses = 2 * torch.cat(losses, dim=2)\n\n    return losses\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import KFold\n\nkf = KFold(n_splits=5)\nX = list(range(len(train)))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"oof_pred = np.zeros_like(train_labels.values)\noof_sigmas = np.zeros_like(train_labels.values)\n\nfor ifold, (train_index, test_index) in enumerate(kf.split(X)):\n    train_x = torch.from_numpy(train[train_index]).unsqueeze(1).float()\n    train_y = torch.from_numpy(train_labels.values[train_index]).float()\n    train_dataset = torch.utils.data.TensorDataset(train_x, train_y)\n    training_loader = torch.utils.data.DataLoader(train_dataset, batch_size=16, shuffle=True)\n\n    val_x = torch.from_numpy(train[test_index]).unsqueeze(1).float()\n    val_y = torch.from_numpy(train_labels.values[test_index]).float()\n    val_dataset = torch.utils.data.TensorDataset(val_x, val_y)\n    validation_loader = torch.utils.data.DataLoader(val_dataset, batch_size=16, shuffle=False)\n    \n    model = ImpModel().cuda()\n\n    best_metric = 0\n    total_train_losses = []\n\n    optimizer = torch.optim.AdamW(model.parameters(), lr=0.001)\n    scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, 200, eta_min=0, last_epoch=-1)\n    for epoch in range(200):\n        ep_losses = []        \n        model.train()\n        for i, data in enumerate(training_loader):\n            inputs, labels = data\n            optimizer.zero_grad()\n\n            #augmentations: disturb channel and reverse the time\n            if epoch > 0:\n                for j in range(inputs.shape[-1]):\n                    if np.random.random() > 0.8:\n                        inputs[...,j] = inputs[...,j] * (1 + np.random.randn() * 0.01)\n                for k in range(inputs.shape[0]):\n                    if np.random.random() > 0.8:\n                        inputs[k] = torch.flip(inputs[k], (1,))\n\n            outputs = model(inputs.cuda()).reshape((inputs.shape[0], 283, 3))\n\n            loss = q_loss([0.1587,0.5,0.8413], outputs, labels.cuda()).mean()\n            loss.backward()\n\n            optimizer.step()\n\n            ep_losses.append(loss.item())\n\n        avg_loss = np.mean(ep_losses)\n        total_train_losses.append(avg_loss)\n\n        model.eval()\n        running_vloss = 0\n        preds = np.zeros((len(val_dataset), 283))\n        sigmas = np.zeros((len(val_dataset), 283))\n        v_offset = 0\n        with torch.no_grad():        \n            for i, vdata in enumerate(validation_loader):\n                vinputs, vlabels = vdata\n                voutputs = model(vinputs.cuda()).reshape((vinputs.shape[0], 283, 3))\n                preds[v_offset:v_offset+len(vinputs)] = voutputs[:,:,1].detach().cpu().numpy()\n                sigmas[v_offset:v_offset+len(vinputs)] = voutputs[:,:,2].detach().cpu().numpy() - voutputs[:,:,0].detach().cpu().numpy()\n                vloss = q_loss([0.1587,0.5,0.8413], voutputs, vlabels.cuda()).mean()\n                running_vloss += vloss\n                v_offset += len(vinputs)\n\n        avg_vloss = running_vloss / (i + 1)\n\n        metric = ariel_score(        \n            train_labels.values[test_index],\n            np.concatenate([preds.clip(0), sigmas.clip(0)], axis=1), \n            train_labels.values[test_index].mean(),\n            train_labels.values[test_index].std(),\n            sigma_true=1e-5)\n        \n        print('fold {} epoch {} train {} valid {} metric {}'.format(ifold, epoch, \n                                                                    round(avg_loss*100,3), \n                                                                    round(avg_vloss.item()*100,3), \n                                                                    round(metric,3) ))\n        scheduler.step()\n        \n    oof_pred[test_index] = preds\n    oof_sigmas[test_index] = sigmas\n    torch.save(model.state_dict(), 'purple_hat_{}'.format(ifold))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ariel_score(        \n            train_labels.values,\n            np.concatenate([oof_pred.clip(0), oof_sigmas.clip(0)], axis=1), \n            train_labels.values.mean(),\n            train_labels.values.std(),\n            sigma_true=1e-5)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.save('purple_hat_oof_pred.npy', oof_pred)\nnp.save('purple_hat_oof_sigmas.npy', oof_sigmas)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xs,ys = [], []\nfor i in range(200):\n    xs.append(i/100.0)\n    ys.append(ariel_score(        \n            train_labels.values,\n            np.concatenate([oof_pred.clip(0), oof_sigmas.clip(0) * (i / 100.0)], axis=1), \n            train_labels.values.mean(),\n            train_labels.values.std(),\n            sigma_true=1e-5))\nplt.plot(xs, np.array(ys).clip(0))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}