{"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":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# preface \nI was inspired by @AmbrosM idea https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/530152#2969648\n\nThe idea of ​​this approach is to predict the average answer for each test sample. Instead of building models, we will search for the correct answer experimentally. We will select for each spectrum such a multiplier that the transit part multiplied by it will \"line up\" with the other points. The line can be either a straight line or a polynomial up to the 5th degree. For selection, we will use the Nelder-Mead method. The found multiplication factor minus one is our answer.\n\nThere are some changes in data preparation here.\n* dt for dark calibration changed in favor https://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data/comments#2964759\n* signal clipped to zero due this https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/530247#2970709\n* I masked hot and dead pixels with NaN in flat and averaging through spatial dimension.","metadata":{}},{"cell_type":"markdown","source":"### **upd for version 2**: it seems I forgot to remove **5x sigma** copied from somewhere. I wonder how this will affect the score.","metadata":{}},{"cell_type":"markdown","source":"### preprocess the data","metadata":{}},{"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 itertools\nfrom scipy.optimize import minimize\nfrom functools import partial\nimport random, os\nfrom astropy.stats import sigma_clip'''\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport itertools\nfrom tqdm import tqdm\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom scipy.optimize import minimize\nfrom functools import partial\nfrom astropy.stats import sigma_clip\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import DataLoader, Dataset","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-09-21T02:39:17.669846Z","iopub.execute_input":"2024-09-21T02:39:17.670932Z","iopub.status.idle":"2024-09-21T02:39:17.680403Z","shell.execute_reply.started":"2024-09-21T02:39:17.670884Z","shell.execute_reply":"2024-09-21T02:39:17.678892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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')\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n","metadata":{"execution":{"iopub.status.busy":"2024-09-21T02:39:17.682866Z","iopub.execute_input":"2024-09-21T02:39:17.683968Z","iopub.status.idle":"2024-09-21T02:39:17.737557Z","shell.execute_reply.started":"2024-09-21T02:39:17.683916Z","shell.execute_reply":"2024-09-21T02:39:17.735919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def apply_linear_corr(linear_corr,clean_signal):\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(\n                range(clean_signal.shape[1]), range(clean_signal.shape[2])\n            ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal\n\ndef clean_dark(signal, dark, dt):\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n    signal -= dark* dt[:, np.newaxis, np.newaxis]\n    return signal\n\ndef preproc(dataset, adc_info, sensor, binning = 15):\n    cut_inf, cut_sup = 39, 321\n    sensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\n    binned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\n    linear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\n    planet_ids = adc_info.index\n    \n    feats = []\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_signal.parquet').to_numpy()\n        dark_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/dark.parquet', engine='pyarrow').to_numpy()\n        dead_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/dead.parquet', engine='pyarrow').to_numpy()\n        flat_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/flat.parquet', engine='pyarrow').to_numpy()\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/linear_corr.parquet').values.astype(np.float64).reshape(linear_corr_dict[sensor])\n\n        signal = signal.reshape(sensor_sizes_dict[sensor][0]) \n        gain = adc_info[f'{sensor}_adc_gain'].values[i]\n        offset = adc_info[f'{sensor}_adc_offset'].values[i]\n        signal = signal / gain + offset\n        \n        hot = sigma_clip(\n            dark_frame, sigma=5, maxiters=5\n        ).mask\n        \n        if sensor != \"FGS1\":\n            signal = signal[:, :, cut_inf:cut_sup] #11250 * 32 * 282\n            #dt = axis_info['AIRS-CH0-integration_time'].dropna().values\n            dt = np.ones(len(signal))*0.1 \n            dt[1::2] += 4.5 #@bilzard idea\n            linear_corr = linear_corr[:, :, cut_inf:cut_sup]\n            dark_frame = dark_frame[:, cut_inf:cut_sup]\n            dead_frame = dead_frame[:, cut_inf:cut_sup]\n            flat_frame = flat_frame[:, cut_inf:cut_sup]\n            hot = hot[:, cut_inf:cut_sup]\n        else:\n            dt = np.ones(len(signal))*0.1\n            dt[1::2] += 0.1\n            \n        signal = signal.clip(0) #@graySnow idea\n        linear_corr_signal = apply_linear_corr(linear_corr, signal)\n        signal = clean_dark(linear_corr_signal, dark_frame, dt)\n        \n        flat = flat_frame.reshape(sensor_sizes_dict[sensor][1])\n        flat[dead_frame.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        flat[hot.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        signal = signal / flat\n        \n        if sensor == \"FGS1\":\n            signal = signal.reshape((sensor_sizes_dict[sensor][0][0], sensor_sizes_dict[sensor][0][1]*sensor_sizes_dict[sensor][0][2]))\n        \n        mean_signal = np.nanmean(signal, axis=1) # mean over the 32*32(FGS1) or 32(CH0) pixels\n        cds_signal = (mean_signal[1::2] - mean_signal[0::2])\n        \n        binned = np.zeros((binned_dict[sensor]))\n        for j in range(cds_signal.shape[0] // binning):\n            binned[j] = cds_signal[j*binning:j*binning+binning].mean(axis=0)\n                   \n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0],1))\n            \n        feats.append(binned)\n        \n    return np.stack(feats)\n    \npre_train = np.concatenate([preproc('test', test_adc_info, \"FGS1\", 30*12), preproc('test', test_adc_info, \"AIRS-CH0\", 30)], axis=2)","metadata":{"execution":{"iopub.status.busy":"2024-09-21T02:39:17.73986Z","iopub.execute_input":"2024-09-21T02:39:17.740613Z","iopub.status.idle":"2024-09-21T02:39:31.559879Z","shell.execute_reply.started":"2024-09-21T02:39:17.74056Z","shell.execute_reply":"2024-09-21T02:39:31.55874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### fit polynoms for each sample","metadata":{}},{"cell_type":"code","source":"'''def phase_detector(signal):\n    phase1, phase2 = None, None\n    best_drop = 0\n    for i in range(50//2,150//2):        \n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase1 = i+(20+5)//2\n            best_drop = t1\n    \n    best_drop = 0\n    for i in range(200//2,250//2):\n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase2 = i-5//2\n            best_drop = t1\n    \n    return phase1, phase2\n\ndef try_s(signal, p1, p2, deg, s):\n    out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n    x, y = out, signal[out].tolist()\n    x = x + list(range(p1,p2))\n\n    y = y + (signal[p1:p2] * (1 + s[0])).tolist()\n    z = np.polyfit(x, y, deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n\n    if s < 1e-4:\n        return q + 1e3\n\n    return q\n    \ndef calibrate_signal(signal):\n    p1,p2 = phase_detector(signal)\n\n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 6):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n        \n        print(deg, q)\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n\n    return s, x, y, p(x)\n\ndef calibrate_train(signal):\n    p1,p2 = phase_detector(signal)\n    \n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 6):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n    \n    return s, p(np.arange(signal.shape[0])), p1, p2\n\n\ntrain = pre_train.copy()\nall_s = []\nfor i in range(len(test_adc_info)):\n    signal = train[i,:,1:].mean(axis=1)\n    s, p, p1, p2 = calibrate_train(pre_train[i,:,1:].mean(axis=1))\n    all_s.append(s)\n        \n#copy answer 283 times because we predict mean value\ntrain_s = np.repeat(np.array(all_s), 283).reshape((len(all_s), 283))        \ntrain_sigma = np.ones_like(train_s) * 0.00016'''\nclass SignalDataset(Dataset):\n    def __init__(self, data):\n        self.data = data\n        \n    def __len__(self):\n        return len(self.data)\n    \n    def __getitem__(self, idx):\n        return torch.tensor(self.data[idx], dtype=torch.float32)\n    \nclass CNNModel(nn.Module):\n    def __init__(self):\n        super(CNNModel, self).__init__()\n        self.conv1 = nn.Conv1d(in_channels=1, out_channels=16, kernel_size=3, padding=1)\n        self.conv2 = nn.Conv1d(in_channels=16, out_channels=32, kernel_size=3, padding=1)\n        self.pool = nn.MaxPool1d(2)\n        self.fc1 = nn.Linear(32 * 141, 64) # ？\n        self.fc2 = nn.Linear(64, 1)  # 结果维度？\n    \n    def forward(self, x):\n        x = self.pool(torch.relu(self.conv1(x)))\n        x = self.pool(torch.relu(self.conv2(x)))\n        x = x.view(-1, 32 * 141)  # 展平\n        x = torch.relu(self.fc1(x))\n        x = self.fc2(x)\n        return x\n    \ndef calculate_gll(y, mu_user, sigma_user):\n    sigma_user = torch.clamp(sigma_user, min=1e-6)\n    gll = -0.5 * (torch.log(2 * np.pi) + torch.log(sigma_user**2) + ((y - mu_user)**2) / (sigma_user**2))\n    return gll.sum()\n\n    \ndef train_model_with_gll(train_data, y_true, num_epochs=10, learning_rate=0.001):\n    model = CNNModel()\n    optimizer = optim.Adam(model.parameters(), lr=learning_rate)\n    train_loader = DataLoader(SignalDataset(train_data), batch_size=32, shuffle=True)\n\n    # 计算 L_ref 和 L_ideal\n    mean_y = np.mean(y_true)  # 真实值的均值\n    std_y = np.std(y_true)  # 真实值的标准差\n    L_ref = mean_y  # 使用均值作为参考\n    L_ideal = mean_y - 0.00001  # 理想情况下的不确定性，假设为均值减去一个小量\n    for epoch in range(num_epochs):\n        model.train()\n        running_loss = 0.0\n        for inputs in train_loader:\n            inputs = inputs.unsqueeze(1)  # 添加 channel 维度\n            optimizer.zero_grad()\n            mu_user = model(inputs)  # 模型输出\n            sigma_user = torch.tensor(np.ones_like(mu_user) * 0.00016, dtype=torch.float32)  # 不确定性\n            L = calculate_gll(y_true, mu_user, sigma_user)\n            loss = -L  \n            loss.backward()\n            optimizer.step()\n            running_loss += loss.item()\n        print(f'Epoch {epoch + 1}, Loss: {running_loss / len(train_loader)}')\n    return model\n\ndef predict(model, test_data):\n    test_loader = DataLoader(SignalDataset(test_data), batch_size=32, shuffle=False)\n    model.eval()\n    preds = []\n    with torch.no_grad():\n        for inputs in test_loader:\n            inputs = inputs.unsqueeze(1)  # 添加 channel 维度\n            outputs = model(inputs)\n            preds.append(outputs.numpy())\n    \n    return np.concatenate(preds, axis=0)\n\ntrain_data = preproc('train', train_adc_info, \"AIRS-CH0\", 30)\ntest_data = preproc('test', test_adc_info, \"AIRS-CH0\", 30)\nmodel = train_model(train_data, num_epochs=10)\npreds = predict(model, test_data)","metadata":{"execution":{"iopub.status.busy":"2024-09-21T02:39:31.56257Z","iopub.execute_input":"2024-09-21T02:39:31.562945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Probably we can accurately estimate sigma from train","metadata":{}},{"cell_type":"code","source":"n = 0\ns, x, y, y_new = calibrate_signal(pre_train[n, :, 1:].mean(axis=1))  # 确保 calibrate_signal 已定义\nplt.scatter(x, y, label='Original')\nplt.scatter(x, y_new, label='Calibrated')\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I call the orange line \"starline\". This is probably what we would see if the planet weren't in the way.","metadata":{}},{"cell_type":"markdown","source":"### Making submission","metadata":{}},{"cell_type":"code","source":"'''ss = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv')\n\npreds = train_s.clip(0)\nsigmas = train_sigma\nsubmission = pd.DataFrame(np.concatenate([preds,sigmas], axis=1), columns=ss.columns[1:])\nsubmission.index = test_adc_info.index\nsubmission.to_csv('submission.csv')'''","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''submission'''","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}