{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.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":101849,"databundleVersionId":13093295,"sourceType":"competition"},{"sourceId":9629432,"sourceType":"datasetVersion","datasetId":5846888},{"sourceId":12533930,"sourceType":"datasetVersion","datasetId":7912519},{"sourceId":12536938,"sourceType":"datasetVersion","datasetId":7914703}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"! pip install --no-index --find-links=/kaggle/input/ariel-2024-pqdm pqdm\n# docker v161 numpy-1.26.4 --> numpy-2.3.1 / scipy-1.15.3 --> scipy-1.16.0\n#! pip install --upgrade --no-index --find-links=/kaggle/input/whl-numpy-docker-v161 numpy==2.3.1\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:15:43.193813Z","iopub.execute_input":"2025-07-22T12:15:43.193941Z","iopub.status.idle":"2025-07-22T12:15:48.230767Z","shell.execute_reply.started":"2025-07-22T12:15:43.193926Z","shell.execute_reply":"2025-07-22T12:15:48.229582Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import dill\nimport pickle\nimport numpy as np; print(f\"numpy: {np.__version__}\")\nimport pandas as pd\nfrom pqdm.processes import pqdm\nimport itertools\nfrom astropy.stats import sigma_clip\nimport scipy; print(f\"scipy: {scipy.__version__}\")\nfrom scipy.optimize import minimize\nfrom scipy.ndimage import gaussian_filter\nfrom scipy.cluster.hierarchy import linkage, dendrogram, fcluster\nfrom scipy.spatial.distance import squareform\n  \n#dataset = \"train\"\ndataset = \"test\"\ntrain_calibration = True  # False\ntest_single_calibration = False\npath_train = '/kaggle/input/ariel-2025-06-01'\npath_input = '/kaggle/input/ariel-data-challenge-2025'\npath_dataset = f\"{path_input}/{dataset}\"\nadc_info = pd.read_csv(f\"{path_input}/adc_info.csv\")\nFGS1_gain, FGS1_offset = adc_info.iloc[0]['FGS1_adc_gain'], adc_info.iloc[0]['FGS1_adc_offset']\nAIRS_CH0_gain, AIRS_CH0_offset = adc_info.iloc[0]['AIRS-CH0_adc_gain'], adc_info.iloc[0]['AIRS-CH0_adc_offset']\ndel adc_info\nprint(f\"AIRS_CH0_gain: {AIRS_CH0_gain}, AIRS_CH0_offset: {AIRS_CH0_offset}, FGS1_adc_gain: {FGS1_gain}, FGS1_adc_offset: {FGS1_offset}\")\nstar_info = pd.read_csv(f\"{path_input}/{dataset}_star_info.csv\")\ntrain_labels = pd.read_csv(f\"{path_input}/train.csv\")\nplanet_df = star_info[['planet_id']].copy()\ntrain_labels.set_index(\"planet_id\", inplace=True)\nplanet_df.set_index(\"planet_id\", inplace=True)\nplanet_ids = planet_df.index.astype(int)\nprint(f\"star_info:{star_info.shape}, train_labels:{train_labels.shape}, planet_ids:{len(planet_ids)}\")\n\ndef save_as_dill(file_name, data, file_ext='.dill'):\n    with open(f\"{file_name}{file_ext}\", \"wb\") as file_handle:\n        dill.dump(data, file_handle, protocol=4)\n\ndef load_from_dill(file_name, file_ext='.dill'):\n    data = None\n    with open(f\"{file_name}{file_ext}\", \"rb\") as file_handle:\n        data = dill.load(file_handle)\n    return data\n\ndef save_as_pickle(file_name, data):\n    with open(f'{file_name}.pkl', 'wb') as file_handle:\n        pickle.dump(data, file_handle)\n\ndef load_from_pickle(file_name):\n    data = None\n    with open(f\"{file_name}.pkl\", \"rb\") as file_handle:\n        data = dill.load(file_handle)\n    return data\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:15:48.232901Z","iopub.execute_input":"2025-07-22T12:15:48.233249Z","iopub.status.idle":"2025-07-22T12:15:49.892696Z","shell.execute_reply.started":"2025-07-22T12:15:48.233204Z","shell.execute_reply":"2025-07-22T12:15:49.891702Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Calibrator:\n    cut_inf = 39\n    cut_sup = 321\n    sensor_to_sizes_dict = {\n        \"AIRS-CH0\": [[11250, 32, 356], [1, 32, cut_sup - cut_inf]],\n        \"FGS1\": [[135000, 32, 32], [1, 32, 32]],\n    }\n    sensor_to_linear_corr_dict = {\"AIRS-CH0\": (6, 32, 356), \"FGS1\": (6, 32, 32)}\n\n    def __init__(self, dataset, planet_id, sensor):\n        self.dataset = dataset\n        self.planet_id = planet_id\n        self.sensor = sensor\n\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(range(clean_signal.shape[1]), range(clean_signal.shape[2])):\n            poli = np.poly1d(linear_corr[:, x, y])\n            clean_signal[:, x, y] = poli(clean_signal[:, x, y])  # The process hangs. Need another numpy.\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 get_calibrated_signal(self, gain, offset, entry=0):\n        path_sensor = f\"{path_input}/{self.dataset}/{self.planet_id}/{self.sensor}\"\n        signal = pd.read_parquet(f\"{path_sensor}_signal_{entry}.parquet\").to_numpy()\n        dark_frame = pd.read_parquet(f\"{path_sensor}_calibration_{entry}/dark.parquet\", engine=\"pyarrow\").to_numpy()\n        dead_frame = pd.read_parquet(f\"{path_sensor}_calibration_{entry}/dead.parquet\", engine=\"pyarrow\").to_numpy()\n        flat_frame = pd.read_parquet(f\"{path_sensor}_calibration_{entry}/flat.parquet\", engine=\"pyarrow\").to_numpy()\n        linear_corr = (\n            pd.read_parquet(f\"{path_sensor}_calibration_{entry}/linear_corr.parquet\")\n            .values.astype(np.float64)\n            .reshape(self.sensor_to_linear_corr_dict[self.sensor])\n        )\n        signal = signal.reshape(self.sensor_to_sizes_dict[self.sensor][0])\n        signal = signal / gain + offset\n        hot = sigma_clip(dark_frame, sigma=5, maxiters=5).mask\n        if self.sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cut_inf : self.cut_sup]\n            dt = np.ones(len(signal)) * 0.1\n            dt[1::2] += 4.5  # @bilzard idea\n            linear_corr = linear_corr[:, :, self.cut_inf : self.cut_sup]\n            dark_frame = dark_frame[:, self.cut_inf : self.cut_sup]\n            dead_frame = dead_frame[:, self.cut_inf : self.cut_sup]\n            flat_frame = flat_frame[:, self.cut_inf : self.cut_sup]\n            hot = hot[:, self.cut_inf : self.cut_sup]\n        elif self.sensor == \"FGS1\":\n            dt = np.ones(len(signal)) * 0.1\n            dt[1::2] += 0.1\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        flat = flat_frame.reshape(self.sensor_to_sizes_dict[self.sensor][1])\n        flat[dead_frame.reshape(self.sensor_to_sizes_dict[self.sensor][1])] = np.nan\n        flat[hot.reshape(self.sensor_to_sizes_dict[self.sensor][1])] = np.nan\n        signal = signal / flat\n        return signal\n\ndef preprocess_signal(dataset, planet_id, sensor, gain, offset):\n    sensor_to_binning = {\"AIRS-CH0\": 30, \"FGS1\": 30 * 12}\n    sensor_to_binned_dict = {\n        \"AIRS-CH0\": [11250 // sensor_to_binning[\"AIRS-CH0\"] // 2, 282],\n        \"FGS1\": [135000 // sensor_to_binning[\"FGS1\"] // 2],\n    }\n    binning = sensor_to_binning[sensor]\n    signal = Calibrator(dataset=dataset, planet_id=planet_id, \n        sensor=sensor).get_calibrated_signal(gain=gain, offset=offset)\n    if sensor == \"AIRS-CH0\":\n        signal = signal[:, 10:22, :]\n    elif sensor == \"FGS1\":\n        signal = signal[:, 10:22, 10:22]\n        signal = signal.reshape(signal.shape[0], signal.shape[1] * signal.shape[2])\n    mean_signal = np.nanmean(signal, axis=1)\n    cds_signal = mean_signal[1::2] - mean_signal[0::2]\n    binned = np.zeros((sensor_to_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    if sensor == \"FGS1\":\n        binned = binned.reshape((binned.shape[0], 1))\n    return binned\n\ndef preprocessor(x):\n    return preprocess_signal(**x)\n\nif test_single_calibration:\n    args_fgs1 = dict(dataset=dataset, planet_id=planet_ids[0], sensor=\"FGS1\", gain=FGS1_gain, offset=FGS1_offset)\n    preprocessed_signal_fgs1 = preprocessor(args_fgs1)\n    print(preprocessed_signal_fgs1.shape)\n    args_airs_ch0 = dict(dataset=dataset, planet_id=planet_ids[0], sensor=\"AIRS-CH0\", gain=AIRS_CH0_gain, offset=AIRS_CH0_offset)\n    preprocessed_signal_airs_ch0 = preprocessor(args_airs_ch0)\n    print(preprocessed_signal_airs_ch0.shape)\n    display(preprocessed_signal_fgs1)\n    display(preprocessed_signal_airs_ch0)\n\nif train_calibration:\n    n_jobs = 4  # os.cpu_count()\n    print(f\"n_jobs: {n_jobs}\")\n    args_fgs1 = [dict(dataset=dataset, planet_id=planet_id, sensor=\"FGS1\", \n                    gain=FGS1_gain, offset=FGS1_offset) for planet_id in planet_ids]\n    preprocessed_signal_fgs1 = pqdm(args_fgs1, preprocessor, n_jobs=n_jobs)\n    preprocessed_signal_fgs1 = np.stack(preprocessed_signal_fgs1)\n    print(f\"preprocessed_signal_fgs1.shape: {preprocessed_signal_fgs1.shape}\")\n    np.save(f\"preprocessed_signal_fgs1.npy\", preprocessed_signal_fgs1, allow_pickle=False)\n\n    args_airs_ch0 = [dict(dataset=dataset, planet_id=planet_id, sensor=\"AIRS-CH0\", \n                        gain=AIRS_CH0_gain, offset=AIRS_CH0_offset) for planet_id in planet_ids]\n    preprocessed_signal_airs_ch0 = pqdm(args_airs_ch0, preprocessor, n_jobs=n_jobs)\n    preprocessed_signal_airs_ch0 = np.stack(preprocessed_signal_airs_ch0)\n    print(f\"preprocessed_signal_airs_ch0.shape:{preprocessed_signal_airs_ch0.shape}\")\n    np.save(f\"preprocessed_signal_airs_ch0.npy\", preprocessed_signal_airs_ch0, allow_pickle=False)\n\n    preprocessed_signal = np.concatenate([preprocessed_signal_fgs1, preprocessed_signal_airs_ch0], axis=2)\n    np.save(f\"preprocessed_signal.npy\", preprocessed_signal, allow_pickle=False)\nelse:\n    print(np.load(f\"{path_train}/preprocessed_signal_fgs1.npy\", allow_pickle=False).shape)\n    print(np.load(f\"{path_train}/preprocessed_signal_airs_ch0.npy\", allow_pickle=False).shape)\n    preprocessed_signal = signal_list = np.load(f\"{path_train}/preprocessed_signal.npy\", allow_pickle=False)\nprint(f\"preprocessed_signal.shape:{preprocessed_signal.shape}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:15:49.896426Z","iopub.execute_input":"2025-07-22T12:15:49.896758Z","execution_failed":"2025-07-22T12:17:02.466Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def predict_spectra(signal, gauss_sigma=False):\n    def phase_detector(signal):\n        MIN = np.argmin(signal[30:140]) + 30\n        signal1 = signal[:MIN]\n        signal2 = signal[MIN:]\n        first_derivative1 = np.gradient(signal1)\n        first_derivative1 /= first_derivative1.max()\n        first_derivative2 = np.gradient(signal2)\n        first_derivative2 /= first_derivative2.max()\n        phase1 = np.argmin(first_derivative1)\n        phase2 = np.argmax(first_derivative2) + MIN\n        return phase1, phase2\n    \n    def objective_to_minimize(s, original_signal, gauss_signal, phase1, phase2, gauss_phase_1, gauss_phase_2):\n        try:\n            delta = 2\n            power = 3\n            x = list(range(signal.shape[0] - delta * 4))\n            if gauss_phase_1 < delta:\n                if dataset == \"train\":\n                    print(f\"gauss_phase_1 < delta: {gauss_phase_1}, {delta}\")\n                gauss_phase_1 = delta\n            y = (\n                gauss_signal[: gauss_phase_1 - delta].tolist()\n                + (gauss_signal[gauss_phase_1 + delta : gauss_phase_2 - delta] * (1 + s)).tolist()\n                + gauss_signal[gauss_phase_2 + delta :].tolist()\n            )\n            if len(x) != len(y):               \n                # Bad First : s=[0.0001], x=179, y=366, signal.shape=(187, 283),                  gauss_phase_1=0,  gauss_phase_2=175\n                # Good First: s=[0.0001], x=179, y=179, signal.shape=(187, 31), gauss_signal=187, gauss_phase_1=65, gauss_phase_2=136\n                if dataset == \"train\":\n                    print(f\"First: s={s}, x={len(x)}, y={len(y)}, signal.shape={signal.shape}, gauss_signal={len(gauss_signal)}, gauss_phase_1={gauss_phase_1}, gauss_phase_2={gauss_phase_2}\")\n                min_len = min(len(x), len(y))\n                if min_len != 179:\n                    raise Exception()\n                x = x[:min_len]\n                y = y[:min_len]\n            z = np.polyfit(x, y, deg=power)\n            p = np.poly1d(z)\n            q = np.abs(p(x) - y).mean()\n            return q\n        except:\n            pass\n        try:\n            delta = 2\n            power = 3\n            x = list(range(signal.shape[0] - delta * 4))\n            if phase1 < delta:\n                if dataset == \"train\":\n                    print(f\"phase1 < delta: {phase1}, {delta}\")\n                phase1 = delta\n            y = (\n                original_signal[: phase1 - delta].tolist()\n                + (original_signal[phase1 + delta : phase2 - delta] * (1 + s)).tolist()\n                + original_signal[phase2 + delta :].tolist()\n            )\n            if len(x) != len(y):\n                # Bad Second : s=[0.0001], x=179, y=366, signal.shape=(187, 283),                      gauss_phase_1=0, gauss_phase_2=175\n                # Good second: s=[0.0001], x=179, y=179, signal.shape=(187, 31),  original_signal=187, gauss_phase_1=65, gauss_phase_2=136\n                if dataset == \"train\":\n                    print(f\"Second: s={s}, x={len(x)}, y={len(y)}, signal.shape={signal.shape}, original_signal={len(original_signal)}, gauss_phase_1={gauss_phase_1}, gauss_phase_2={gauss_phase_2}\")\n                min_len = min(len(x), len(y))\n                if min_len != 179:\n                    raise Exception()\n                x = x[:min_len]\n                y = y[:min_len]                    \n            z = np.polyfit(x, y, deg=power)\n            p = np.poly1d(z)\n            q = np.abs(p(x) - y).mean()\n            return q\n        except:\n            return 1.0\n    \n    def gaussian_denoise(signal, sigma=1):\n        return gaussian_filter(signal, sigma=sigma)\n\n    mean_signal = signal[:, :].mean(axis=1)\n    phase1, phase2 = phase_detector(mean_signal)\n    gauss_signal = gaussian_denoise(mean_signal, sigma=gauss_sigma)\n    gauss_phase_1, gauss_phase_2 = phase_detector(gauss_signal)\n    s = minimize(\n            fun=lambda s: objective_to_minimize(s, mean_signal, gauss_signal, phase1, phase2, gauss_phase_1, gauss_phase_2), \n            x0=[0.0001], \n            method=\"Nelder-Mead\"\n        ).x[0]\n    return s\n\ndef get_clusters(train_labels, threshold=0.002):\n    train_labels = train_labels.values[:,1:]\n    #train_labels = np.concatenate([train_labels[:,0:1], train_labels[:, 1:][:, ::-1]], axis=1)\n    correlation_matrix = np.corrcoef(train_labels, rowvar=False)\n    distance_matrix = 1 - correlation_matrix\n    condensed_distance_matrix = squareform(distance_matrix, checks=False)\n    Z = linkage(condensed_distance_matrix, method='complete')\n    clusters = fcluster(Z, threshold, criterion='distance')\n    return clusters\n\ndef _get_cluster_indizes(cluster_number, clusters):\n    columns_in_cluster = np.where(clusters == cluster_number)[0]\n    return columns_in_cluster\n\nclusters = get_clusters(train_labels)\nunique_clusters = np.unique(clusters)\ncluster_dict = {}\nfor cluster_n in unique_clusters:\n    cluster_dict[int(cluster_n)] = len(_get_cluster_indizes(cluster_n, clusters))\nprint(f\"Clusters:{cluster_dict}\")\ncluster_placeholder = np.zeros((planet_df.shape[0], 283))\nfor i in range(len(planet_ids)):\n    base_signal = preprocessed_signal[i,:,:]\n    for cluster_n in unique_clusters:\n        relevant_indizes = _get_cluster_indizes(cluster_n, clusters)\n        signal = base_signal[:,relevant_indizes]\n        for gauss_sigm in [2.0, 2.5, 3.0]:\n            try:\n                spectra = predict_spectra(signal, gauss_sigma=2.0)\n                cluster_placeholder[i, relevant_indizes] = spectra\n                break\n            except:\n                continue\n\ncluster_sigmas = np.ones_like(cluster_placeholder) * cluster_placeholder.clip(0).std(axis=1).reshape(cluster_placeholder.shape[0],-1)\n","metadata":{"trusted":true,"execution":{"execution_failed":"2025-07-22T12:17:02.466Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"predictions_spectra = np.array([predict_spectra(preprocessed_signal[i]) for i in range(len(preprocessed_signal))])\n","metadata":{"trusted":true,"execution":{"execution_failed":"2025-07-22T12:17:02.467Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\npredictions_spectra = np.repeat(np.array(predictions_spectra), 283).reshape((len(predictions_spectra), 283))\npredictions_spectra = predictions_spectra.clip(0)\n\nsubmission = pd.DataFrame(np.concatenate([predictions_spectra, cluster_sigmas], axis=1), columns=sample_submission.columns)\nsubmission.index = planet_ids  # sample_submission.index\n\nsubmission.to_csv(\"submission.csv\")\nsubmission\n","metadata":{"trusted":true,"execution":{"execution_failed":"2025-07-22T12:17:02.467Z"}},"outputs":[],"execution_count":null}]}