{"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":"gpu","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30786,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nfrom tqdm.notebook import tqdm\nimport numpy as np\nimport pandas as pd\nimport random\nimport scipy.signal as ss\nimport scipy\nimport joblib\nimport pickle\nimport numba\nfrom scipy.optimize import minimize\nfrom functools import partial\n\nfrom scipy import signal\nimport seaborn as sns\nimport itertools\n# from astropy.stats import sigma_clip\nimport time","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!mkdir /kaggle/avg_public","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def norm(v):\n    return v / v.sum()\n\ndef conv_rect_freq(x, r=4):\n    y = x.copy()\n    for i in range(r, x.shape[1] - r):\n        y[:, i] = x[:, i-r:i+r+1].mean(axis=1)\n    return y\n\ndef conv_hann_freq(x, r=4):\n    w = norm(ss.windows.hann(2*r+1, sym=True))\n    y = x.copy()\n    for i in range(r, x.shape[1] - r):\n        y[:, i] = (x[:, i-r:i+r+1] * w).sum(axis=1)\n    return y","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def format_mu(planet_codes, mu, wavelengths):\n    return pd.DataFrame(mu.clip(0, None), index=planet_codes, columns=wavelengths.columns)\n\ndef format_sigma(planet_codes, sigma):\n    return pd.DataFrame(sigma, index=planet_codes, columns=[f\"sigma_{i}\" for i in range(1, 284)])\n\ndef format_to_submission(planet_codes, mus, sigmas):\n    mus = pd.concat(mus)\n    sigmas = pd.concat(sigmas)\n    df_index = pd.DataFrame({'planet_id': list(planet_codes)}, index=planet_codes)\n    df_submission = df_index.join(mus).join(sigmas)\n    return df_submission.reset_index(drop=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SignalCalibratorPublic:\n    def apply_linear_corr(self, 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\n    def clean_dark(self, 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    \n    def calibrate(self, planet_id):\n        cut_inf, cut_sup = 0, 356\n\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/{planet_id}/AIRS-CH0_signal.parquet').to_numpy()\n        dark_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/' + str(planet_id) + '/' + 'AIRS-CH0' + '_calibration/dark.parquet').to_numpy()\n        dead_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/' + str(planet_id) + '/' + 'AIRS-CH0' + '_calibration/dead.parquet').to_numpy()\n        flat_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/' + str(planet_id) + '/' + 'AIRS-CH0' + '_calibration/flat.parquet').to_numpy()\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/test/' + str(planet_id) + '/' + 'AIRS-CH0' + '_calibration/linear_corr.parquet').values.astype(np.float64).reshape((6, 32, 356))\n        signal = signal.reshape([11250, 32, 356]) \n        gain = adc_info[adc_info['planet_id'] == planet_id][f'AIRS-CH0_adc_gain'].iloc[0]\n        offset = adc_info[adc_info['planet_id'] == planet_id][f'AIRS-CH0_adc_offset'].iloc[0]\n        signal = signal / gain + offset\n        \n        # hot = sigma_clip(\n        #     dark_frame, sigma=5, maxiters=5\n        # ).mask\n        \n        signal = signal[:, :, cut_inf:cut_sup] \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            \n        signal = signal.clip(0) #@graySnow idea\n        linear_corr_signal = self.apply_linear_corr(linear_corr, signal)\n        signal = self.clean_dark(linear_corr_signal, dark_frame, dt)\n        \n        flat = flat_frame.reshape([1, 32, cut_sup-cut_inf])\n        flat[dead_frame.reshape([1, 32, cut_sup-cut_inf])] = np.nan\n        # flat[hot.reshape([1, 32, cut_sup-cut_inf])] = np.nan\n        signal = signal / flat\n\n        signal = signal[:, 10:22 , :]\n        mean_signal = np.nanmean(signal, axis=1)\n        cds_signal = (mean_signal[1::2] - mean_signal[0::2])\n        return cds_signal\n\n    def calibrate_parallel(self, adc_info):\n        planet_ids = adc_info['planet_id']\n        feats = []\n        \n        feats = joblib.Parallel(n_jobs=-1, backend=\"loky\")(\n            joblib.delayed(self.calibrate)(planet_id) for planet_id in tqdm(list(planet_ids))\n        )\n            \n        return np.stack(feats)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SignalProcessor:\n    def __init__(self, transit_search_radius=500, transit_search_diff=150, drop_search_radius=300):\n        transit_start = 1820\n        transit_end = 3800\n        self.transit_search_radius = transit_search_radius\n        self.transit_start_area = [\n            transit_start - transit_search_radius,\n            transit_start + transit_search_radius,\n        ]\n        self.transit_end_area = [\n            transit_end - transit_search_radius,\n            transit_end + transit_search_radius,\n        ]\n        self.transit_search_diff = transit_search_diff\n        self.drop_search_radius = drop_search_radius\n\n        self.low_filter = ss.butter(4, 0.01, btype=\"low\", output=\"sos\")\n        self.high_filter = ss.butter(4, 0.01, btype=\"high\", output=\"sos\")\n\n    def denoise(self, x, axis=0, lpf=None):\n        if lpf is None:\n            lpf = self.low_filter\n        length = x.shape[0]\n        x_direct = ss.sosfiltfilt(lpf, x, axis=axis)\n        x_reversed = ss.sosfiltfilt(lpf, x[::-1], axis=axis)[::-1]\n        x_denoised = np.concatenate([x_reversed[:length // 2], x_direct[length // 2:]], axis=0)\n        return x_denoised\n    \n    def denoise_gauss(self, x, sigma, radius, axis=0):\n        x_denoised = gaussian_filter1d(x, sigma=sigma, axis=axis, radius=radius, mode='mirror')\n        return x_denoised\n    \n    def highpass(self, x, axis=0):\n        x = ss.sosfiltfilt(self.high_filter, x, axis=axis)\n        return x\n\n    def _transit_center_search(self, x):        \n        x_denoised = self.denoise(x)\n        x_denoised[:100] = x[:100].mean(axis=0)\n        x = x_denoised\n        \n        diff_vals = x[0 : -self.transit_search_diff] - x[self.transit_search_diff :]\n        transit_limit_left = (\n            diff_vals[self.transit_start_area[0] : self.transit_start_area[1]].argmax()\n            + self.transit_start_area[0]\n            + self.transit_search_diff // 2\n        )\n        transit_limit_right = (\n            diff_vals[self.transit_end_area[0] : self.transit_end_area[1]].argmin()\n            + self.transit_end_area[0]\n            + self.transit_search_diff // 2\n        )\n        return transit_limit_left, transit_limit_right\n    \n    def transit_search(self, x, ):\n        if x.ndim > 1:\n            x = x.mean(axis=tuple(range(x.ndim)[1:]))\n        x_denoised = self.denoise(x) \n        x_denoised[:100] = x[:100].mean(axis=0)\n        x = x_denoised\n\n        left_center, right_center = self._transit_center_search(x)\n\n        dy2 = np.diff(np.diff(x, 1), 1)\n        # left drop\n        dy2_left = dy2[left_center - self.drop_search_radius : left_center + self.drop_search_radius]\n        left_high = dy2_left.argmin() + left_center - self.drop_search_radius\n        left_low = dy2_left.argmax() + left_center - self.drop_search_radius\n        # right drop\n        dy2_right = dy2[right_center - self.drop_search_radius : right_center + self.drop_search_radius]\n        right_low = dy2_right.argmax() + right_center - self.drop_search_radius\n        right_high = dy2_right.argmin() + right_center - self.drop_search_radius\n\n        return np.hstack([left_high, left_low, right_low, right_high])\n\n    @staticmethod\n    def split(x, split_idxs, return_idxs=False):\n        xs = []\n        idxs = []\n        split_idxs = [0] + list(split_idxs) + [len(x)]\n        for i in range(len(split_idxs) - 1):\n            xs.append(x[split_idxs[i] : split_idxs[i+1]])\n            idxs.append(np.array(range(split_idxs[i], split_idxs[i+1])))\n        if not return_idxs:\n            return xs\n        else:\n            return xs, idxs\n        \n    # ======================================== EDGE FEATURES ==========================================\n\n    def edge_features(self, *x_list, split_idxs=None, drop_avg_guard=20, drop_avg_window=90):\n        x = np.concatenate(x_list, axis=1)\n        x_left_out, _, x_transit, _, x_right_out = self.split(x, split_idxs)\n        \n        left_high = x_left_out[- drop_avg_guard - drop_avg_window : - drop_avg_guard].mean(axis=0)\n        left_low = x_transit[drop_avg_guard : drop_avg_guard + drop_avg_window].mean(axis=0)\n        \n        right_high = x_right_out[drop_avg_guard : drop_avg_guard + drop_avg_window].mean(axis=0)\n        right_low = x_transit[- drop_avg_guard - drop_avg_window : - drop_avg_guard].mean(axis=0)\n\n        rel_drop_left = (left_high - left_low) / left_high\n        rel_drop_right = (right_high - right_low) / right_high\n\n        rel_drop = (rel_drop_left + rel_drop_right) / 2\n\n        return np.hstack([rel_drop])\n    \n    # ======================================= POLYNOM FEATURES ==========================================\n    \n    def poly_rel_drop_loss(self, x_out, idxs_out, x_transit, idxs_transit, poly_deg, rel_drop):\n        x = np.hstack([idxs_out, idxs_transit])\n        y = np.hstack([x_out, x_transit * (1 + rel_drop[0])])\n        \n        poly_coeffs = np.polyfit(x, y, poly_deg)\n        polynom_func = np.poly1d(poly_coeffs)\n        poly_mae = np.abs(polynom_func(x) - y).mean()\n        return poly_mae\n\n    def poly_rel_drop_optimizer(self, x_out, idxs_out, x_transit, idxs_transit, max_poly_deg=4):\n        rel_drop_buffer, loss_value_buffer = np.zeros(max_poly_deg-1), np.zeros(max_poly_deg-1)\n        for poly_deg in range(2, max_poly_deg + 1):\n            f_loss = partial(self.poly_rel_drop_loss, x_out, idxs_out, x_transit, idxs_transit, poly_deg)\n            rel_drop =  minimize(f_loss, [0.001], method='Powell').x[0]\n            loss_value = f_loss([rel_drop])\n            \n            rel_drop_buffer[poly_deg - 2] = rel_drop\n            loss_value_buffer[poly_deg - 2] = loss_value\n        rel_drop = rel_drop_buffer[loss_value_buffer.argmin()]\n        return np.hstack([rel_drop])\n    \n    def poly_features_binned(self, *x_list, split_idxs=None, b=10, max_poly_deg=4):\n        x_array = np.concatenate(x_list, axis=1)\n        features = []\n        for i in range(x_array.shape[1]):\n            x = x_array[:, i]\n\n            dist = 75\n\n            xs, idxs = self.split(x, split_idxs, return_idxs=True)\n            x_left_out, _, x_transit, _, x_right_out = xs\n            idx_left_out, _, idx_transit, _, idx_right_out = idxs\n\n            x_out = np.hstack([x_left_out[:-dist], x_right_out[dist:]])\n            x_out = x_out[len(x_out) % b :]\n            x_out = x_out.reshape(-1, b).mean(axis=1)\n            x_transit = x_transit[dist:-dist]\n            x_transit = x_transit[len(x_transit) % b :]\n            x_transit = x_transit.reshape(-1, b).mean(axis=1)\n\n            idx_out = np.hstack([idx_left_out[:-dist], idx_right_out[dist:]])\n            idx_out = idx_out[len(idx_out) % b :]\n            idx_out = idx_out.reshape(-1, b).mean(axis=1)\n            idx_transit = idx_transit[dist:-dist]\n            idx_transit = idx_transit[len(idx_transit) % b :]\n            idx_transit = idx_transit.reshape(-1, b).mean(axis=1)\n\n            rel_drop = self.poly_rel_drop_optimizer(x_out, idx_out, x_transit, idx_transit, max_poly_deg)\n            features.append(rel_drop)\n        return np.hstack(features)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_folder = '/kaggle/input/ariel-data-challenge-2024/test'\nadc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal_calibrator = SignalCalibratorPublic()\nsignal_processor = SignalProcessor(transit_search_radius=1200,\n                                   transit_search_diff=250,\n                                   drop_search_radius=500)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data preprocessing","metadata":{}},{"cell_type":"code","source":"s = signal_calibrator.calibrate_parallel(adc_info)\nfor i, pc in enumerate(adc_info['planet_id']):\n    np.save(f'/kaggle/avg_public/{pc}_AIRS-CH0.npy', s[i])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_X(planet_code, processor, folder, cut=(39, 321), b=10, max_poly_deg=2):\n    cutl, cutr = cut\n\n    x_airs_base = np.load(f'/kaggle/{folder}/{planet_code}_AIRS-CH0.npy')\n    x_airs_conv = conv_hann_freq(x_airs_base, r=12) #  CONV\n    x_airs_base = x_airs_base[:, cutl:cutr]\n    x_airs_base = x_airs_base[:, ::-1]\n    split_idxs = processor.transit_search(x_airs_base.mean(axis=1, keepdims=True))\n    \n    x_airs_conv = x_airs_conv[:, cutl:cutr]\n    x_airs_conv = x_airs_conv[:, ::-1]\n\n    def spectral_division(x, chunk_size, chunk_num, b=10, max_poly_deg=2):\n        chunk_idxs = [(i*chunk_size, i*chunk_size + chunk_size) for i in range(chunk_num)]\n        h_edge, h_poly = [], []\n        for l, r in chunk_idxs:\n            x_chunk = x[:, l:r].mean(axis=1, keepdims=True)\n            X_edge = processor.edge_features(processor.denoise(x_chunk), split_idxs=split_idxs, drop_avg_guard=20, drop_avg_window=90).mean()\n            X_poly = processor.poly_features_binned(processor.denoise(x_chunk), split_idxs=split_idxs, b=b, max_poly_deg=max_poly_deg).mean()\n            h_edge.append(X_edge)\n            h_poly.append(X_poly)\n        h_edge, h_poly = np.hstack(h_edge), np.hstack(h_poly)\n        return h_edge, h_poly\n\n    x_airs = x_airs_base[:, 0:282].mean(axis=1, keepdims=True)\n    X_edge = processor.edge_features(x_airs, split_idxs=split_idxs).mean()\n    X_poly = processor.poly_features_binned(signal_processor.denoise(x_airs), split_idxs=split_idxs, b=1, max_poly_deg=4).mean()\n    h_1_e, h_1_p = np.hstack([X_edge]), np.hstack([X_poly])\n    h_2_e, h_2_p = spectral_division(x_airs_base, chunk_size=141, chunk_num=2, b=b, max_poly_deg=max_poly_deg)\n    h_4_e, h_4_p = spectral_division(x_airs_base, chunk_size=70, chunk_num=4, b=b, max_poly_deg=max_poly_deg)\n    h_47_e, h_47_p = spectral_division(x_airs_conv, chunk_size=6, chunk_num=47, b=b, max_poly_deg=max_poly_deg)\n\n    return h_1_e, h_1_p, h_2_e, h_2_p, h_4_e, h_4_p, h_47_e, h_47_p\n\nh_1_e, h_1_p, h_2_e, h_2_p, h_4_e, h_4_p, h_47_e, h_47_p = create_X(adc_info['planet_id'].iloc[0], signal_processor, 'avg_public', cut=(39, 321), b=10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = time.time()\n\n  ###################################\n ###          STAR 0             ###\n###################################\n\nXXs = []\ncutls = [48] #  [36, 38, 40, 42, 44, 46, 48, 50, 52, 54]\nfor l in cutls:\n    X = joblib.Parallel(n_jobs=-1, backend=\"loky\")(\n        joblib.delayed(create_X)(pc, signal_processor, 'avg_public', (l, l+282), 10, 2) for pc in  tqdm(adc_info['planet_id'][adc_info['star'] == 0])\n    )\n    XXs.append(list(zip(*X)))\n\nXX_0 = []\nfor n in range(len(XXs[0])):\n    XX_0.append(sum([np.hstack([XXs[i][n]]) for i in range(len(cutls))]) / len(cutls))\n    \n  ###################################\n ###          STAR 1             ###\n###################################\n    \nXXs = []\ncutls = [48]\nfor l in cutls:\n    X = joblib.Parallel(n_jobs=-1, backend=\"loky\")(\n        joblib.delayed(create_X)(pc, signal_processor, 'avg_public', (l, l+282), 10, 3) for pc in  tqdm(adc_info['planet_id'][adc_info['star'] == 1])\n    )\n    XXs.append(list(zip(*X)))\n\nXX_1 = []\nfor n in range(len(XXs[0])):\n    XX_1.append(sum([np.hstack([XXs[i][n]]) for i in range(len(cutls))]) / len(cutls))\n    \n  ###################################\n ###          STAR 2             ###\n###################################   \n    \nXXs = []\ncutls = [41]\nfor l in cutls:\n    X = joblib.Parallel(n_jobs=-1, backend=\"loky\")(\n        joblib.delayed(create_X)(pc, signal_processor, 'avg_public', (l, l+282), 10, 3) for pc in  tqdm(adc_info['planet_id'][adc_info['star'] == 2])\n    )\n    XXs.append(list(zip(*X)))\n\nXX_2 = []\nfor n in range(len(XXs[0])):\n    XX_2.append(sum([np.hstack([XXs[i][n]]) for i in range(len(cutls))]) / len(cutls))\n    \n  ###################################\n ###          STAR 3             ###\n###################################\n    \nXXs = []\ncutls = [46]\nfor l in cutls:\n    X = joblib.Parallel(n_jobs=-1, backend=\"loky\")(\n        joblib.delayed(create_X)(pc, signal_processor, 'avg_public', (l, l+282), 10, 3) for pc in  tqdm(adc_info['planet_id'][adc_info['star'] == 3])\n    )\n    XXs.append(list(zip(*X)))\n\nXX_3 = []\nfor n in range(len(XXs[0])):\n    XX_3.append(sum([np.hstack([XXs[i][n]]) for i in range(len(cutls))]) / len(cutls))\n\nprint('time:', time.time() - s)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class XModel:\n    def __init__(self, std_threshold=1):\n        self.std_threshold = std_threshold\n\n    def predict_mu(self, h_1_e, h_1_p, h_2_e, h_2_p, h_4_e, h_4_p, h_47_e, h_47_p):\n\n        h_1 = (h_1_e * 0.2 + h_1_p) / 1.2\n        h_2 = (h_2_e * 0.3 + h_2_p) / 1.3\n        h_4 = (h_4_e * 0.5 + h_4_p) / 1.5\n        h_47 = (h_47_e * 0.5 + h_47_p) / 1.5\n\n        X_const = (h_1 + 0.2 * h_2.mean(axis=1, keepdims=True)) / 1.2\n\n        X_diff = h_47\n\n        X_std = np.hstack([h_4[:, :3]])\n        \n        mus = []\n        for i in range(len(X_const)):\n            \n            mu_const = (X_const[i]).sum()\n            mu_const = np.tile(mu_const, 283)\n            \n            mu_diff = signal_processor.denoise(np.repeat(X_diff[i], 6), axis=0, lpf=ss.bessel(2, 0.03, btype=\"low\", output=\"sos\"))\n            mu_diff = np.hstack([mu_diff[0], mu_diff[:-24], np.full(24, mu_const[0])])\n\n            if X_std[i].std() < self.std_threshold:\n                mu = (mu_const + mu_diff * (X_std[i].std() / self.std_threshold)**0.6) / (1 + (X_std[i].std() / self.std_threshold)**0.6)\n            else:\n                mu = mu_diff\n            \n            mus.append(mu)\n        pred = np.stack(mus)\n        return pred\n    \n    def predict_sigma(self, h_1_e, h_1_p, h_2_e, h_2_p, h_4_e, h_4_p, h_47_e, h_47_p):\n\n        h_1 = (h_1_e * 0.2 + h_1_p) / 1.2\n        h_4 = (h_4_e * 0.5 + h_4_p) / 1.5\n\n        X_std_const = np.hstack([h_1, h_1, h_1,\n                                 h_4])\n\n        sigmas = []\n        for i in range(len(X_std_const)):\n            if X_std_const[i].std() < self.std_threshold:\n                sigma = X_std_const[i].std()\n            else:\n                sigma = X_std_const[i].std()\n            sigma = np.tile(sigma, 283)\n            sigmas.append(sigma)\n        pred = np.stack(sigmas)\n        return pred","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xmodel = XModel(std_threshold=100e-6)\n\n  ###################################\n ###             Mu              ###\n###################################\n\nmu_0 = xmodel.predict_mu(*XX_0)\nmu_0 = format_mu(adc_info['planet_id'][adc_info['star'] == 0], mu_0, wavelengths)\n\nmu_1 = xmodel.predict_mu(*XX_1)\nmu_1 = format_mu(adc_info['planet_id'][adc_info['star'] == 1], mu_1, wavelengths)\n\nmu_2 = xmodel.predict_mu(*XX_2)\nmu_2 = format_mu(adc_info['planet_id'][adc_info['star'] == 2], mu_2, wavelengths)\n\nmu_3 = xmodel.predict_mu(*XX_3)\nmu_3 = format_mu(adc_info['planet_id'][adc_info['star'] == 3], mu_3, wavelengths)\n\n  ###################################\n ###            Sigma            ###\n###################################\n\nsigma_0 = xmodel.predict_sigma(*XX_0)\nsigma_0 = np.clip(sigma_0, a_min=3e-5, a_max=np.inf)\nsigma_0 = format_sigma(adc_info['planet_id'][adc_info['star'] == 0], sigma_0)\n\nsigma_1 = xmodel.predict_sigma(*XX_1)\nsigma_1 = np.clip(sigma_1, a_min=3e-5, a_max=np.inf)\nsigma_1 = format_sigma(adc_info['planet_id'][adc_info['star'] == 1], sigma_1)\n\nsigma_2 = xmodel.predict_sigma(*XX_2)\nsigma_2 = np.clip(sigma_2, a_min=5e-5, a_max=np.inf)\nsigma_2 = format_sigma(adc_info['planet_id'][adc_info['star'] == 2], sigma_2)\n\nsigma_3 = xmodel.predict_sigma(*XX_3)\nsigma_3 = np.clip(sigma_3, a_min=5e-5, a_max=np.inf)\nsigma_3 = format_sigma(adc_info['planet_id'][adc_info['star'] == 3], sigma_3)\n\n\n\nsubmission_df = format_to_submission(adc_info['planet_id'], \n                                     mus=[mu_0, mu_1, mu_2, mu_3],\n                                     sigmas=[sigma_0, sigma_1, sigma_2, sigma_3])\nsubmission_df.to_csv('submission.csv', index=False)\n\npd.read_csv('submission.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}