{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"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"},{"sourceId":121871,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":102548,"modelId":126770}],"dockerImageVersionId":30746,"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\nimport itertools\nimport os\nimport glob \nfrom astropy.stats import sigma_clip\n\nfrom tqdm import tqdm\n\ndef ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal\ndef mask_hot_dead(signal, dead, dark):\n    hot = sigma_clip(\n        dark, sigma=5, maxiters=5\n    ).mask\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal)\n    return signal\ndef 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\ndef clean_dark(signal, dead, dark, dt):\n    dark = np.ma.masked_where(dead, dark)\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n\n    signal -= dark* dt[:, np.newaxis, np.newaxis]\n    return signal\ndef get_cds(signal):\n    cds = signal[:,1::2,:,:] - signal[:,::2,:,:]\n    return cds\n\ndef bin_obs(cds_signal,binning):\n    cds_transposed = cds_signal.transpose(0,1,3,2)\n    cds_binned = np.zeros((cds_transposed.shape[0], cds_transposed.shape[1]//binning, cds_transposed.shape[2], cds_transposed.shape[3]))\n    for i in range(cds_transposed.shape[1]//binning):\n        cds_binned[:,i,:,:] = np.sum(cds_transposed[:,i*binning:(i+1)*binning,:,:], axis=1)\n    return cds_binned\ndef correct_flat_field(flat,dead, signal):\n    flat = flat.transpose(1, 0)\n    dead = dead.transpose(1, 0)\n    flat = np.ma.masked_where(dead, flat)\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    signal = signal / flat\n    return signal\ndef get_index(files, CHUNKS_SIZE):\n    index = []\n    for file in files:\n        # 使用 os.path.split 代替手动的 '/' 分割，兼容不同操作系统\n        file_name = os.path.basename(file)\n        if file_name.startswith('AIRS-CH0') and file_name.endswith('signal.parquet'):\n            file_index = os.path.basename(os.path.dirname(file))  # 获取父目录名作为行星索引\n            # 确保文件存在\n            signal_file = os.path.join(os.path.dirname(file), 'AIRS-CH0_signal.parquet')\n            if os.path.exists(signal_file):\n                index.append(int(file_index))\n    \n    index = np.array(index)\n    index = np.sort(index)  # 对索引进行排序\n    # 将索引数组按 CHUNKS_SIZE 分割\n    index = np.array_split(index, len(index) // CHUNKS_SIZE)\n    \n    return index\n\n","metadata":{"execution":{"iopub.status.busy":"2024-09-27T14:30:42.267017Z","iopub.execute_input":"2024-09-27T14:30:42.267523Z","iopub.status.idle":"2024-09-27T14:30:43.31209Z","shell.execute_reply.started":"2024-09-27T14:30:42.267474Z","shell.execute_reply":"2024-09-27T14:30:43.310636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm  # 进度条库\nimport glob\n\n# 声明全局变量来存储处理结果\nglobal_air_clean_data = None\nglobal_fgs_clean_data = None\n\ndef process_planet_data(path_folder, dataset_type='train', num_planets=None, CHUNKS_SIZE=10, DO_MASK=True, DO_THE_NL_CORR=False, DO_DARK=True, DO_FLAT=True, TIME_BINNING=True):\n    \"\"\"\n    处理行星数据的主函数。根据提供的参数，访问 `train` 或 `test` 文件夹中的行星数据。\n    \n    参数：\n    - path_folder: 数据集的根目录路径。\n    - dataset_type: 数据集类型 ('train' 或 'test')。\n    - num_planets: 要访问的行星数量。如果为 None，则处理全部行星。\n    - CHUNKS_SIZE: 每个数据块的大小。\n    - DO_MASK: 是否应用掩膜。\n    - DO_THE_NL_CORR: 是否进行非线性校正。\n    - DO_DARK: 是否进行暗信号校正。\n    - DO_FLAT: 是否进行平场校正。\n    - TIME_BINNING: 是否进行时间合并（减少空间占用）。\n    \"\"\"\n    \n    global global_air_clean_data, global_fgs_clean_data  # 使用全局变量\n    \n    # 动态构建文件路径\n    files = glob.glob(os.path.join(path_folder, f'{dataset_type}/', '*/*'))\n\n    # 如果未定义 num_planets，则处理所有行星\n    if num_planets is None:\n        num_planets = len(files)\n\n    # 获取指定数量的行星数据索引\n    index = get_index(files[:num_planets], CHUNKS_SIZE)\n\n    # 加载行星相关的辅助信息\n    train_adc_info = pd.read_csv(os.path.join(path_folder, f'{dataset_type}_adc_info.csv'))\n    train_adc_info = train_adc_info.set_index('planet_id')\n    axis_info = pd.read_parquet(os.path.join(path_folder, 'axis_info.parquet'))\n\n    # 定义裁剪区间\n    cut_inf, cut_sup = 39, 321\n    l = cut_sup - cut_inf\n\n    all_air_data = []  # 存储全部处理后的 AIRS 数据\n    all_fgs_data = []  # 存储全部处理后的 FGS1 数据\n\n    for n, index_chunk in enumerate(tqdm(index)):\n        AIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, l))\n        FGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))\n        \n        for i in range(CHUNKS_SIZE):\n            try:\n                # 读取并处理 AIRS-CH0 信号\n                df = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/AIRS-CH0_signal.parquet'))\n                signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 356))\n                gain = train_adc_info['AIRS-CH0_adc_gain'].loc[index_chunk[i]]\n                offset = train_adc_info['AIRS-CH0_adc_offset'].loc[index_chunk[i]]\n                signal = ADC_convert(signal, gain, offset)\n                dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n                dt_airs[1::2] += 0.1\n                chopped_signal = signal[:, :, cut_inf:cut_sup]\n                del signal, df\n\n                # 清理 AIRS 数据\n                flat = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n                dark = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/AIRS-CH0_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n                dead_airs = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/AIRS-CH0_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n                linear_corr = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/AIRS-CH0_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n                if DO_MASK:\n                    chopped_signal = mask_hot_dead(chopped_signal, dead_airs, dark)\n                    AIRS_CH0_clean[i] = chopped_signal\n                else:\n                    AIRS_CH0_clean[i] = chopped_signal\n\n                if DO_THE_NL_CORR:\n                    linear_corr_signal = apply_linear_corr(linear_corr, AIRS_CH0_clean[i])\n                    AIRS_CH0_clean[i, :, :, :] = linear_corr_signal\n                del linear_corr\n\n                if DO_DARK:\n                    cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark, dt_airs)\n                    AIRS_CH0_clean[i] = cleaned_signal\n                del dark\n\n                # 读取并处理 FGS1 信号\n                df = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/FGS1_signal.parquet'))\n                fgs_signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 32))\n\n                FGS1_gain = train_adc_info['FGS1_adc_gain'].loc[index_chunk[i]]\n                FGS1_offset = train_adc_info['FGS1_adc_offset'].loc[index_chunk[i]]\n\n                fgs_signal = ADC_convert(fgs_signal, FGS1_gain, FGS1_offset)\n                dt_fgs1 = np.ones(len(fgs_signal)) * 0.1\n                dt_fgs1[1::2] += 0.1\n                chopped_FGS1 = fgs_signal\n                del fgs_signal, df\n\n                # 清理 FGS1 数据\n                flat = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/FGS1_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n                dark = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/FGS1_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 32))\n                dead_fgs1 = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/FGS1_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 32))\n                linear_corr = pd.read_parquet(os.path.join(path_folder, f'{dataset_type}/{index_chunk[i]}/FGS1_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 32))\n\n                if DO_MASK:\n                    chopped_FGS1 = mask_hot_dead(chopped_FGS1, dead_fgs1, dark)\n                    FGS1_clean[i] = chopped_FGS1\n                else:\n                    FGS1_clean[i] = chopped_FGS1\n\n                if DO_THE_NL_CORR:\n                    linear_corr_signal = apply_linear_corr(linear_corr, FGS1_clean[i])\n                    FGS1_clean[i, :, :, :] = linear_corr_signal\n                del linear_corr\n\n                if DO_DARK:\n                    cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark, dt_fgs1)\n                    FGS1_clean[i] = cleaned_signal\n                del dark\n\n            except FileNotFoundError as e:\n                print(f\"File not found: {e}. Skipping this planet.\")\n                continue  # 跳过当前循环，继续处理下一个行星\n\n        # 保存 AIRS 和 FGS1 数据到列表中\n        AIRS_cds = get_cds(AIRS_CH0_clean)\n        FGS1_cds = get_cds(FGS1_clean)\n        \n        if TIME_BINNING:\n            AIRS_cds_binned = bin_obs(AIRS_cds, binning=30)\n            FGS1_cds_binned = bin_obs(FGS1_cds, binning=30 * 12)\n        else:\n            AIRS_cds = AIRS_cds.transpose(0, 1, 3, 2)  # 确保平场一致性\n            AIRS_cds_binned = AIRS_cds\n            FGS1_cds = FGS1_cds.transpose(0, 1, 3, 2)\n            FGS1_cds_binned = FGS1_cds\n\n        # 将处理后的数据存储到列表中\n        all_air_data.append(AIRS_cds_binned)\n        all_fgs_data.append(FGS1_cds_binned)\n\n        del AIRS_cds_binned\n        del FGS1_cds_binned\n\n    # 将处理后的全部数据存储到全局变量中\n    global_air_clean_data = np.concatenate(all_air_data, axis=0)\n    global_fgs_clean_data = np.concatenate(all_fgs_data, axis=0)\n\n    print(f\"所有 {num_planets} 颗行星的数据处理完成，存储到全局变量 global_air_clean_data 和 global_fgs_clean_data 中。\")                                                          \n       \n","metadata":{"execution":{"iopub.status.busy":"2024-09-27T14:30:43.31956Z","iopub.execute_input":"2024-09-27T14:30:43.319999Z","iopub.status.idle":"2024-09-27T14:30:43.361343Z","shell.execute_reply.started":"2024-09-27T14:30:43.319963Z","shell.execute_reply":"2024-09-27T14:30:43.35986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom scipy.signal import find_peaks\nfrom scipy.fftpack import fft\nfrom scipy.ndimage import gaussian_filter1d\nfrom scipy.signal import savgol_filter\n\ndef process_spectral_data(data_airs_ch0, data_fgs1, axis_info, method='savgol'):\n    # 数据预处理函数 (标准化)\n    def preprocess_data(data):\n        num_observations = data.shape[0]\n        preprocessed_data = []\n        for i in range(num_observations):\n            observation = data[i]\n            observation_normalized = (observation - np.mean(observation)) / np.std(observation)\n            preprocessed_data.append(observation_normalized)\n        return np.array(preprocessed_data)\n\n    # 去噪函数\n    def denoise_spectra(data, method='savgol'):\n        num_observations, num_time_points, _, _ = data.shape\n        denoised_data = np.zeros_like(data)\n        for i in range(num_observations):\n            for j in range(num_time_points):\n                if method == 'savgol':\n                    denoised_data[i, j] = savgol_filter(data[i, j], window_length=11, polyorder=2, axis=0)\n                elif method == 'gaussian':\n                    denoised_data[i, j] = gaussian_filter1d(data[i, j], sigma=2, axis=0)\n        return denoised_data\n\n    # 提取光曲线函数\n    def extract_light_curve(data):\n        num_observations, num_time_points = data.shape[:2]\n        light_curves = []\n        for i in range(num_observations):\n            light_curve = np.sum(data[i], axis=(1, 2))\n            light_curves.append(light_curve)\n        return np.array(light_curves)\n\n    # 提取时间序列特征的函数\n    def extract_time_series_features(light_curves, axis_info):\n        features = []\n        for i in range(light_curves.shape[0]):\n            observation_features = {}\n            time_series = light_curves[i].flatten()\n\n            # 基本统计特征\n            observation_features['mean'] = np.mean(time_series)\n            observation_features['std'] = np.std(time_series)\n            observation_features['max'] = np.max(time_series)\n            observation_features['min'] = np.min(time_series)\n\n            # 傅里叶变换特征\n            fft_features = np.abs(fft(time_series))\n            observation_features['fft_mean'] = np.mean(fft_features)\n            observation_features['fft_max'] = np.max(fft_features)\n\n            # 滑动窗口统计特征\n            window_size = 10\n            rolling_mean = pd.Series(time_series).rolling(window=window_size).mean().mean()\n            rolling_std = pd.Series(time_series).rolling(window=window_size).std().mean()\n            observation_features['rolling_mean'] = rolling_mean\n            observation_features['rolling_std'] = rolling_std\n\n            # 峰值特征\n            peaks, _ = find_peaks(time_series, height=0)\n            observation_features['num_peaks'] = len(peaks)\n            observation_features['mean_peak_height'] = np.mean(time_series[peaks]) if len(peaks) > 0 else 0\n\n            # 结合 axis_info 中的特征\n            observation_features['AIRS-CH0-axis0-h'] = axis_info['AIRS-CH0-axis0-h'].iloc[i]\n            observation_features['AIRS-CH0-axis2-um'] = axis_info['AIRS-CH0-axis2-um'].iloc[i]\n            observation_features['AIRS-CH0-integration_time'] = axis_info['AIRS-CH0-integration_time'].iloc[i]\n            observation_features['FGS1-axis0-h'] = axis_info['FGS1-axis0-h'].iloc[i]\n\n            features.append(observation_features)\n        \n        return pd.DataFrame(features)\n\n    # 预处理数据\n    preprocessed_airs_ch0 = preprocess_data(data_airs_ch0)\n    preprocessed_fgs1 = preprocess_data(data_fgs1)\n    \n    # 去噪\n    denoised_airs_ch0 = denoise_spectra(preprocessed_airs_ch0, method=method)\n    denoised_fgs1 = denoise_spectra(preprocessed_fgs1, method=method)\n    \n    # 提取光曲线\n    light_curve_airs_ch0 = extract_light_curve(denoised_airs_ch0)\n    light_curve_fgs1 = extract_light_curve(denoised_fgs1)\n    \n    # 提取时间序列特征\n    features_airs_ch0 = extract_time_series_features(light_curve_airs_ch0, axis_info)\n    features_fgs1 = extract_time_series_features(light_curve_fgs1, axis_info)\n    \n    # 合并两个数据集的特征\n    test_features = pd.concat([features_airs_ch0, features_fgs1], axis=1)\n    test_features = test_features.fillna(method='ffill')\n    \n    return test_features","metadata":{"execution":{"iopub.status.busy":"2024-09-27T14:30:43.363986Z","iopub.execute_input":"2024-09-27T14:30:43.364696Z","iopub.status.idle":"2024-09-27T14:30:44.092203Z","shell.execute_reply.started":"2024-09-27T14:30:43.364649Z","shell.execute_reply":"2024-09-27T14:30:44.090831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestRegressor\nimport numpy as np\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Conv1D, Dense, Flatten, MaxPooling1D\nfrom sklearn.metrics import mean_squared_error, r2_score\n\n\ndef add_gaussian_noise(data, noise_factor=0.00021):\n    noise = np.random.normal(0, noise_factor, data.shape)\n    return data + noise\n    from tensorflow.keras.models import load_model\nfrom tensorflow.keras.models import load_model\n\ndef custom_loss(y_true, y_pred):\n    mse = tf.reduce_mean(tf.square(y_true - y_pred), axis=-1)\n    penalty = tf.reduce_sum(tf.maximum(-y_pred, 0.0), axis=-1)\n    return mse + penalty * 0.01\n\ndef load_ensemble(n_models, path_prefix, custom_objects=None):\n    models = []\n    for i in range(n_models):\n        model = load_model(f'{path_prefix}/model_{i}.h5', custom_objects=custom_objects)\n        models.append(model)\n    return models\n\ndef predict_ensemble(models, x_test):\n    predictions = np.array([model.predict(x_test).astype(np.float64) for model in models])\n    mu_pred = np.mean(predictions, axis=0).astype(np.float64)\n    return mu_pred\ndef postprocessing(predictions, index, sigma_pred):\n    processed_df = pd.DataFrame({\n        \n        'planet_id': index,\n        'predictions': predictions.flatten(),\n        'sigma_pred': sigma_pred.flatten()\n    })\n    \n    # 你可以根据需要在这里进行更多处理\n    processed_df = processed_df.sort_values(by='planet_id')\n    return processed_df\n","metadata":{"execution":{"iopub.status.busy":"2024-09-27T14:30:44.095614Z","iopub.execute_input":"2024-09-27T14:30:44.096054Z","iopub.status.idle":"2024-09-27T14:30:49.515454Z","shell.execute_reply.started":"2024-09-27T14:30:44.096018Z","shell.execute_reply":"2024-09-27T14:30:49.514195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH = '/kaggle/input/ariel-data-challenge-2024'\npath_prefix = '/kaggle/input/11111/tensorflow2/default/1'\n# 生成提交文件\ntest_adc_info = pd.read_csv(f'{PATH}/test_adc_info.csv', index_col='planet_id')\nsample_submission = pd.read_csv(f'{PATH}/sample_submission.csv', index_col='planet_id')\naxis_info = pd.read_parquet(f'{PATH}/axis_info.parquet')\ntest_adc_info_index = test_adc_info.index\nprocess_planet_data(\n    path_folder='/kaggle/input/ariel-data-challenge-2024/',\n    dataset_type='test',  # 也可以是 'test' \n    CHUNKS_SIZE=1,   # 每个块的大小为 20\n    DO_MASK=True, \n    DO_DARK=True, \n    DO_FLAT=True, \n    TIME_BINNING=True\n)\na_raw_test = global_air_clean_data\nf_raw_test = global_fgs_clean_data\nfirst_dimension_length = a_raw_test.shape[0]\n# 提取与 first_dimension_length 一样长度的索引\ntest_adc_info_index = test_adc_info.index[:first_dimension_length]\ntest_features = process_spectral_data(a_raw_test, f_raw_test, axis_info, method='savgol')\n\nn_models = 5\ncustom_objects = {'custom_loss': custom_loss}\nloaded_models = load_ensemble(n_models, path_prefix, custom_objects=custom_objects)    \ntest_features = np.expand_dims(test_features, axis=-1)\ntest_pred = predict_ensemble(loaded_models, test_features)\n#wl\npred = pd.DataFrame(test_pred, columns=[f'wl{i}' for i in range(1, 284)])\npred.insert(0, 'planet_id', test_adc_info_index)\n#unc\nnum_rows = test_pred.shape[0]\ncolumn_names = [f'unc_{i+1}' for i in range(283)]\ndf = pd.DataFrame(np.full((num_rows, 283), 0.00021, dtype='float64'), columns=column_names)\nsubmission = pd.concat([pred,df], axis=1)\nsubmission.to_csv('submission.csv',index = False)","metadata":{"execution":{"iopub.status.busy":"2024-09-27T14:33:16.507257Z","iopub.execute_input":"2024-09-27T14:33:16.507961Z","iopub.status.idle":"2024-09-27T14:33:39.296273Z","shell.execute_reply.started":"2024-09-27T14:33:16.507916Z","shell.execute_reply":"2024-09-27T14:33:39.29494Z"},"trusted":true},"execution_count":null,"outputs":[]}]}