{"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"}],"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\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nimport random\nfrom tqdm import tqdm\nimport pickle\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-11-30T07:43:45.094427Z","iopub.execute_input":"2024-11-30T07:43:45.094836Z","iopub.status.idle":"2024-11-30T07:43:48.016727Z","shell.execute_reply.started":"2024-11-30T07:43:45.0948Z","shell.execute_reply":"2024-11-30T07:43:48.015525Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 任务分析：\n\n**训练数据**：\n* 行星数量：673颗行星。\n* 预测label：需要预测的目标数量为283个(光谱分解)——多输出回归。\n\n**测试数据**：\n* 约800颗隐藏的行星用于测试。","metadata":{}},{"cell_type":"code","source":"train_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv',index_col='planet_id')\n# test_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',index_col='planet_id')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')","metadata":{"execution":{"iopub.status.busy":"2024-11-30T07:43:55.294053Z","iopub.execute_input":"2024-11-30T07:43:55.294611Z","iopub.status.idle":"2024-11-30T07:43:55.616231Z","shell.execute_reply.started":"2024-11-30T07:43:55.294571Z","shell.execute_reply":"2024-11-30T07:43:55.615009Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## FGS1数据观察（暂时不用校准文件）\n\n**数据描述**\n\n* 每个文件包含135,000行图像，图像以0.1秒的时间间隔拍摄。每行是一个32x32的单波长图像。","metadata":{}},{"cell_type":"markdown","source":"取100468857号行星的FGS1数据进行观察","metadata":{}},{"cell_type":"code","source":"# 获取所有行星 ID 列表\nplanet_ids = train_adc_info.index.tolist()\n\n# 随机从中选择 9 个行星 ID\nrandom_ids = random.sample(planet_ids, 9)\nprint(\"从现有行星 ID 中随机生成的 9 个 ID：\", random_ids)\n\n# 使用 random_ids 遍历处理\nfor planet_id in random_ids:\n    try:\n        f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n        print(f\"行星 {planet_id} 的数据：\")\n        print(f_signal)\n    except FileNotFoundError:\n        print(f\"行星 {planet_id} 的数据文件未找到！\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:43:59.257607Z","iopub.execute_input":"2024-11-30T07:43:59.25806Z","iopub.status.idle":"2024-11-30T07:44:12.087301Z","shell.execute_reply.started":"2024-11-30T07:43:59.258021Z","shell.execute_reply":"2024-11-30T07:44:12.08609Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"将每一行数据恢复成32*32的像素点观察","metadata":{}},{"cell_type":"code","source":"# 遍历每个随机行星\nfor planet_id in random_ids:\n    try:\n        # 读取对应行星的 FGS1 数据\n        f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n        \n        # 可视化第 0 时刻和第 1 时刻的数据\n        _, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n        sns.heatmap(f_signal.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000)\n        ax1.set_aspect('equal')\n        ax1.set_title(f\"Planet {planet_id} - Time 0\")\n\n        sns.heatmap(f_signal.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000)\n        ax2.set_aspect('equal')\n        ax2.set_title(f\"Planet {planet_id} - Time 1\")\n\n        plt.suptitle(f'Comparison of FGS1 Data for Planet {planet_id}')\n        plt.show()\n        \n    except FileNotFoundError:\n        print(f\"行星 {planet_id} 的数据文件未找到！\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:44:23.010706Z","iopub.execute_input":"2024-11-30T07:44:23.011129Z","iopub.status.idle":"2024-11-30T07:44:34.844062Z","shell.execute_reply.started":"2024-11-30T07:44:23.011092Z","shell.execute_reply":"2024-11-30T07:44:34.842977Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n\nplanet_id = 100468857\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\nf_signal","metadata":{"execution":{"iopub.status.busy":"2024-11-27T11:43:20.616057Z","iopub.execute_input":"2024-11-27T11:43:20.616457Z","iopub.status.idle":"2024-11-27T11:43:21.099704Z","shell.execute_reply.started":"2024-11-27T11:43:20.616423Z","shell.execute_reply":"2024-11-27T11:43:21.098533Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"将每一行数据恢复成32\\*32的像素点观察","metadata":{}},{"cell_type":"code","source":"# 取100468857号行星0时刻和1时刻的FGS1数据进行比较\n_, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\nsns.heatmap(f_signal.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000)\nax1.set_aspect('equal')\nsns.heatmap(f_signal.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000)\nax2.set_aspect('equal')\nplt.suptitle('A pair of FGS1 images')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-27T11:43:21.101219Z","iopub.execute_input":"2024-11-27T11:43:21.101554Z","iopub.status.idle":"2024-11-27T11:43:21.870605Z","shell.execute_reply.started":"2024-11-27T11:43:21.101523Z","shell.execute_reply":"2024-11-27T11:43:21.86935Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"观察时序变化","metadata":{}},{"cell_type":"code","source":"# 定义滑动窗口平滑函数\ndef smooth_signal(signal, window=800):\n    \"\"\"\n    使用滑动窗口平滑信号。\n    参数：\n    - signal: ndarray, 累积信号。\n    - window: int, 滑动窗口大小。\n    返回：\n    - smooth_signal: ndarray, 平滑后的信号。\n    \"\"\"\n    return (signal[window:] - signal[:-window]) / window\n\n# 定义处理单个行星信号的函数\ndef process_planet_signal(planet_id):\n    \"\"\"\n    加载并处理指定行星的 FGS1 数据。\n    参数：\n    - planet_id: int, 行星 ID。\n    返回：\n    - net_signal: ndarray, 奇偶帧差分信号。\n    - smooth_signal: ndarray, 平滑后的信号。\n    \"\"\"\n    f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n    mean_signal = f_signal.values.mean(axis=1)  # 直接对每帧图像取平均值\n    net_signal = mean_signal[1::2] - mean_signal[0::2]  # 奇数帧 - 偶数帧\n    cum_signal = net_signal.cumsum()  # 累积信号\n    smoothed_signal = smooth_signal(cum_signal)  # 平滑信号\n    return net_signal, smoothed_signal\n\n\n# 创建子图，根据行星数量动态调整布局\nnum_planets = len(random_ids)\nfig, axes = plt.subplots(num_planets, 2, figsize=(12, 4 * num_planets))\nif num_planets == 1:\n    axes = [axes]  # 保证 axes 是二维列表形式，适配后续循环\n\n# 遍历每个随机行星 ID\nfor i, planet_id in enumerate(random_ids):\n    ax_signal, ax_smooth = axes[i]\n\n    # 处理每个行星的信号\n    net_signal, smoothed_signal = process_planet_signal(planet_id)\n\n    # 绘制原始信号\n    ax_signal.set_title(f'FGS1: time series of planet {planet_id} (raw)')\n    ax_signal.plot(net_signal, label='raw signal', alpha=0.7)\n    ax_signal.legend()\n\n    # 绘制平滑信号\n    ax_smooth.set_title(f'FGS1: time series of planet {planet_id} (smoothed)')\n    ax_smooth.plot(smoothed_signal, color='c', label='smoothened signal')\n    ax_smooth.legend()\n    ax_smooth.set_xlabel('time step')\n\n    # 可选：标注特定的时间点（这里根据 net_signal 自动计算关键点）\n    for time_step in [20500, 23500, 44000, 47000]:  # 示例特定时间点\n        ax_smooth.axvline(time_step, color='gray', linestyle='--', alpha=0.7)\n\n# 调整布局并显示图像\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:44:47.86876Z","iopub.execute_input":"2024-11-30T07:44:47.86914Z","iopub.status.idle":"2024-11-30T07:45:03.294746Z","shell.execute_reply.started":"2024-11-30T07:44:47.869108Z","shell.execute_reply":"2024-11-30T07:45:03.293243Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 4))\n\n#变化比较明显的一个行星\nplanet_id = 100468857\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\n#直接对每行图像取平均？是否合理？怎么改进？\nmean_signal = f_signal.values.mean(axis=1)\n# 奇数帧-偶数帧（观察数据好像两帧之间会有跳变）\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\n\n#滑动窗口平滑数据\nwindow=800\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax1.set_title('FGS1: time series of planet with strong signal')\nax1.plot(net_signal, label='raw signal')\nax1.legend()\nax3.plot(smooth_signal, color='c', label='smoothened signal')\nax3.legend()\nax3.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax3.axvline(time_step, color='gray')\n\n#变化没有那么明显的一个行星\nplanet_id = 4249337798\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\nmean_signal = f_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=800\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax2.set_title('FGS1: time series of planet with weak signal')\nax2.plot(net_signal, label='raw signal')\nax2.legend()\nax4.plot(smooth_signal, color='c', label='smoothened signal')\nax4.legend()\nax4.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax4.axvline(time_step, color='gray')\n\n# plt.suptitle('FGS1 time series', y=0.96)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-27T11:43:38.008493Z","iopub.execute_input":"2024-11-27T11:43:38.008844Z","iopub.status.idle":"2024-11-27T11:43:41.000319Z","shell.execute_reply.started":"2024-11-27T11:43:38.008809Z","shell.execute_reply":"2024-11-27T11:43:40.999126Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## AIRS 数据观察（暂时未校准）\n\n**数据描述**\n\n* 每个文件包含11,250行图像，图像以0.1秒的时间间隔拍摄。每行是一个32 x 356的单波长图像。","metadata":{}},{"cell_type":"markdown","source":"还是来观察100468857","metadata":{}},{"cell_type":"code","source":"planet_id = 100468857\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\na_signal","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:41.002398Z","iopub.execute_input":"2024-11-27T11:43:41.002875Z","iopub.status.idle":"2024-11-27T11:43:42.143828Z","shell.execute_reply.started":"2024-11-27T11:43:41.002805Z","shell.execute_reply":"2024-11-27T11:43:42.142409Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"将观察random_ids 的数据","metadata":{}},{"cell_type":"code","source":"# 遍历每个行星 ID\nfor planet_id in random_ids:\n    try:\n        # 加载 A 信号\n        a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n        print(f\"Planet ID: {planet_id}\")\n        print(a_signal)  # 输出数据内容或结构信息\n    except Exception as e:\n        print(f\"Failed to process Planet ID {planet_id}: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:47:21.874954Z","iopub.execute_input":"2024-11-30T07:47:21.875511Z","iopub.status.idle":"2024-11-30T07:47:40.470677Z","shell.execute_reply.started":"2024-11-30T07:47:21.875372Z","shell.execute_reply":"2024-11-30T07:47:40.469349Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 取100468857号行星0时刻和1时刻的AIRS数据进行比较\n_, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10))\n\nsns.heatmap(a_signal.iloc[0].values.reshape(32, 356), ax=ax1, vmin=0, vmax=52000)\nax1.set_title('AIRS Image at Time 0')\nax1.set_aspect('equal')\nax1.set_ylim(32, 0)\nax1.set_aspect('auto')\n\nsns.heatmap(a_signal.iloc[1].values.reshape(32, 356), ax=ax2, vmin=0, vmax=52000)\nax2.set_title('AIRS Image at Time 1')\nax2.set_aspect('equal')\nax2.set_ylim(32, 0)\nax2.set_aspect('auto') \n\nplt.suptitle('A Pair of AIRS Images')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-11-27T11:43:42.145154Z","iopub.execute_input":"2024-11-27T11:43:42.14551Z","iopub.status.idle":"2024-11-27T11:43:43.693344Z","shell.execute_reply.started":"2024-11-27T11:43:42.145476Z","shell.execute_reply":"2024-11-27T11:43:43.692003Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"对random_ids进行操作","metadata":{}},{"cell_type":"code","source":"# 遍历每个行星 ID\nfor planet_id in random_ids:\n    try:\n        # 加载 A 信号数据\n        a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n\n        # 创建热力图进行比较\n        _, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10))\n\n        sns.heatmap(a_signal.iloc[0].values.reshape(32, 356), ax=ax1, vmin=0, vmax=52000)\n        ax1.set_title(f'AIRS Image at Time 0 (Planet ID: {planet_id})')\n        ax1.set_ylim(32, 0)\n        ax1.set_aspect('auto')\n\n        sns.heatmap(a_signal.iloc[1].values.reshape(32, 356), ax=ax2, vmin=0, vmax=52000)\n        ax2.set_title(f'AIRS Image at Time 1 (Planet ID: {planet_id})')\n        ax2.set_ylim(32, 0)\n        ax2.set_aspect('auto')\n\n        plt.suptitle(f'A Pair of AIRS Images for Planet ID: {planet_id}')\n        plt.show()\n        \n    except Exception as e:\n        print(f\"Failed to process Planet ID {planet_id}: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:47:50.252315Z","iopub.execute_input":"2024-11-30T07:47:50.252694Z","iopub.status.idle":"2024-11-30T07:48:13.973233Z","shell.execute_reply.started":"2024-11-30T07:47:50.252662Z","shell.execute_reply":"2024-11-30T07:48:13.97199Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import InputLayer, Conv2D\nfrom astropy.stats import sigma_clip\n# 构建 CNN 模型，用于对死点进行插值\ndef build_cnn_model():\n    model = Sequential()\n    model.add(InputLayer(input_shape=(None, None, 1)))\n    model.add(Conv2D(32, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(64, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(1, kernel_size=(3, 3), padding='same'))\n    model.compile(optimizer='adam', loss='mean_squared_error')\n    return model\n\ndef mask_hot_dead_single_frame(frame, dead, dark):\n    \"\"\"\n    针对单个帧（即二维图像数据）进行热点和死点的掩盖与插值。\n\n    参数：\n    - frame: ndarray, 单个帧数据，形状为 (32, 356) 或类似\n    - dead: ndarray, 死点掩码，形状与 frame 相同\n    - dark: ndarray, 暗噪声数据，形状与 frame 相同\n    \n    返回值：\n    - 填充插值后的帧数据，形状与输入帧一致\n    \"\"\"\n    # 识别热点和死点，并将其转换为布尔类型\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask.astype(bool)  # 将 `hot` 转换为布尔类型\n    dead = dead.astype(bool)  # 将 `dead` 转换为布尔类型\n\n    # 合并热点与死点掩码\n    combined_mask = dead\n\n    # 使用掩码遮盖信号中的热点和死点位置\n    frame_masked = np.ma.masked_where(combined_mask, frame).filled(np.nan)  # 用 `NaN` 表示缺失值\n\n    # 使用周围点均值进行插值\n    filled_frame = np.copy(frame_masked)\n    nan_indices = np.argwhere(np.isnan(filled_frame))  # 找出所有的 `NaN` 位置\n\n    for (i, j) in nan_indices:\n        # 找出邻近的四个点并计算均值（跳过超出边界的点）\n        neighbors = []\n        if j - 1 >= 0:  # 左\n            neighbors.append(filled_frame[i, j - 1])\n        if j + 1 < frame.shape[1]:  # 右\n            neighbors.append(filled_frame[i, j + 1])\n        \n        # 计算邻居均值，排除 `NaN` 值\n        neighbors = [val for val in neighbors if not np.isnan(val)]\n        if len_neighbors := len(neighbors):  # 如果邻居存在非 NaN 的值\n            filled_frame[i, j] = sum(neighbors) / len_neighbors\n\n    return filled_frame","metadata":{"execution":{"iopub.status.busy":"2024-11-27T11:43:43.694998Z","iopub.execute_input":"2024-11-27T11:43:43.695474Z","iopub.status.idle":"2024-11-27T11:43:58.155822Z","shell.execute_reply.started":"2024-11-27T11:43:43.695426Z","shell.execute_reply":"2024-11-27T11:43:58.15448Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.157806Z","iopub.execute_input":"2024-11-27T11:43:58.158552Z","iopub.status.idle":"2024-11-27T11:43:58.753582Z","shell.execute_reply.started":"2024-11-27T11:43:58.158514Z","shell.execute_reply":"2024-11-27T11:43:58.752385Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/100468857/AIRS-CH0_calibration/dark.parquet')\ndead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/100468857/AIRS-CH0_calibration/dead.parquet')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.755171Z","iopub.execute_input":"2024-11-27T11:43:58.75563Z","iopub.status.idle":"2024-11-27T11:43:58.829253Z","shell.execute_reply.started":"2024-11-27T11:43:58.755582Z","shell.execute_reply":"2024-11-27T11:43:58.827995Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dark","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.83453Z","iopub.execute_input":"2024-11-27T11:43:58.834913Z","iopub.status.idle":"2024-11-27T11:43:58.872584Z","shell.execute_reply.started":"2024-11-27T11:43:58.834879Z","shell.execute_reply":"2024-11-27T11:43:58.871416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dead","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.874054Z","iopub.execute_input":"2024-11-27T11:43:58.874462Z","iopub.status.idle":"2024-11-27T11:43:58.910664Z","shell.execute_reply.started":"2024-11-27T11:43:58.874427Z","shell.execute_reply":"2024-11-27T11:43:58.909478Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dead.iloc[14:17,210:230]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.911974Z","iopub.execute_input":"2024-11-27T11:43:58.912367Z","iopub.status.idle":"2024-11-27T11:43:58.930405Z","shell.execute_reply.started":"2024-11-27T11:43:58.912334Z","shell.execute_reply":"2024-11-27T11:43:58.929319Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"p1=a_signal[1:2].values.reshape(32, 356)\np1 = p1.astype(float)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.931552Z","iopub.execute_input":"2024-11-27T11:43:58.931876Z","iopub.status.idle":"2024-11-27T11:43:58.945546Z","shell.execute_reply.started":"2024-11-27T11:43:58.931845Z","shell.execute_reply":"2024-11-27T11:43:58.944285Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from astropy.stats import sigma_clip\np1=mask_hot_dead_single_frame(p1, dead, dark)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.94701Z","iopub.execute_input":"2024-11-27T11:43:58.947386Z","iopub.status.idle":"2024-11-27T11:43:58.96253Z","shell.execute_reply.started":"2024-11-27T11:43:58.947352Z","shell.execute_reply":"2024-11-27T11:43:58.961206Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"对random_ids进行相同操作","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import InputLayer, Conv2D\nfrom astropy.stats import sigma_clip\n# 构建 CNN 模型，用于对死点进行插值\ndef build_cnn_model():\n    model = Sequential()\n    model.add(InputLayer(input_shape=(None, None, 1)))\n    model.add(Conv2D(32, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(64, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(1, kernel_size=(3, 3), padding='same'))\n    model.compile(optimizer='adam', loss='mean_squared_error')\n    return model\n\ndef mask_hot_dead_single_frame(frame, dead, dark):\n    \"\"\"\n    针对单个帧（即二维图像数据）进行热点和死点的掩盖与插值。\n\n    参数：\n    - frame: ndarray, 单个帧数据，形状为 (32, 356) 或类似\n    - dead: ndarray, 死点掩码，形状与 frame 相同\n    - dark: ndarray, 暗噪声数据，形状与 frame 相同\n    \n    返回值：\n    - 填充插值后的帧数据，形状与输入帧一致\n    \"\"\"\n    # 识别热点和死点，并将其转换为布尔类型\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask.astype(bool)  # 将 `hot` 转换为布尔类型\n    dead = dead.astype(bool)  # 将 `dead` 转换为布尔类型\n\n    # 合并热点与死点掩码\n    combined_mask = dead\n\n    # 使用掩码遮盖信号中的热点和死点位置\n    frame_masked = np.ma.masked_where(combined_mask, frame).filled(np.nan)  # 用 `NaN` 表示缺失值\n\n    # 使用周围点均值进行插值\n    filled_frame = np.copy(frame_masked)\n    nan_indices = np.argwhere(np.isnan(filled_frame))  # 找出所有的 `NaN` 位置\n\n    for (i, j) in nan_indices:\n        # 找出邻近的四个点并计算均值（跳过超出边界的点）\n        neighbors = []\n        if j - 1 >= 0:  # 左\n            neighbors.append(filled_frame[i, j - 1])\n        if j + 1 < frame.shape[1]:  # 右\n            neighbors.append(filled_frame[i, j + 1])\n        \n        # 计算邻居均值，排除 `NaN` 值\n        neighbors = [val for val in neighbors if not np.isnan(val)]\n        if len_neighbors := len(neighbors):  # 如果邻居存在非 NaN 的值\n            filled_frame[i, j] = sum(neighbors) / len_neighbors\n\n    return filled_frame\n    # 遍历所有行星 ID\nfor planet_id in random_ids:\n    try:\n        f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n        dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_calibration/dark.parquet')\n        dead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_calibration/dead.parquet')\n        p1=a_signal[1:2].values.reshape(32, 356)\n        p1 = p1.astype(float)\n        p1=mask_hot_dead_single_frame(p1, dead, dark)\n        # 取100468857号行星0时刻和1时刻的AIRS数据进行比较\n        _, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10))\n\n        sns.heatmap(p1.reshape(32, 356), ax=ax1, vmin=0, vmax=52000)\n        ax1.set_title(f'AIRS Image at Time 0 (Planet ID: {planet_id})')\n        ax1.set_aspect('equal')\n        ax1.set_ylim(32, 0)\n        ax1.set_aspect('auto')\n\n        sns.heatmap(p1.reshape(32, 356), ax=ax2, vmin=0, vmax=52000)\n        ax2.set_title(f'AIRS Image at Time 1 (Planet ID: {planet_id})')\n        ax2.set_aspect('equal')\n        ax2.set_ylim(32, 0)\n        ax2.set_aspect('auto') \n\n        \n        plt.suptitle(f'A Pair of AIRS Images for Planet ID: {planet_id}')\n        plt.show()\n    except Exception as e:\n        print(f\"Failed to process Planet ID {planet_id}: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T07:48:54.566939Z","iopub.execute_input":"2024-11-30T07:48:54.567315Z","iopub.status.idle":"2024-11-30T07:49:28.548918Z","shell.execute_reply.started":"2024-11-30T07:48:54.567283Z","shell.execute_reply":"2024-11-30T07:49:28.547488Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 取100468857号行星0时刻和1时刻的AIRS数据进行比较\n_, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10))\n\nsns.heatmap(p1.reshape(32, 356), ax=ax1, vmin=0, vmax=52000)\nax1.set_title('AIRS Image at Time 0')\nax1.set_aspect('equal')\nax1.set_ylim(32, 0)\nax1.set_aspect('auto')\n\nsns.heatmap(p1.reshape(32, 356), ax=ax2, vmin=0, vmax=52000)\nax2.set_title('AIRS Image at Time 1')\nax2.set_aspect('equal')\nax2.set_ylim(32, 0)\nax2.set_aspect('auto') \n\nplt.suptitle('A Pair of AIRS Images')\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:43:58.96389Z","iopub.execute_input":"2024-11-27T11:43:58.964295Z","iopub.status.idle":"2024-11-27T11:44:00.576421Z","shell.execute_reply.started":"2024-11-27T11:43:58.96424Z","shell.execute_reply":"2024-11-27T11:44:00.575171Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"这个图里面似乎就有一个坏点（）","metadata":{}},{"cell_type":"markdown","source":"生成点线图","metadata":{}},{"cell_type":"markdown","source":"但是缺数据，没有没有遮掩的数据","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nfrom astropy.stats import sigma_clip\n\n# 构建 CNN 模型，用于对死点进行插值\ndef build_cnn_model():\n    model = Sequential()\n    model.add(InputLayer(input_shape=(None, None, 1)))\n    model.add(Conv2D(32, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(64, kernel_size=(3, 3), padding='same', activation='relu'))\n    model.add(Conv2D(1, kernel_size=(3, 3), padding='same'))\n    model.compile(optimizer='adam', loss='mean_squared_error')\n    return model\n\n# 掩盖热点与死点，并进行插值的函数\ndef mask_hot_dead_single_frame(frame, dead, dark):\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask.astype(bool)\n    dead = dead.astype(bool)\n    \n    combined_mask = dead\n    frame_masked = np.ma.masked_where(combined_mask, frame).filled(np.nan)\n    \n    filled_frame = np.copy(frame_masked)\n    nan_indices = np.argwhere(np.isnan(filled_frame))\n    \n    for (i, j) in nan_indices:\n        neighbors = []\n        if j - 1 >= 0:  # 左\n            neighbors.append(filled_frame[i, j - 1])\n        if j + 1 < frame.shape[1]:  # 右\n            neighbors.append(filled_frame[i, j + 1])\n        \n        # 计算邻居均值，排除 NaN 值\n        neighbors = [val for val in neighbors if not np.isnan(val)]\n        if len(neighbors) > 0:  # 如果邻居存在非 NaN 的值\n            filled_frame[i, j] = sum(neighbors) / len(neighbors)\n\n    return filled_frame\n\n# 循环处理 random_ids 中的每个行星\nfor planet_id in random_ids:\n    # 读取该行星的 AIRC 和 AIRS 信号数据\n    f_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n    a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n    dark = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_calibration/dark.parquet')\n    dead = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_calibration/dead.parquet')\n    \n    # 假设我们选择第1帧进行处理\n    p1 = a_signal.iloc[1:2].values.reshape(32, 356)\n    p1 = p1.astype(float)\n    \n    # 处理热点和死点\n    p1 = mask_hot_dead_single_frame(p1, dead, dark)\n    \n    # 计算亮度变化\n    frame_brightness = np.sum(p1, axis=(0, 1))  # 计算该帧的亮度总和\n    \n    # 假设有一个基准亮度 baseline_flux（没有行星遮挡时的亮度），这里我们做一个简化处理,但是出问题了\n    baseline_flux = np.sum(a_signal.iloc[0].values.reshape(32, 356))  # 使用第一帧作为基准亮度\n\n    print(f\"frame_bringhtness = {frame_brightness}, baseline_flux = {baseline_flux}\")\n    \n    # 计算亮度变化\n    flux_change = baseline_flux - frame_brightness  # 计算基准帧和当前帧的亮度差\n    \n    if flux_change > 0 and baseline_flux > 0:\n        Rp_Rs = np.sqrt(flux_change / baseline_flux)\n    else:\n        Rp_Rs = np.nan  # 或者其他处理方式\n        print(f\"Invalid flux values: flux_change = {flux_change}, baseline_flux = {baseline_flux}\")\n    \n    # 根据亮度变化计算 Rp/Rs\n    Rp_Rs = np.sqrt(flux_change / baseline_flux)  # 根据透过率模型估算行星半径与恒星半径之比\n    \n    # 打印计算的 Rp/Rs 值\n    print(f\"Planet ID: {planet_id}, Rp/Rs: {Rp_Rs}\")\n    \n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-30T08:17:35.470186Z","iopub.execute_input":"2024-11-30T08:17:35.471164Z","iopub.status.idle":"2024-11-30T08:17:51.409598Z","shell.execute_reply.started":"2024-11-30T08:17:35.471114Z","shell.execute_reply":"2024-11-30T08:17:51.407841Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"和上面fgs1的数据处理基本一致，观察时序图","metadata":{}},{"cell_type":"code","source":"_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 4))\n\n#还是取上面两个行星\nplanet_id = 100468857\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n\nmean_signal = a_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=80\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax1.set_title('AIRS-CH0 time series of planet 100468857')\nax1.plot(net_signal, label='raw signal')\nax1.legend()\nax3.plot(smooth_signal, color='c', label='smoothened signal')\nax3.legend()\nax3.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax3.axvline(time_step* 11250 // 135000, color='gray')\n    \n\nplanet_id = 4249337798\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n\nmean_signal = a_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=80\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax2.set_title('AIRS-CH0 time series of planet 4249337798')\nax2.plot(net_signal, label='raw signal')\nax2.legend()\nax4.plot(smooth_signal, color='c', label='smoothened signal')\nax4.legend()\nax4.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax4.axvline(time_step* 11250 // 135000, color='gray')\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:44:00.578245Z","iopub.execute_input":"2024-11-27T11:44:00.578639Z","iopub.status.idle":"2024-11-27T11:44:05.158863Z","shell.execute_reply.started":"2024-11-27T11:44:00.578605Z","shell.execute_reply":"2024-11-27T11:44:05.157594Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"与上文一样，用random_ids来替换","metadata":{}},{"cell_type":"code","source":"# 创建子图，根据行星数量动态调整布局\nnum_planets = len(random_ids)\nfig, axes = plt.subplots(num_planets, 2, figsize=(12, 4 * num_planets))\nif num_planets == 1:\n    axes = [axes]  # 保证 axes 是二维列表形式，适配后续循环\n\n# 遍历每个随机行星 ID\nfor i, planet_id in enumerate(random_ids):\n    ax_signal, ax_smooth = axes[i]\n\n    a_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n\n    mean_signal = a_signal.values.mean(axis=1)\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n    cum_signal = net_signal.cumsum()\n    window=80\n    smooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax1.set_title('AIRS-CH0 time series of planet 100468857')\nax1.plot(net_signal, label='raw signal')\nax1.legend()\nax3.plot(smooth_signal, color='c', label='smoothened signal')\nax3.legend()\nax3.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax3.axvline(time_step* 11250 // 135000, color='gray')\n    \n\nplanet_id = 4249337798\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\n\nmean_signal = a_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow=80\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\nax2.set_title('AIRS-CH0 time series of planet 4249337798')\nax2.plot(net_signal, label='raw signal')\nax2.legend()\nax4.plot(smooth_signal, color='c', label='smoothened signal')\nax4.legend()\nax4.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax4.axvline(time_step* 11250 // 135000, color='gray')\n\nplt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"为所有673个训练行星读取FGS1数据和AIRS-CH0数据。\n\n由于数据集无法完全放入RAM，我们仅保留每个行星的两个**一维时间序列**。即：\n\n* 从FGS1数据中提取的每个行星67500步的时间序列\n\n* 从AIRS-CH0数据中提取的每个行星5625步的时间序列\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport polars as pl\nfrom tqdm import tqdm\n\ndef f_read_and_preprocess(dataset, adc_info, planet_ids):\n    \n#     读取所有行星ID的FGS1文件并提取时间序列。\n\n#     参数：\n#     dataset：'train' 或 'test'\n#     adc_info：元数据数据框，可能是 train_adc_info 或 test_adc_info\n#     planet_ids：行星ID列表\n\n#     返回：\n#     每个行星ID对应一行的 数据框，每行包含67500个值\n\n    f_raw_train = np.full((len(planet_ids), 67500), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        f_signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_signal.parquet')\n        mean_signal = f_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / 1024 \n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        f_raw_train[i] = net_signal\n    return f_raw_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:44:05.160193Z","iopub.execute_input":"2024-11-27T11:44:05.160527Z","iopub.status.idle":"2024-11-27T11:44:05.168582Z","shell.execute_reply.started":"2024-11-27T11:44:05.160495Z","shell.execute_reply":"2024-11-27T11:44:05.167252Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nf_raw_train = f_read_and_preprocess('train', train_adc_info, train_labels.index)\nwith open('f_raw_train.pickle', 'wb') as f:\n    pickle.dump(f_raw_train, f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:44:05.171139Z","iopub.execute_input":"2024-11-27T11:44:05.172001Z","iopub.status.idle":"2024-11-27T11:54:09.65876Z","shell.execute_reply.started":"2024-11-27T11:44:05.171905Z","shell.execute_reply":"2024-11-27T11:54:09.656302Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport polars as pl\nfrom tqdm import tqdm\n\ndef a_read_and_preprocess(dataset, adc_info, planet_ids):\n    \n#     读取所有行星ID的AIRS-CH0文件并提取时间序列。\n#     参数：\n#     dataset：'train' 或 'test'\n#     adc_info：元数据数据框，可能是 train_adc_info 或 test_adc_info\n#     planet_ids：行星ID列表\n\n#     返回：\n#     每个行星ID对应一行的 数据框，每行包含5625个值\n\n    a_raw_train = np.full((len(planet_ids), 5625), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet')\n        mean_signal = signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / (32*356) # 对32*356个像素求均值\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        a_raw_train[i] = net_signal\n    return a_raw_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:54:09.661625Z","iopub.execute_input":"2024-11-27T11:54:09.662191Z","iopub.status.idle":"2024-11-27T11:54:09.673647Z","shell.execute_reply.started":"2024-11-27T11:54:09.662129Z","shell.execute_reply":"2024-11-27T11:54:09.672448Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\na_raw_train = a_read_and_preprocess('train', train_adc_info, train_labels.index)\nwith open('a_raw_train.pickle', 'wb') as f:\n    pickle.dump(a_raw_train, f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T11:54:09.675055Z","iopub.execute_input":"2024-11-27T11:54:09.675472Z","iopub.status.idle":"2024-11-27T12:10:36.039536Z","shell.execute_reply.started":"2024-11-27T11:54:09.675439Z","shell.execute_reply":"2024-11-27T12:10:36.03785Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"全数据观察","metadata":{}},{"cell_type":"code","source":"f_raw_train.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T12:10:36.04277Z","iopub.execute_input":"2024-11-27T12:10:36.043233Z","iopub.status.idle":"2024-11-27T12:10:36.05261Z","shell.execute_reply.started":"2024-11-27T12:10:36.043189Z","shell.execute_reply":"2024-11-27T12:10:36.051186Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(6, 2))\nplt.plot(f_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step, color='gray')\nplt.xlabel('time step')\nplt.title('FGS1: Overall mean')\nplt.show()\n\nplt.figure(figsize=(6, 2))\nplt.plot(a_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step * 11250 // 135000, color='gray')\nplt.xlabel('time step')\nplt.title('AIRS-CH0: Overall mean')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-27T12:10:36.054599Z","iopub.execute_input":"2024-11-27T12:10:36.054976Z","iopub.status.idle":"2024-11-27T12:10:36.658594Z","shell.execute_reply.started":"2024-11-27T12:10:36.054942Z","shell.execute_reply":"2024-11-27T12:10:36.657191Z"}},"outputs":[],"execution_count":null}]}