{"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":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":204366594,"sourceType":"kernelVersion"}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom astropy.stats import sigma_clip\nfrom scipy.signal import savgol_filter\nfrom scipy.fft import fft, ifft\nimport pywt\nfrom sklearn.tree import DecisionTreeRegressor\nimport pickle","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-31T10:08:42.181306Z","iopub.execute_input":"2024-10-31T10:08:42.181807Z","iopub.status.idle":"2024-10-31T10:08:43.607538Z","shell.execute_reply.started":"2024-10-31T10:08:42.181754Z","shell.execute_reply":"2024-10-31T10:08:43.606147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = '/kaggle/input/ariel-data-challenge-2024/'\ntest_df = pd.read_csv(path + 'test_adc_info.csv')\nplanets = test_df['planet_id']\naxis_info = pd.read_parquet(path + 'axis_info.parquet')\nrmse = pd.read_csv('/kaggle/input/ariel-training/rmse.csv').values\nslope = pd.read_csv('/kaggle/input/ariel-training/slope.csv').values[0]\nstars = test_df['star'].values","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:08:43.609293Z","iopub.execute_input":"2024-10-31T10:08:43.609821Z","iopub.status.idle":"2024-10-31T10:08:43.665539Z","shell.execute_reply.started":"2024-10-31T10:08:43.609777Z","shell.execute_reply":"2024-10-31T10:08:43.664314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_dt = np.ones(135000) * 0.1\nf_dt[1::2] += 0.1\nf_dt = f_dt[1::2]\na_dt = axis_info['AIRS-CH0-integration_time'].dropna().values\na_dt[1::2] += 0.1\na_dt = a_dt[1::2]\n\ndef double_sampling(signal):\n    #signal_low = signal[0::2, :, :]\n    signal_high = signal[1::2, :, :]\n    #signal = signal_high - signal_low\n    return signal_high\n\ndef adc_convertion(signal, gain, offset):\n    signal = signal / gain + offset\n    return signal\n\ndef masking(signal, dead, dark):\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal)\n    return signal\n\ndef clipping(signal):\n    signal = np.ma.clip(signal, 0, None)\n    return signal\n\ndef poly_correction(signal, poly):\n    poly = np.tile(poly, (signal.shape[0], 1, 1, 1))\n    signal = poly[:, 0] + signal * poly[:, 1] + signal ** 2 * poly[:, 2] + signal ** 3 * poly[:, 3] + signal ** 4 * poly[:, 4]\\\n        + signal ** 5 * poly[:, 5]\n    return signal\n\ndef dark_subtraction(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 flat_correction(signal, flat):\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    signal = signal / flat\n    return signal\n\ndef fft_filter(signal, sampling_rate, cutoff_freq):\n    for i in range(signal.shape[1]):\n        signal_1d = signal[:, i]\n        freqs = np.fft.fftfreq(len(signal_1d), d=1/sampling_rate)\n        signal_fft = fft(signal_1d)\n        filter_mask = np.abs(freqs) <= cutoff_freq\n        filtered_fft = signal_fft * filter_mask\n        filtered_signal = ifft(filtered_fft).real\n        signal[:, i] = filtered_signal\n    return signal","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:08:43.667164Z","iopub.execute_input":"2024-10-31T10:08:43.667533Z","iopub.status.idle":"2024-10-31T10:08:43.68721Z","shell.execute_reply.started":"2024-10-31T10:08:43.667496Z","shell.execute_reply":"2024-10-31T10:08:43.685955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gaussian_curve(x, mean, sigma, multiplier, constant):\n    f_x = np.exp(-((x - mean) ** 2) / (2 * sigma ** 2)) * multiplier + constant\n    return f_x\n\ndef smoothen(x, m, v, k):\n    return m * np.exp(-v * x) + k\n\ndef phase_detection_f(f_signal, window, sgwindow, shift, whiskers):\n    \n    mov_avg = f_signal[0].rolling(window=window).mean()\n    savgol = savgol_filter(mov_avg, window_length=sgwindow, polyorder=1)\n    df = pd.DataFrame(savgol).rename(columns={0: 'mean'}).join(pd.DataFrame(savgol).rename(columns={0: 'mean'}).shift(shift), rsuffix='2')\n    df['diff'] = df['mean'] - df['mean2']\n    argmin, argmax = df.iloc[10000:33750]['diff'].argmin() + 10000, df.iloc[33750:57500]['diff'].argmax() + 33750\n    df['diff2'] = np.where(((df.index > argmin - shift) & (df.index < argmin + shift)) \\\n        | ((df.index > argmax - shift) & (df.index < argmax + shift)), np.nan, df['diff'])\n    df1, df2 = df.iloc[argmin-shift-whiskers:argmin+shift+whiskers], df.iloc[argmax-shift-whiskers:argmax+shift+whiskers]\n    z = np.full(len(df), np.nan)\n    a1, b = np.polyfit(df1['diff2'].dropna().index, df1['diff2'].dropna(), 1)\n    z[argmin-shift-whiskers:argmin+shift+whiskers] = a1 * df1.index + b\n    a2, b = np.polyfit(df2['diff2'].dropna().index, df2['diff2'].dropna(), 1)\n    z[argmax-shift-whiskers:argmax+shift+whiskers] = a2 * df2.index + b\n    df['diff'] = df['diff'] - 0.6 * np.nan_to_num(z, nan=0)\n    argmin, argmax = df.iloc[argmin-shift:argmin+shift]['diff'].argmin() + argmin-shift, df.iloc[argmax-shift:argmax+shift]['diff'].argmax() + argmax-shift\n    s_min, s_max = df.loc[argmin, 'diff'], df.loc[argmax, 'diff'] #df.loc[argmin-25:argmin+25, 'diff'].mean(), df.loc[argmax-25:argmax+25, 'diff'].mean()\n    signal_decrease, signal_increase = abs(s_min) / df.loc[argmin, 'mean2'], abs(s_max) / df.loc[argmax, 'mean']\n    confidence_ratio = abs(a2/a1 * z[argmax]/z[argmin]) ** 0.5\n    if confidence_ratio > 1:\n        smooth = smoothen(confidence_ratio - 1, m=-0.15, v=0.9, k=0.65)\n        signal_delta_f = smooth * signal_decrease + (1 - smooth) * signal_increase\n        z_ratio = smooth * abs(z[argmin]/s_min) + (1 - smooth) * abs(z[argmax]/s_max)\n        z_slope = abs(a1 * smooth) + abs(a2 * (1 - smooth))\n    else:\n        smooth = smoothen(confidence_ratio ** -1 - 1, m=-0.15, v=0.9, k=0.65)\n        signal_delta_f = (1 - smooth) * signal_decrease + smooth * signal_increase\n        z_ratio = (1 - smooth) * abs(z[argmin]/s_min) + smooth * abs(z[argmax]/s_max)\n        z_slope = abs(a1 * (1 - smooth)) + abs(a2 * smooth)\n    signal_ratio = max(signal_increase, signal_decrease) / min(signal_increase, signal_decrease)\n    noise = f_signal[0].std()/f_signal[0].mean() * 10000\n    m_noise = f_signal[0].rolling(window=window).std().median()/f_signal[0].mean() * 10000\n    ma_noise = mov_avg.rolling(window=window*2).std().median()/mov_avg.mean() * 10000\n    savgol_d = savgol_filter(mov_avg, window_length=sgwindow, polyorder=1, deriv=1)\n    savgol_d_shift = np.squeeze(pd.DataFrame(savgol_d).shift(1).values)\n    savgol_oscillations = np.sum(np.where((savgol_d < 0) & (savgol_d_shift > 0) | (savgol_d > 0) & (savgol_d_shift < 0), 1, 0))\n    \n    return signal_delta_f, signal_ratio, noise, m_noise, ma_noise, z_ratio, z_slope, savgol_oscillations\n\ndef phase_detection_a(a_signal, normalizer, window, sgwindow, shift, whiskers):\n    \n    a_signal = a_signal / a_signal.mean() ** normalizer\n    a_signal[0] = a_signal.iloc[:, 34:316].mean(axis=1)\n    a_signal[0] = a_signal[0] / a_signal[0].mean()\n    mov_avg = a_signal[0].rolling(window=window).mean()\n    savgol = savgol_filter(mov_avg, window_length=sgwindow, polyorder=1) # Savitzky–Golay filter\n    df = pd.DataFrame(savgol).rename(columns={0: 'mean'}).join(pd.DataFrame(savgol).rename(columns={0: 'mean'}).shift(shift), rsuffix='2')\n    df['diff'] = df['mean'] - df['mean2']\n    argmin, argmax = df.iloc[625:2812]['diff'].argmin() + 625, df.iloc[2812:5000]['diff'].argmax() + 2812\n    df['diff2'] = np.where(((df.index > argmin - shift) & (df.index < argmin + shift)) \\\n        | ((df.index > argmax - shift) & (df.index < argmax + shift)), np.nan, df['diff'])\n    df1, df2 = df.iloc[argmin-shift-1000:argmin+shift+whiskers], df.iloc[argmax-shift-whiskers:argmax+shift+1000]\n    z = np.full(len(df), np.nan)\n    a1, b = np.polyfit(df1['diff2'].dropna().index, df1['diff2'].dropna(), 1)\n    z[argmin-shift-1000:argmin+shift+whiskers] = a1 * df1.index + b\n    a2, b = np.polyfit(df2['diff2'].dropna().index, df2['diff2'].dropna(), 1)\n    z[argmax-shift-whiskers:argmax+shift+1000] = a2 * df2.index + b\n    df['diff'] = df['diff'] - 0.93 * np.nan_to_num(z, nan=0)\n    argmin, argmax = df.iloc[argmin-shift:argmin+shift]['diff'].argmin() + argmin-shift, df.iloc[argmax-shift:argmax+shift]['diff'].argmax() + argmax-shift\n    s_min, s_max = df.loc[argmin, 'diff'], df.loc[argmax, 'diff']\n    signal_decrease, signal_increase = abs(s_min) / df.loc[argmin, 'mean2'], abs(s_max) / df.loc[argmax, 'mean']\n    confidence_ratio = abs(a2/a1 * z[argmax]/z[argmin]) ** 0.5\n    if confidence_ratio > 1:\n        smooth = smoothen(confidence_ratio - 1, m=-0.20, v=0.08, k=0.70)\n    else:\n        smooth = 1 - smoothen(confidence_ratio ** -1 - 1, m=-0.20, v=0.08, k=0.70)\n    signal_delta_a = smooth * signal_decrease + (1 - smooth) * signal_increase\n        \n    err = np.ones(5)\n    for deg in range(1, 6):\n        df['adjusted'] = np.where(((df.index > argmin - shift) & (df.index < argmin)) | ((df.index > argmax - shift) & (df.index < argmax)), np.nan, np.where( \\\n            (df.index >= argmin) & (df.index <= argmax - shift), df['mean'] * (1 + signal_delta_a), df['mean']))\n        coefs = np.polyfit(df['adjusted'].dropna().index, df['adjusted'].dropna(), deg)\n        df['poly'] = np.polyval(coefs, df.index)\n        err[deg-1] = ((df['adjusted'] - df['poly']) ** 2).sum() * 1.02 ** (deg - 1)\n    best_deg = err.argmin() + 1\n\n    df['adjusted'] = np.where(((df.index > argmin - shift) & (df.index < argmin)) | ((df.index > argmax - shift) & (df.index < argmax)), np.nan, np.where( \\\n        (df.index >= argmin) & (df.index <= argmax - shift), df['mean'] * (1 + signal_delta_a), df['mean']))\n    coefs = np.polyfit(df['adjusted'].dropna().index, df['adjusted'].dropna(), best_deg)\n    df['poly'] = np.polyval(coefs, df.index)\n    df['mean'] = df['mean'] - df['poly'] + df['poly'].mean()\n    \n    signal_delta_a = 1 - df.iloc[argmin:argmax-shift]['mean'].mean()\n    signal_delta_a += smooth * (df.iloc[0:argmin-shift]['mean'].mean() - 1) + (1 - smooth) * (df.iloc[argmax:5625]['mean'].mean() - 1)\n    \n    return signal_delta_a, argmin, argmax, df, df['poly'], best_deg\n\ndef phase_detection_a_sensor(a_signal, j, borders, sigma, multiplier, argmin, argmax, window, sgwindow, shift, whiskers):\n    \n    w = np.round(gaussian_curve(np.arange(282+2*borders), mean=j+borders, sigma=sigma, multiplier=1, constant=0))\n    w_mean = (a_signal.iloc[:, 34-borders:316+borders] * w).mean(axis=1)\n    m_noise = w_mean.rolling(window=window).std().median()/w_mean.mean() * 10000\n    w = np.round(gaussian_curve(np.arange(282+2*borders), mean=j+borders, sigma=sigma, multiplier=multiplier/(0.1+m_noise*1.8), constant=1) + gaussian_curve(np.arange(282+2*borders), j+borders, 33, 3, 0), 2)\n    w_mean = (a_signal.iloc[:, 34-borders:316+borders] * w).mean(axis=1)\n    w_mean = w_mean / w_mean.mean()\n    mov_avg = w_mean.rolling(window=window).mean()\n    savgol = savgol_filter(mov_avg, window_length=sgwindow, polyorder=1) - poly[i] + poly[i].mean()\n    df = pd.DataFrame(savgol).rename(columns={0: 'mean'}).join(pd.DataFrame(savgol).rename(columns={0: 'mean'}).shift(shift), rsuffix='2')\n    df['diff'] = df['mean'] - df['mean2']\n    df['diff2'] = np.where(((df.index > argmin - shift) & (df.index < argmin + shift)) | ((df.index > argmax - shift) & (df.index < argmax + shift)), np.nan, df['diff'])\n    df1, df2 = df.copy().iloc[argmin-shift-1000:argmin+shift+whiskers], df.copy().iloc[argmax-shift-whiskers:argmax+shift+1000]\n    z = np.full(len(df), np.nan)\n    a1, b = np.polyfit(df1['diff2'].dropna().index, df1['diff2'].dropna(), 1)\n    z[argmin-shift-1000:argmin+shift+whiskers] = a1 * df1.index + b\n    a2, b = np.polyfit(df2['diff2'].dropna().index, df2['diff2'].dropna(), 1)\n    z[argmax-shift-whiskers:argmax+shift+1000] = a2 * df2.index + b\n    df['diff'] = df['diff'] - 0.93 * np.nan_to_num(z, nan=0)\n    s_min, s_max = df.loc[argmin, 'diff'], df.loc[argmax, 'diff']\n    signal_decrease, signal_increase = abs(s_min) / df.loc[argmin, 'mean2'], abs(s_max) / df.loc[argmax, 'mean']\n    confidence_ratio = abs(a2/a1 * z[argmax]/z[argmin]) ** 0.5\n    if confidence_ratio > 1:\n        smooth = smoothen(confidence_ratio - 1, m=-0.20, v=0.08, k=0.70)\n    else:\n        smooth = 1 - smoothen(confidence_ratio ** -1 - 1, m=-0.20, v=0.08, k=0.70)\n    z_ratio = smooth * abs(z[argmin]/s_min) + (1 - smooth) * abs(z[argmax]/s_max)\n    z_slope = abs(a1 * smooth) + abs(a2 * (1 - smooth))\n    signal_delta = smooth * signal_decrease + (1 - smooth) * signal_increase\n    \n    signal_ratio = max(signal_increase, signal_decrease) / min(signal_increase, signal_decrease)\n    noise = a_signal.iloc[:, j+34].std()/a_signal.iloc[:, j].mean() * 10000\n    ma_noise = mov_avg.rolling(window=window*2).std().median()/mov_avg.mean() * 10000\n    savgol_d = savgol_filter(mov_avg, window_length=sgwindow, polyorder=1, deriv=1)\n    savgol_d_shift = np.squeeze(pd.DataFrame(savgol_d).shift(1).values)\n    savgol_oscillations = np.sum(np.where((savgol_d < 0) & (savgol_d_shift > 0) | (savgol_d > 0) & (savgol_d_shift < 0), 1, 0))\n    \n    err = np.ones(3)\n    for deg in range(1, 4):\n        df['adjusted'] = np.where(((df.index > argmin - shift) & (df.index < argmin)) | ((df.index > argmax - shift) & (df.index < argmax)), np.nan, np.where( \\\n            (df.index >= argmin) & (df.index <= argmax - shift), df['mean'] * (1 + signal_delta), df['mean']))\n        coefs = np.polyfit(df['adjusted'].dropna().index, df['adjusted'].dropna(), deg)\n        df['poly'] = np.polyval(coefs, df.index)\n        err[deg-1] = ((df['adjusted'] - df['poly']) ** 2).sum() * 1.02 ** (deg - 1)\n    best_deg = err.argmin() + 1\n    df['adjusted'] = np.where(((df.index > argmin - shift) & (df.index < argmin)) | ((df.index > argmax - shift) & (df.index < argmax)), np.nan, np.where( \\\n        (df.index >= argmin) & (df.index <= argmax - shift), df['mean'] * (1 + signal_delta), df['mean']))\n    coefs = np.polyfit(df['adjusted'].dropna().index, df['adjusted'].dropna(), best_deg)\n    df['poly'] = np.polyval(coefs, df.index)\n    df['mean'] = df['mean'] - df['poly'] + df['poly'].mean()\n    \n    signal_delta = 1 - df.iloc[argmin:argmax-shift]['mean'].mean()\n    signal_delta += smooth * (df.iloc[0:argmin-shift]['mean'].mean() - 1) + (1 - smooth) * (df.iloc[argmax:5625]['mean'].mean() - 1)\n    \n    return df, savgol, signal_delta, signal_ratio, noise, m_noise, ma_noise, z_ratio, z_slope, savgol_oscillations","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:08:43.69129Z","iopub.execute_input":"2024-10-31T10:08:43.691692Z","iopub.status.idle":"2024-10-31T10:08:43.751329Z","shell.execute_reply.started":"2024-10-31T10:08:43.691652Z","shell.execute_reply":"2024-10-31T10:08:43.750085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CustomRandomForestRegressor:\n    def __init__(self, n_estimators=10, max_depth=None, max_features=None, min_samples_leaf=1, min_samples_split=2, random_state=None):\n        self.n_estimators = n_estimators\n        self.max_depth = max_depth\n        self.max_features = max_features\n        self.min_samples_leaf = min_samples_leaf\n        self.min_samples_split = min_samples_split\n        self.random_state = random_state\n        self.trees = []\n        \n    def fit(self, x, y):\n        self.trees = []\n        for i in range(self.n_estimators):\n            np.random.seed(self.random_state + i)\n            indices = np.random.choice(x.shape[0], np.random.choice(x.shape[0]), replace=True) # Bootstrap sampling\n            x_sample = x[indices]\n            y_sample = y[indices]\n            tree = DecisionTreeRegressor(max_depth=self.max_depth, max_features=self.max_features, min_samples_leaf=self.min_samples_leaf, \n                                         min_samples_split=self.min_samples_split, random_state=self.random_state)\n            tree.fit(x_sample, y_sample)\n            self.trees.append(tree)\n    \n    def predict(self, x):\n        tree_predictions = np.array([tree.predict(x) for tree in self.trees])\n        predictions = np.mean(tree_predictions ** 2, axis=0) ** 0.5\n        \n        return predictions","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:08:43.752917Z","iopub.execute_input":"2024-10-31T10:08:43.753423Z","iopub.status.idle":"2024-10-31T10:08:43.76877Z","shell.execute_reply.started":"2024-10-31T10:08:43.753373Z","shell.execute_reply":"2024-10-31T10:08:43.767578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"signal_delta_f, signal_delta_a = np.zeros(len(planets)), np.zeros(len(planets))\nargmin, argmax, shift, whiskers = np.zeros(len(planets)).astype(int), np.zeros(len(planets)).astype(int), np.zeros(len(planets)).astype(int), np.zeros(len(planets)).astype(int)\nsignal_ratio, noise, m_noise, ma_noise = np.zeros((len(planets), 283)), np.zeros((len(planets), 283)), np.zeros((len(planets), 283)), np.zeros((len(planets), 283))\nsavgol_oscillations, z_ratio, z_slope = np.zeros((len(planets), 283)), np.zeros((len(planets), 283)), np.zeros((len(planets), 283))\nx, x_meanratio = np.zeros((len(planets), 283)), np.zeros((len(planets), 283))\npoly = np.zeros((len(planets), 5625))\ndeg = np.zeros(len(planets))\n\nfor i in range(0, len(planets)):\n    \n    f_signal = pd.read_parquet(path + 'test/' + str(planets[i]) + '/FGS1_signal.parquet').values.astype(np.float32).reshape(135000, 32, 32)[:, 10:22, 10:22]\n    a_signal = pd.read_parquet(path + 'test/' + str(planets[i]) + '/AIRS-CH0_signal.parquet').values.astype(np.float32).reshape(11250, 32, 356)[:, 10:22, :]\n    f_dead = pd.read_parquet(path + 'test/' + str(planets[i]) + '/FGS1_calibration/dead.parquet').values.astype(np.float32).reshape(32, 32)[10:22, 10:22]\n    f_dark = pd.read_parquet(path + 'test/' + str(planets[i]) + '/FGS1_calibration/dark.parquet').values.astype(np.float32).reshape(32, 32)[10:22, 10:22]\n    f_flat = pd.read_parquet(path + 'test/' + str(planets[i]) + '/FGS1_calibration/flat.parquet').values.astype(np.float32).reshape(32, 32)[10:22, 10:22]\n    a_dead = pd.read_parquet(path + 'test/' + str(planets[i]) + '/AIRS-CH0_calibration/dead.parquet').values.astype(np.float32).reshape(32, 356)[10:22, :]\n    a_dark = pd.read_parquet(path + 'test/' + str(planets[i]) + '/AIRS-CH0_calibration/dark.parquet').values.astype(np.float32).reshape(32, 356)[10:22, :]\n    a_flat = pd.read_parquet(path + 'test/' + str(planets[i]) + '/AIRS-CH0_calibration/flat.parquet').values.astype(np.float32).reshape(32, 356)[10:22, :]\n    f_offset, f_gain = test_df.iloc[i]['FGS1_adc_offset'], test_df.iloc[i]['FGS1_adc_gain']\n    a_offset, a_gain = test_df.iloc[i]['AIRS-CH0_adc_offset'], test_df.iloc[i]['AIRS-CH0_adc_gain']\n    f_poly = pd.read_parquet(path + 'test/' + str(planets[i]) + '/FGS1_calibration/linear_corr.parquet').values.astype(np.float32).reshape(6, 32, 32)[:, 10:22, 10:22]\n    a_poly = pd.read_parquet(path + 'test/' + str(planets[i]) + '/AIRS-CH0_calibration/linear_corr.parquet').values.astype(np.float32).reshape(6, 32, 356)[:, 10:22, :]\n    \n    f_signal = double_sampling(f_signal)\n    a_signal = double_sampling(a_signal)\n    f_signal = adc_convertion(f_signal, f_gain, f_offset)\n    a_signal = adc_convertion(a_signal, a_gain, a_offset)\n    f_signal = masking(f_signal, f_dead, f_dark)\n    a_signal = masking(a_signal, a_dead, a_dark)\n    f_signal = clipping(f_signal)\n    a_signal = clipping(a_signal)\n    f_signal = poly_correction(f_signal, f_poly)\n    a_signal = poly_correction(a_signal, a_poly)\n    f_signal = dark_subtraction(f_signal, f_dark, f_dt)\n    a_signal = dark_subtraction(a_signal, a_dark, a_dt)\n    f_signal = clipping(f_signal)\n    a_signal = clipping(a_signal)\n    f_signal = flat_correction(f_signal, f_flat)\n    a_signal = flat_correction(a_signal, a_flat)\n    f_signal = f_signal.reshape(67500, 144).mean(axis=1).reshape(67500, 1)\n    a_signal = a_signal.mean(axis=1)\n    \n    f_signal = fft_filter(f_signal, sampling_rate=1/0.2, cutoff_freq=0.02)\n    a_signal = fft_filter(a_signal, sampling_rate=1/4.6, cutoff_freq=0.006)\n    f_signal = pd.DataFrame(f_signal)\n    a_signal = pd.DataFrame(a_signal).iloc[:, ::-1]\n    \n    signal_delta_a[i], argmin[i], argmax[i], df, _, _ = phase_detection_a(a_signal, normalizer=0.7, window=387, sgwindow=39, shift=572, whiskers=400)\n    shift[i] = 375 + np.clip(int(round(signal_delta_a[i] * 41092, 0)), 375, 600)\n    whiskers[i] = np.clip(argmax[i] - argmin[i] - 2 * shift[i], 2, 1000)\n    signal_delta_a[i], argmin[i], argmax[i], df, poly[i], deg[i] = phase_detection_a(a_signal, normalizer=0.7, window=357, sgwindow=47, shift=shift[i], whiskers=whiskers[i])\n    whiskers[i] = np.clip(argmax[i] - argmin[i] - 2 * shift[i], 2, 1000)\n    shift_f = 2280 + int(round(signal_delta_a[i] * 174000, 0))\n    signal_delta_f[i], signal_ratio[i, 0], noise[i, 0], m_noise[i, 0], ma_noise[i, 0], z_ratio[i, 0], z_slope[i, 0], savgol_oscillations[i, 0] = phase_detection_f(f_signal, \n        window=1920, sgwindow=850, shift=shift_f, whiskers=2080)\n\n    signal_delta = np.zeros(282)\n    a_signal_normalized = a_signal / a_signal.mean() ** 0.3\n    _, _, signal_delta[0], _, _, _, _, _, _, _ = phase_detection_a_sensor(a_signal_normalized, \n            j=0, borders=34, sigma=9.5, multiplier=110, argmin=argmin[i], argmax=argmax[i], window=357, sgwindow=47, shift=shift[i], whiskers=whiskers[i])\n    x_meanratio[i, 0] = abs(signal_delta[0] - signal_delta_a[0]) / signal_delta_a[0]\n    for j in range(282):\n        multiplier = 110*(1+x_meanratio[i, j]*2)**10/signal_ratio[i, j]**8/(1+z_ratio[i, j])**6\n        df, _, signal_delta[j], signal_ratio[i, j+1], noise[i, j+1], m_noise[i, j+1], ma_noise[i, j+1], z_ratio[i, j+1], z_slope[i, j+1], savgol_oscillations[i, j+1] = phase_detection_a_sensor(a_signal_normalized, \n            j, borders=34, sigma=9.5, multiplier=multiplier, argmin=argmin[i], argmax=argmax[i], window=357, sgwindow=47, shift=shift[i], whiskers=whiskers[i])\n        x_meanratio[i, j+1] = abs(signal_delta[j] - signal_delta_a[i]) / signal_delta_a[i]\n    \n    x[i, 1:] = signal_delta.copy()\n    \n    print(f'Progress: {i+1}/{len(planets)}', end='\\r')\n    \nx[:, 0] = signal_delta_f","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:08:43.770764Z","iopub.execute_input":"2024-10-31T10:08:43.771167Z","iopub.status.idle":"2024-10-31T10:09:04.667365Z","shell.execute_reply.started":"2024-10-31T10:08:43.771125Z","shell.execute_reply":"2024-10-31T10:09:04.666154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_sigma = pickle.load(open('/kaggle/input/ariel-training/model_sigma.pkl', 'rb'))\n\ny_pred = x * slope","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:09:04.668837Z","iopub.execute_input":"2024-10-31T10:09:04.669222Z","iopub.status.idle":"2024-10-31T10:09:09.454845Z","shell.execute_reply.started":"2024-10-31T10:09:04.669174Z","shell.execute_reply":"2024-10-31T10:09:09.453072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_ratio = (x - np.squeeze(np.tile(x.mean(axis=1), (283, 1))).T) / np.squeeze(np.tile(x.mean(axis=1), (283, 1))).T ** 0.4\nsensor = np.tile(np.arange(0,283), (len(planets), 1))\nx_std = np.squeeze(np.tile(x.std(axis=1), (283, 1))).T.reshape(len(planets), 283)\ny_pred_std = np.squeeze(np.tile(y_pred.std(axis=1), (283, 1))).T.reshape(len(planets), 283)\ny_pred_maxmin = np.squeeze(np.tile(y_pred.max(axis=1)/y_pred.min(axis=1), (283, 1))).T.reshape(len(planets), 283)\ny_pred_ratio = y_pred / np.squeeze(np.tile(y_pred.mean(axis=1), (283, 1))).T.reshape(len(planets), 283)\nx_pred_ratio = x / np.squeeze(np.tile(x.mean(axis=1), (283, 1))).T.reshape(len(planets), 283)\n\nx2 = np.stack((sensor, x_ratio, y_pred, y_pred_std, y_pred_ratio, y_pred_maxmin, m_noise, savgol_oscillations, x_meanratio, noise, signal_ratio), axis=-1)\n\ny2_pred = model_sigma.predict(x2.reshape(283*len(planets), np.shape(x2)[2])).reshape(len(planets), 283)","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:09:09.456451Z","iopub.execute_input":"2024-10-31T10:09:09.456821Z","iopub.status.idle":"2024-10-31T10:09:09.510399Z","shell.execute_reply.started":"2024-10-31T10:09:09.456785Z","shell.execute_reply":"2024-10-31T10:09:09.509021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = test_df.copy()[['planet_id', 'star']].set_index('planet_id')\n\nwl_cols = {}\nsigma_cols = {}\nfor i in range(283):\n    wl_cols['wl_' + str(i+1)] = y_pred[:, i]\n    sigma_cols['sigma_' + str(i+1)] = np.where(submission['star'].isin([0, 1]), y2_pred[:, i] * 1,\n        np.where(signal_delta_a < 0.007, y2_pred[:, i] * 1, y2_pred[:, i] * 1))\nsubmission = pd.concat([submission, pd.DataFrame(wl_cols, index=submission.index), pd.DataFrame(sigma_cols, index=submission.index)], axis=1)\nsubmission.drop(columns=['star'], inplace=True)\n\nsubmission.to_csv('submission.csv')\nsubmission","metadata":{"execution":{"iopub.status.busy":"2024-10-31T10:09:09.511713Z","iopub.execute_input":"2024-10-31T10:09:09.512154Z","iopub.status.idle":"2024-10-31T10:09:09.629741Z","shell.execute_reply.started":"2024-10-31T10:09:09.512111Z","shell.execute_reply":"2024-10-31T10:09:09.628065Z"},"trusted":true},"execution_count":null,"outputs":[]}]}