{"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":"none","dataSources":[{"sourceType":"competition","sourceId":70367,"databundleVersionId":9188054},{"sourceType":"datasetVersion","sourceId":9629432,"datasetId":5846888,"databundleVersionId":9854288},{"sourceType":"datasetVersion","sourceId":9774241,"datasetId":5941595,"databundleVersionId":10016265}],"dockerImageVersionId":30761,"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","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:49:53.298883Z","iopub.execute_input":"2024-10-31T17:49:53.299546Z","iopub.status.idle":"2024-10-31T17:50:07.294371Z","shell.execute_reply.started":"2024-10-31T17:49:53.299504Z","shell.execute_reply":"2024-10-31T17:50:07.29283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Librairies","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport pandas.api.types\nimport scipy.stats\nimport cv2\n\nfrom tqdm import tqdm\nfrom pqdm.processes import pqdm\n\nimport itertools\n\nfrom scipy.optimize import minimize\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.linear_model import Ridge, LinearRegression\n\nimport plotly.express as px\n\nfrom astropy.stats import sigma_clip\n\nfrom scipy.optimize import curve_fit\nfrom scipy.ndimage import gaussian_filter1d","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-31T17:50:07.296704Z","iopub.execute_input":"2024-10-31T17:50:07.297074Z","iopub.status.idle":"2024-10-31T17:50:10.058076Z","shell.execute_reply.started":"2024-10-31T17:50:07.297035Z","shell.execute_reply":"2024-10-31T17:50:10.056963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Signal Preprocessing","metadata":{}},{"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(\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 get_calibrated_signal(self):\n        signal = pd.read_parquet(\n            f\"/kaggle/input/ariel-data-challenge-2024/{self.dataset}/{self.planet_id}/{self.sensor}_signal.parquet\"\n        ).to_numpy()\n        dark_frame = pd.read_parquet(\n            f\"/kaggle/input/ariel-data-challenge-2024/{self.dataset}/{self.planet_id}/{self.sensor}_calibration/dark.parquet\",\n            engine=\"pyarrow\",\n        ).to_numpy()\n        dead_frame = pd.read_parquet(\n            f\"/kaggle/input/ariel-data-challenge-2024/{self.dataset}/{self.planet_id}/{self.sensor}_calibration/dead.parquet\",\n            engine=\"pyarrow\",\n        ).to_numpy()\n        flat_frame = pd.read_parquet(\n            f\"/kaggle/input/ariel-data-challenge-2024/{self.dataset}/{self.planet_id}/{self.sensor}_calibration/flat.parquet\",\n            engine=\"pyarrow\",\n        ).to_numpy()\n        linear_corr = (\n            pd.read_parquet(\n                f\"/kaggle/input/ariel-data-challenge-2024/{self.dataset}/{self.planet_id}/{self.sensor}_calibration/linear_corr.parquet\"\n            )\n            .values.astype(np.float64)\n            .reshape(self.sensor_to_linear_corr_dict[self.sensor])\n        )\n\n        signal = signal.reshape(self.sensor_to_sizes_dict[self.sensor][0])\n        gain = adc_info.loc[self.planet_id, f\"{self.sensor}_adc_gain\"]\n        offset = adc_info.loc[self.planet_id, f\"{self.sensor}_adc_offset\"]\n        signal = signal / gain + offset\n\n        hot = sigma_clip(dark_frame, sigma=5, maxiters=5).mask\n\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\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(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\n\nclass Preprocessor:\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\n    def __init__(self, dataset, planet_id, sensor):\n        self.dataset = dataset\n        self.planet_id = planet_id\n        self.sensor = sensor\n        self.binning = self.sensor_to_binning[sensor]\n\n    def preprocess_signal(self):\n        signal = Calibrator(\n            dataset=self.dataset, planet_id=self.planet_id, sensor=self.sensor\n        ).get_calibrated_signal()\n\n        if self.sensor == \"AIRS-CH0\":\n            signal = signal[:, 10:22, :]\n        elif self.sensor == \"FGS1\":\n            signal = signal[:, 10:22, 10:22]\n            signal = signal.reshape(\n                signal.shape[0], signal.shape[1] * signal.shape[2]\n            )\n\n        mean_signal = np.nanmean(signal, axis=1)\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        binned = np.zeros((self.sensor_to_binned_dict[self.sensor]))\n        for j in range(cds_signal.shape[0] // self.binning):\n            binned[j] = cds_signal[\n                j * self.binning : j * self.binning + self.binning\n            ].mean(axis=0)\n\n        if self.sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        return binned\n\n\ndef preprocessor(x):\n    return Preprocessor(**x).preprocess_signal()\n","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:10.059807Z","iopub.execute_input":"2024-10-31T17:50:10.060365Z","iopub.status.idle":"2024-10-31T17:50:10.084347Z","shell.execute_reply.started":"2024-10-31T17:50:10.060317Z","shell.execute_reply":"2024-10-31T17:50:10.083186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = \"test\"\nadc_info = pd.read_csv(\n    \"/kaggle/input/ariel-data-challenge-2024/\" + f\"{dataset}_adc_info.csv\",\n    index_col=\"planet_id\",\n)\naxis_info = pd.read_parquet(\"/kaggle/input/ariel-data-challenge-2024/axis_info.parquet\")\nplanet_ids = adc_info.index","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:10.087772Z","iopub.execute_input":"2024-10-31T17:50:10.088368Z","iopub.status.idle":"2024-10-31T17:50:10.280417Z","shell.execute_reply.started":"2024-10-31T17:50:10.088301Z","shell.execute_reply":"2024-10-31T17:50:10.279359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"args_fgs1 = [\n    dict(dataset=dataset, planet_id=planet_id, sensor=\"FGS1\")\n    for planet_id in planet_ids\n]\npreprocessed_signal_fgs1 = pqdm(args_fgs1, preprocessor, n_jobs=4)\n\nargs_airs_ch0 = [\n    dict(dataset=dataset, planet_id=planet_id, sensor=\"AIRS-CH0\")\n    for planet_id in planet_ids\n]\npreprocessed_signal_airs_ch0 = pqdm(args_airs_ch0, preprocessor, n_jobs=4)\n\npreprocessed_signal = np.concatenate(\n    [np.stack(preprocessed_signal_fgs1), np.stack(preprocessed_signal_airs_ch0)], axis=2\n)\npreprocessed_signal.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:10.282017Z","iopub.execute_input":"2024-10-31T17:50:10.282563Z","iopub.status.idle":"2024-10-31T17:50:21.79137Z","shell.execute_reply.started":"2024-10-31T17:50:10.282509Z","shell.execute_reply":"2024-10-31T17:50:21.790239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modelization","metadata":{}},{"cell_type":"code","source":"def phase_detector(signal):\n    MIN = np.argmin(signal[30:140]) + 30\n    signal1 = signal[:MIN]\n    signal2 = signal[MIN:]\n\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\n    phase1 = np.argmin(first_derivative1)\n    phase2 = np.argmax(first_derivative2) + MIN\n\n    return phase1, phase2","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.793091Z","iopub.execute_input":"2024-10-31T17:50:21.793468Z","iopub.status.idle":"2024-10-31T17:50:21.800293Z","shell.execute_reply.started":"2024-10-31T17:50:21.793432Z","shell.execute_reply":"2024-10-31T17:50:21.799144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def predict_spectra(signal):\n    def objective_to_minimize(s):\n        delta = 2\n        power = 3\n        \n        x = list(range(signal.shape[0] - delta * 4))\n        #x = list(range(signal.shape[0]))\n        #x = x[: phase1 - delta] + x[phase1 + delta : phase2 - delta] + x[phase2 + delta :]\n        \n        y = (\n            signal[: phase1 - delta].tolist()\n            + (signal[phase1 + delta : phase2 - delta] * (1 + s)).tolist()\n            + signal[phase2 + delta :].tolist()\n        )\n\n        z = np.polyfit(x, y, deg=power)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        # q = ((p(x) - y) ** 2).mean()\n        return q, z  # zを返す\n\n    signal_for_phase_detector = signal.mean(axis=1)\n    phase1, phase2 = phase_detector(signal_for_phase_detector)\n    \n    scale = signal.mean()\n    signal = signal / signal.mean(0, keepdims=True) * scale\n    signal = signal.mean(1)\n    # zも一緒に受け取るようにする\n    result = minimize(fun=lambda s: objective_to_minimize(s)[0], x0=[0.0001], method=\"Nelder-Mead\")\n    s = result.x[0]\n\n    # objective_to_minimizeからzも得る\n    q, z = objective_to_minimize(s)\n\n    return s, z, phase1, phase2","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.801627Z","iopub.execute_input":"2024-10-31T17:50:21.80196Z","iopub.status.idle":"2024-10-31T17:50:21.816613Z","shell.execute_reply.started":"2024-10-31T17:50:21.801918Z","shell.execute_reply":"2024-10-31T17:50:21.815254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def predict_spectra_with_cache(signal, phase1, phase2, x0):\n    def objective_to_minimize(s):\n        delta = 2\n        power = 3\n        \n        x = list(range(signal.shape[0] - delta * 4))\n        #x = list(range(signal.shape[0]))\n        #x = x[: phase1 - delta] + x[phase1 + delta : phase2 - delta] + x[phase2 + delta :]\n        \n        y = (\n            signal[: phase1 - delta].tolist()\n            + (signal[phase1 + delta : phase2 - delta] * (1 + s)).tolist()\n            + signal[phase2 + delta :].tolist()\n        )\n\n        z = np.polyfit(x, y, deg=power)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        # q = ((p(x) - y) ** 2).mean()\n        return q, z  # zを返す\n\n    result = minimize(fun=lambda s: objective_to_minimize(s)[0], x0=[x0], method=\"Nelder-Mead\")\n    s = result.x[0]\n\n    # objective_to_minimizeからzも得る\n    q, z = objective_to_minimize(s)\n\n    return s, z, phase1, phase2","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.818668Z","iopub.execute_input":"2024-10-31T17:50:21.81905Z","iopub.status.idle":"2024-10-31T17:50:21.830993Z","shell.execute_reply.started":"2024-10-31T17:50:21.819014Z","shell.execute_reply":"2024-10-31T17:50:21.830003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def poly3d(xy, b00, b01, b10, b11, b20, b21, b30, b31):\n    x, y = xy\n    return (b00 + b01*x) + (b10 + b11*x)*y + (b20 + b21*x)*y**2 + (b30 + b31*x)*y**3\n\n# フィッティングする関数（任意のサイズの画像に対応）\ndef fit_image(Z_image):\n    # 画像の高さと幅を取得\n    height, width = Z_image.shape\n    Z_image = Z_image.copy()\n    Z_image = cv2.GaussianBlur(Z_image, (1, 5), 0)\n    \n    # グリッド作成\n    x = np.linspace(0, 1, width)\n    y = np.linspace(0, 1, height)\n    X, Y = np.meshgrid(x, y)\n    \n    # 画像データを1次元化\n    x_flat = X.flatten()\n    y_flat = Y.flatten()\n    z_flat = Z_image.flatten()\n    \n    # フィッティング\n    popt, pcov = curve_fit(poly3d, (x_flat, y_flat), z_flat)\n    \n    # フィッティング結果をグリッドに適用して新しいz値を計算\n    Z_fitted = poly3d((X, Y), *popt)\n\n    return Z_fitted, popt, pcov","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.832482Z","iopub.execute_input":"2024-10-31T17:50:21.833363Z","iopub.status.idle":"2024-10-31T17:50:21.84598Z","shell.execute_reply.started":"2024-10-31T17:50:21.833292Z","shell.execute_reply":"2024-10-31T17:50:21.844898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_each_freq(img):\n    s, z, phase1, phase2 = predict_spectra(img)\n    ori_s = s\n    \n    img = np.concatenate([img[:, 0:1], np.flip(img[:, 1:], axis=1)], axis=1)\n    \n    img = img / img.mean(0, keepdims=True) * img.mean()\n    img = gaussian_filter1d(img, 15)\n    \n    result = []\n    \n    for i in range(img.shape[1]):\n        s, z, phase1, phase2 = predict_spectra_with_cache(img[:, i], phase1, phase2, ori_s)\n        result.append(s)\n        \n    result = np.stack(result)\n    return result","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.849139Z","iopub.execute_input":"2024-10-31T17:50:21.849857Z","iopub.status.idle":"2024-10-31T17:50:21.858065Z","shell.execute_reply.started":"2024-10-31T17:50:21.849814Z","shell.execute_reply":"2024-10-31T17:50:21.857014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_seq_stat(signal):\n    return np.array([\n        signal.mean(),\n        signal.std(),\n        signal.max(),\n        signal.min(),\n        signal[0],\n        signal[-1],\n        signal[:len(signal) // 2].mean(),\n        signal[len(signal) // 2:].mean(),\n        signal[:len(signal) // 2].std(),\n        signal[len(signal) // 2:].std(),\n        signal[:len(signal) // 4].mean(),\n        signal[len(signal) // 4:len(signal) // 2].mean(),\n        signal[len(signal) // 2:len(signal) // 4 * 3].mean(),\n        signal[len(signal) // 4 * 3:].mean(),\n        signal[:len(signal) // 4].std(),\n        signal[len(signal) // 4:len(signal) // 2].std(),\n        signal[len(signal) // 2:len(signal) // 4 * 3].std(),\n        signal[len(signal) // 4 * 3:].std(),\n    ])","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.859505Z","iopub.execute_input":"2024-10-31T17:50:21.859901Z","iopub.status.idle":"2024-10-31T17:50:21.873237Z","shell.execute_reply.started":"2024-10-31T17:50:21.859859Z","shell.execute_reply":"2024-10-31T17:50:21.871912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission = pd.read_csv(\n    \"/kaggle/input/ariel-data-challenge-2024/sample_submission.csv\",\n    index_col=\"planet_id\",\n)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.874786Z","iopub.execute_input":"2024-10-31T17:50:21.875875Z","iopub.status.idle":"2024-10-31T17:50:21.903045Z","shell.execute_reply.started":"2024-10-31T17:50:21.875831Z","shell.execute_reply":"2024-10-31T17:50:21.901902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv\")\ndf = df.set_index(\"planet_id\")\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:21.904671Z","iopub.execute_input":"2024-10-31T17:50:21.90506Z","iopub.status.idle":"2024-10-31T17:50:21.926599Z","shell.execute_reply.started":"2024-10-31T17:50:21.905021Z","shell.execute_reply":"2024-10-31T17:50:21.925448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = np.load(\"/kaggle/input/ariel2024-stats/sigmas_and_stats4.npz\")\ntrain_stats = d[\"stats\"]\ntrain_sigmas = d[\"sigmas\"]\nsig_model_2 = RandomForestRegressor()\nsig_model_2.fit(train_stats, train_sigmas)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:25.157348Z","iopub.execute_input":"2024-10-31T17:50:25.157777Z","iopub.status.idle":"2024-10-31T17:50:29.648461Z","shell.execute_reply.started":"2024-10-31T17:50:25.157737Z","shell.execute_reply":"2024-10-31T17:50:29.647354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = np.load(\"/kaggle/input/ariel2024-stats/sigmas_and_stats3.npz\")\ntrain_stats = d[\"stats\"]\ntrain_sigmas = d[\"sigmas\"]\nsig_model = RandomForestRegressor()\nsig_model.fit(train_stats, train_sigmas)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:33.344799Z","iopub.execute_input":"2024-10-31T17:50:33.345207Z","iopub.status.idle":"2024-10-31T17:50:38.706294Z","shell.execute_reply.started":"2024-10-31T17:50:33.345168Z","shell.execute_reply":"2024-10-31T17:50:38.705351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = np.load(\"/kaggle/input/ariel2024-stats/model1_feat.npz\")\nmodel1_feat = d[\"model1_feat\"]\ntargets = d[\"targets\"]\nmodel1 = LinearRegression()\nmodel1.fit(model1_feat, targets)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:38.708199Z","iopub.execute_input":"2024-10-31T17:50:38.708544Z","iopub.status.idle":"2024-10-31T17:50:38.830948Z","shell.execute_reply.started":"2024-10-31T17:50:38.70851Z","shell.execute_reply":"2024-10-31T17:50:38.829937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = np.load(\"/kaggle/input/ariel2024-stats/model2_feat.npz\")\nmodel2_feat = d[\"model2_feat\"]\ntargets = d[\"targets\"]\nmodel2 = LinearRegression()\nmodel2.fit(model2_feat, targets)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:50:38.832025Z","iopub.execute_input":"2024-10-31T17:50:38.832433Z","iopub.status.idle":"2024-10-31T17:50:38.966653Z","shell.execute_reply.started":"2024-10-31T17:50:38.832384Z","shell.execute_reply":"2024-10-31T17:50:38.965166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_spectra = []\nsigmas = []\n\nfor planet_id, img in zip(sample_submission.index, preprocessed_signal):\n    star_id = df.loc[planet_id, \"star\"]\n\n    result2 = fit_each_freq(img)\n    \n    s, z, phase1, phase2 = predict_spectra(img)\n    ori_s = s\n    \n    img = np.concatenate([img[:, 0:1], np.flip(img[:, 1:], axis=1)], axis=1)\n\n    delta = 2\n    new_img = np.concatenate([\n        img[:phase1 - delta],\n        img[phase1 + delta : phase2 - delta] * (1 + s),\n        img[phase2 + delta :]\n    ])\n    \n    scale = new_img.mean()\n    new_img = new_img / new_img.mean(0, keepdims=True) * scale\n\n    fitted, popt, pcov = fit_image(new_img)\n    \n    a = img[:phase1 - delta]\n    b = img[phase1 + delta:phase2 - delta]\n    c = img[phase2 + delta:]\n    d = np.concatenate([a, b, c])\n    scale = d.mean()\n    d = d / d.mean(0, keepdims=True) * scale\n    \n    e = d / fitted\n    f = np.concatenate([e[:phase1 - delta], e[phase2 - delta * 3:]]).mean(0)\n    g = e[phase1 - delta:phase2 - delta * 3].mean(0)\n    h = (f - g) / f\n    #s = np.ones_like(h) * s\n    \n    h = h * s.mean() / h.mean()\n    ori_h = h.copy()\n    \n    h = gaussian_filter1d(h, 10)\n    h[len(h) // 2:] = gaussian_filter1d(h[len(h) // 2:], 24)\n    \n    result1 = h\n    \n    if star_id > 1:\n        result_ave = 0.5 * result1 + 0.5 * result2\n        ensemble_results = result_ave.copy()\n        ensemble_results[:30] = result2[:30]\n        ensemble_results[190:] = result1[190:]\n        ensemble_results[230:] = s.mean()\n        ensemble_results[180-3:187+3] = ori_h[180-3:187+3].mean(keepdims=True)\n        predictions_spectra.append(ensemble_results)\n\n        h_stat = get_seq_stat(h)\n        stat = np.concatenate([popt.flatten(), pcov.flatten(), h_stat])\n        \n        si_pred = sig_model_2.predict([stat])[0]\n    \n    else:\n        x1 = np.concatenate([result1[40:150], result1[180:187]])\n        pred1 = model1.predict([x1])[0]\n        x2 = np.concatenate([result2[:40], result1[180:187]])\n        pred2 = model2.predict([x2])[0]\n        pred2[67:139] = pred1[67:139]\n\n        predictions_spectra.append(pred2)\n\n\n        h_stat = get_seq_stat(h)\n        #stat = np.concatenate([popt.flatten(), pcov.flatten(), h_stat])\n        #stats.append(stat)\n        s1 = get_seq_stat(pred2)\n        #s2 = get_seq_stat(result1[40:150])\n        #s3 = get_seq_stat(result2[:40])\n        stat = np.concatenate([popt.flatten(), pcov.flatten(), s1, h_stat])    \n\n        si_pred = sig_model.predict([stat])[0]\n\n        \n    if star_id > 1:\n        si_pred = si_pred * 1.2\n    else:\n        si_pred = si_pred * 1.2\n    \n    sigmas.append([si_pred])\n    \npredictions_spectra = np.array(predictions_spectra)\nsigmas = np.ones_like(predictions_spectra) * np.array(sigmas)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:54:53.454191Z","iopub.execute_input":"2024-10-31T17:54:53.454846Z","iopub.status.idle":"2024-10-31T17:54:54.519834Z","shell.execute_reply.started":"2024-10-31T17:54:53.454791Z","shell.execute_reply":"2024-10-31T17:54:54.518626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"predictions_spectra = predictions_spectra.clip(0)\n\nsubmission = pd.DataFrame(\n    np.concatenate([predictions_spectra, sigmas], axis=1),\n    columns=sample_submission.columns,\n)\nsubmission.index = sample_submission.index","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:54:54.521836Z","iopub.execute_input":"2024-10-31T17:54:54.52235Z","iopub.status.idle":"2024-10-31T17:54:54.528845Z","shell.execute_reply.started":"2024-10-31T17:54:54.522274Z","shell.execute_reply":"2024-10-31T17:54:54.527732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.to_csv(\"submission.csv\")\nsubmission","metadata":{"execution":{"iopub.status.busy":"2024-10-31T17:54:54.556277Z","iopub.execute_input":"2024-10-31T17:54:54.557272Z","iopub.status.idle":"2024-10-31T17:54:54.585505Z","shell.execute_reply.started":"2024-10-31T17:54:54.557228Z","shell.execute_reply":"2024-10-31T17:54:54.584355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}