{"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":60182,"databundleVersionId":6787572,"sourceType":"competition"}],"dockerImageVersionId":30804,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:41:21.938617Z","iopub.execute_input":"2025-12-04T04:41:21.939103Z","iopub.status.idle":"2025-12-04T04:41:22.734484Z","shell.execute_reply.started":"2025-12-04T04:41:21.93905Z","shell.execute_reply":"2025-12-04T04:41:22.733144Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df=pd.read_csv(\"/kaggle/input/earthquake-prediction/train.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:41:22.736674Z","iopub.execute_input":"2025-12-04T04:41:22.737216Z","iopub.status.idle":"2025-12-04T04:44:05.053217Z","shell.execute_reply.started":"2025-12-04T04:41:22.737177Z","shell.execute_reply":"2025-12-04T04:44:05.051549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df.head(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:05.055193Z","iopub.execute_input":"2025-12-04T04:44:05.056128Z","iopub.status.idle":"2025-12-04T04:44:05.10686Z","shell.execute_reply.started":"2025-12-04T04:44:05.056035Z","shell.execute_reply":"2025-12-04T04:44:05.105567Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Eksik değerlerin kontrolü\nprint(df.isnull().sum())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:05.108368Z","iopub.execute_input":"2025-12-04T04:44:05.108791Z","iopub.status.idle":"2025-12-04T04:44:05.916162Z","shell.execute_reply.started":"2025-12-04T04:44:05.108732Z","shell.execute_reply":"2025-12-04T04:44:05.914679Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Temel istatistiksel bilgiler için\nprint(df.describe())\n\n# TTF'nin histogramı\nimport matplotlib.pyplot as plt\n\nplt.hist(df['ttf'], bins=50, alpha=0.7)\nplt.xlabel('Time to Failure (ttf)')\nplt.ylabel('Frequency')\nplt.title('TTF Distribution')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:05.919397Z","iopub.execute_input":"2025-12-04T04:44:05.920061Z","iopub.status.idle":"2025-12-04T04:44:52.086923Z","shell.execute_reply.started":"2025-12-04T04:44:05.919993Z","shell.execute_reply":"2025-12-04T04:44:52.085701Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(df.columns)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:52.088271Z","iopub.execute_input":"2025-12-04T04:44:52.088599Z","iopub.status.idle":"2025-12-04T04:44:52.095916Z","shell.execute_reply.started":"2025-12-04T04:44:52.088568Z","shell.execute_reply":"2025-12-04T04:44:52.094618Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ttf ve EP_disp'in türevlerini (değişim hızları) hesaplayalım\ndf['EP_disp_derivative'] = df['EP_disp'].diff()  # Yer değiştirme türevi\ndf['ttf_derivative'] = df['ttf'].diff()  # Zaman türevi\n\n# Artış oranlarını hesaplama\ndf['EP_disp_growth_rate'] = df['EP_disp'].pct_change()  # Yüzde değişim oranı\ndf['ttf_growth_rate'] = df['ttf'].pct_change()\n\n# Trendlerin ilk birkaç satırını kontrol edelim\nprint(df[['EP_disp', 'EP_disp_derivative', 'EP_disp_growth_rate', \n          'ttf', 'ttf_derivative', 'ttf_growth_rate']].head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:52.09711Z","iopub.execute_input":"2025-12-04T04:44:52.097505Z","iopub.status.idle":"2025-12-04T04:44:52.142431Z","shell.execute_reply.started":"2025-12-04T04:44:52.097393Z","shell.execute_reply":"2025-12-04T04:44:52.140658Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# EP_disp ve türevini aynı grafikte çizme\nplt.figure(figsize=(12, 6))\nplt.plot(df.index, df['EP_disp'], label='EP_disp', linewidth=2)\nplt.plot(df.index, df['EP_disp_derivative'], label='EP_disp Derivative', linestyle='--')\nplt.title('EP_disp and its Derivative')\nplt.xlabel('Index')\nplt.ylabel('Value')\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:52.144523Z","iopub.execute_input":"2025-12-04T04:44:52.144925Z","iopub.status.idle":"2025-12-04T04:44:52.778635Z","shell.execute_reply.started":"2025-12-04T04:44:52.144888Z","shell.execute_reply":"2025-12-04T04:44:52.777393Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ttf ve türevini aynı grafikte çizme\nplt.figure(figsize=(12, 6))\nplt.plot(df.index, df['ttf'], label='ttf', linewidth=2)\nplt.plot(df.index, df['ttf_derivative'], label='ttf Derivative', linestyle='--')\nplt.title('TTF and its Derivative')\nplt.xlabel('Index')\nplt.ylabel('Value')\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:52.780216Z","iopub.execute_input":"2025-12-04T04:44:52.780663Z","iopub.status.idle":"2025-12-04T04:44:53.375404Z","shell.execute_reply.started":"2025-12-04T04:44:52.780612Z","shell.execute_reply":"2025-12-04T04:44:53.374011Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(df[['EP_disp', 'ttf']].describe())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T04:44:53.377018Z","iopub.execute_input":"2025-12-04T04:44:53.377395Z","iopub.status.idle":"2025-12-04T04:44:53.404967Z","shell.execute_reply.started":"2025-12-04T04:44:53.377357Z","shell.execute_reply":"2025-12-04T04:44:53.403652Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 完整的地震预测竞赛解决方案\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.preprocessing import StandardScaler, RobustScaler\nfrom sklearn.model_selection import TimeSeriesSplit\nfrom sklearn.metrics import mean_absolute_error\nimport tensorflow as tf\nfrom tensorflow.keras import layers, Model, callbacks\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# 设置随机种子以确保可重复性\ntf.random.set_seed(42)\nnp.random.seed(42)\n\nprint(\"=== 地震预测竞赛完整解决方案 ===\")\n\n# 1. 数据加载和基础EDA（基于你的代码）\ndf = pd.read_csv(\"/kaggle/input/earthquake-prediction/train.csv\")\nprint(f\"训练集形状: {df.shape}\")\nprint(f\"列名: {df.columns.tolist()}\")\n\n# 基础EDA（你的代码）\nprint(\"\\n=== 基础数据探索 ===\")\nprint(df.isnull().sum())\nprint(df[['ttf', 'EP_disp', 'EP_disp_FUT_1', 'EP_disp_FUT_2', 'EP_disp_FUT_3']].describe())\n\n# 2. 高级特征工程\nclass AdvancedFeatureEngineer:\n    def __init__(self):\n        self.scaler = RobustScaler()\n        \n    def extract_features(self, df):\n        \"\"\"从6000个AE信号中提取高级特征\"\"\"\n        features_list = []\n        \n        # 对每个样本的6000个AE读数进行特征提取\n        for idx in range(len(df)):\n            ae_data = df.iloc[idx, 1:6001].values  # ae_1到ae_6000\n            \n            # 时域特征\n            time_features = self._extract_time_domain_features(ae_data)\n            \n            # 频域特征\n            freq_features = self._extract_frequency_domain_features(ae_data)\n            \n            # 非线性特征\n            nonlinear_features = self._extract_nonlinear_features(ae_data)\n            \n            # 统计特征\n            stat_features = self._extract_statistical_features(ae_data)\n            \n            # 组合所有特征\n            all_features = np.concatenate([\n                time_features, freq_features, nonlinear_features, stat_features\n            ])\n            \n            features_list.append(all_features)\n        \n        return np.array(features_list)\n    \n    def _extract_time_domain_features(self, signal):\n        \"\"\"提取时域特征\"\"\"\n        features = []\n        # RMS\n        features.append(np.sqrt(np.mean(signal**2)))\n        # 峰值\n        features.append(np.max(np.abs(signal)))\n        # 峰峰值\n        features.append(np.max(signal) - np.min(signal))\n        # 波形因子\n        features.append(np.sqrt(np.mean(signal**2)) / np.mean(np.abs(signal)))\n        # 脉冲因子\n        features.append(np.max(np.abs(signal)) / np.mean(np.abs(signal)))\n        # 裕度因子\n        features.append(np.max(np.abs(signal)) / (np.mean(np.sqrt(np.abs(signal))))**2)\n        \n        return np.array(features)\n    \n    def _extract_frequency_domain_features(self, signal):\n        \"\"\"提取频域特征\"\"\"\n        fft_vals = np.fft.fft(signal)\n        fft_abs = np.abs(fft_vals[:len(fft_vals)//2])\n        \n        features = []\n        if len(fft_abs) > 0:\n            # 频谱重心\n            features.append(np.sum(fft_abs * np.arange(len(fft_abs))) / np.sum(fft_abs))\n            # 频谱宽度\n            freq_centroid = np.sum(fft_abs * np.arange(len(fft_abs))) / np.sum(fft_abs)\n            features.append(np.sqrt(np.sum(fft_abs * (np.arange(len(fft_abs)) - freq_centroid)**2) / np.sum(fft_abs)))\n            # 主频\n            features.append(np.argmax(fft_abs))\n        else:\n            features.extend([0, 0, 0])\n            \n        return np.array(features)\n    \n    def _extract_nonlinear_features(self, signal):\n        \"\"\"提取非线性特征\"\"\"\n        # 近似熵（简化版）\n        def approximate_entropy(signal, m=2, r=0.2):\n            if len(signal) < m + 1:\n                return 0\n            # 简化计算\n            return np.std(signal)\n        \n        features = [\n            approximate_entropy(signal),\n            np.mean(np.diff(signal)**2)  # 变化率平方均值\n        ]\n        return np.array(features)\n    \n    def _extract_statistical_features(self, signal):\n        \"\"\"提取统计特征\"\"\"\n        features = [\n            np.mean(signal),\n            np.std(signal),\n            np.median(signal),\n            np.min(signal),\n            np.max(signal),\n            np.percentile(signal, 25),\n            np.percentile(signal, 75),\n            np.mean(np.abs(signal)),  # 绝对均值\n            signal[-1] - signal[0]   # 首尾差值\n        ]\n        return np.array(features)\n\n# 3. 构建高效的多任务深度学习模型\nclass EarthquakeMultiTaskModel(Model):\n    def __init__(self, input_dim, num_features=19):\n        super(EarthquakeMultiTaskModel, self).__init__()\n        self.input_dim = input_dim\n        self.num_features = num_features\n        \n        # 特征编码器\n        self.feature_encoder = tf.keras.Sequential([\n            layers.Dense(512, activation='relu'),\n            layers.BatchNormalization(),\n            layers.Dropout(0.3),\n            layers.Dense(256, activation='relu'),\n            layers.BatchNormalization(),\n            layers.Dropout(0.3),\n            layers.Dense(128, activation='relu')\n        ])\n        \n        # 注意力机制\n        self.attention = layers.MultiHeadAttention(num_heads=4, key_dim=32)\n        self.attention_dense = layers.Dense(128, activation='relu')\n        \n        # 共享特征层\n        self.shared_dense1 = layers.Dense(256, activation='relu')\n        self.shared_dense2 = layers.Dense(128, activation='relu')\n        \n        # 实时预测分支\n        self.real_time_branch = tf.keras.Sequential([\n            layers.Dense(64, activation='relu'),\n            layers.Dense(32, activation='relu'),\n            layers.Dense(2, activation='linear')  # ttf + EP_disp\n        ])\n        \n        # 未来预测分支\n        self.future_branch = tf.keras.Sequential([\n            layers.Dense(64, activation='relu'),\n            layers.Dense(32, activation='relu'),\n            layers.Dense(3, activation='linear')  # EP_disp_FUT_1,2,3\n        ])\n        \n    def call(self, inputs):\n        # 特征编码\n        encoded = self.feature_encoder(inputs)\n        \n        # 自注意力（添加序列维度）\n        encoded_expanded = tf.expand_dims(encoded, axis=1)\n        attended = self.attention(encoded_expanded, encoded_expanded)\n        attended = tf.squeeze(attended, axis=1)\n        attended = self.attention_dense(attended)\n        \n        # 合并特征\n        combined = layers.concatenate([encoded, attended])\n        \n        # 共享特征学习\n        shared = self.shared_dense1(combined)\n        shared = self.shared_dense2(shared)\n        \n        # 多任务输出\n        real_time_output = self.real_time_branch(shared)\n        future_output = self.future_branch(shared)\n        \n        # 合并所有输出\n        outputs = tf.concat([real_time_output, future_output], axis=1)\n        return outputs\n\n# 4. 数据预处理和特征工程\nprint(\"\\n=== 特征工程开始 ===\")\nfeature_engineer = AdvancedFeatureEngineer()\n\n# 提取高级特征\nprint(\"正在提取高级特征...\")\nX_features = feature_engineer.extract_features(df)\n\n# 目标变量\ny_columns = ['ttf', 'EP_disp', 'EP_disp_FUT_1', 'EP_disp_FUT_2', 'EP_disp_FUT_3']\ny = df[y_columns].values\n\nprint(f\"特征矩阵形状: {X_features.shape}\")\nprint(f\"目标变量形状: {y.shape}\")\n\n# 数据标准化\nscaler_X = StandardScaler()\nX_scaled = scaler_X.fit_transform(X_features)\n\nscaler_y = StandardScaler()\ny_scaled = scaler_y.fit_transform(y)\n\n# 时间序列分割\ntscv = TimeSeriesSplit(n_splits=5)\ntrain_val_splits = list(tscv.split(X_scaled))\n\n# 5. 模型训练配置\ndef create_model():\n    model = EarthquakeMultiTaskModel(input_dim=X_scaled.shape[1])\n    \n    # 自定义损失函数（加权MAE）\n    def weighted_mae(y_true, y_pred):\n        # 为不同目标分配不同权重\n        weights = tf.constant([1.0, 0.8, 0.6, 0.4, 0.2])  # ttf权重最高\n        mae = tf.reduce_mean(tf.abs(y_true - y_pred) * weights, axis=-1)\n        return mae\n    \n    model.compile(\n        optimizer=tf.keras.optimizers.AdamW(learning_rate=0.001, weight_decay=0.0001),\n        loss=weighted_mae,\n        metrics=['mae']\n    )\n    return model\n\n# 6. 训练函数\ndef train_model_with_cv(X, y, n_splits=3):\n    \"\"\"使用时间序列交叉验证训练模型\"\"\"\n    tscv = TimeSeriesSplit(n_splits=n_splits)\n    models = []\n    scores = []\n    \n    for fold, (train_idx, val_idx) in enumerate(tscv.split(X)):\n        print(f\"\\n=== 训练第 {fold+1} 折 ===\")\n        \n        X_train, X_val = X[train_idx], X[val_idx]\n        y_train, y_val = y[train_idx], y[val_idx]\n        \n        model = create_model()\n        \n        # 回调函数\n        callbacks_list = [\n            callbacks.EarlyStopping(patience=15, restore_best_weights=True),\n            callbacks.ReduceLROnPlateau(factor=0.5, patience=10),\n            callbacks.ModelCheckpoint(f'best_model_fold{fold}.keras', save_best_only=True)\n        ]\n        \n        # 训练\n        history = model.fit(\n            X_train, y_train,\n            validation_data=(X_val, y_val),\n            epochs=100,\n            batch_size=64,\n            callbacks=callbacks_list,\n            verbose=1\n        )\n        \n        # 评估\n        val_score = model.evaluate(X_val, y_val, verbose=0)\n        scores.append(val_score[0])\n        models.append(model)\n        \n        print(f\"第 {fold+1} 折验证MAE: {val_score[0]:.4f}\")\n    \n    return models, scores\n\n# 7. 开始训练\nprint(\"\\n=== 开始模型训练 ===\")\nmodels, scores = train_model_with_cv(X_scaled, y_scaled, n_splits=3)\n\nprint(f\"\\n平均交叉验证MAE: {np.mean(scores):.4f} ± {np.std(scores):.4f}\")\n\n# 8. 模型集成预测\ndef ensemble_predict(models, X):\n    \"\"\"使用模型集成进行预测\"\"\"\n    predictions = []\n    for model in models:\n        pred = model.predict(X, verbose=0)\n        predictions.append(pred)\n    \n    # 平均预测结果\n    ensemble_pred = np.mean(predictions, axis=0)\n    return ensemble_pred\n\n# 9. 加载测试数据并预测\nprint(\"\\n=== 加载测试数据 ===\")\ntest_df = pd.read_csv(\"/kaggle/input/earthquake-prediction/test.csv\")\nprint(f\"测试集形状: {test_df.shape}\")\n\n# 为测试集提取特征\nprint(\"为测试集提取特征...\")\nX_test_features = feature_engineer.extract_features(test_df)\nX_test_scaled = scaler_X.transform(X_test_features)\n\n# 集成预测\nprint(\"进行集成预测...\")\ny_test_pred_scaled = ensemble_predict(models, X_test_scaled)\n\n# 反标准化\ny_test_pred = scaler_y.inverse_transform(y_test_pred_scaled)\n\n# 10. 生成提交文件\nprint(\"\\n=== 生成提交文件 ===\")\nsubmission = pd.DataFrame({\n    'index': range(1, len(test_df) + 1),\n    'time_stamp': test_df['time_stamp'],\n    'ttf': y_test_pred[:, 0],\n    'EP_disp': y_test_pred[:, 1],\n    'EP_disp_FUT_1': y_test_pred[:, 2],\n    'EP_disp_FUT_2': y_test_pred[:, 3],\n    'EP_disp_FUT_3': y_test_pred[:, 4]\n})\n\n# 确保预测值非负（物理约束）\nsubmission['ttf'] = np.maximum(submission['ttf'], 0)\nfor col in ['EP_disp', 'EP_disp_FUT_1', 'EP_disp_FUT_2', 'EP_disp_FUT_3']:\n    submission[col] = np.maximum(submission[col], 0)\n\n# 保存提交文件\nsubmission.to_csv('submission.csv', index=False)\nprint(\"提交文件已保存为 'submission.csv'\")\n\n# 11. 结果可视化\nprint(\"\\n=== 预测结果统计 ===\")\nprint(submission[['ttf', 'EP_disp', 'EP_disp_FUT_1', 'EP_disp_FUT_2', 'EP_disp_FUT_3']].describe())\n\n# 绘制预测分布\nplt.figure(figsize=(15, 10))\nfor i, col in enumerate(['ttf', 'EP_disp', 'EP_disp_FUT_1', 'EP_disp_FUT_2', 'EP_disp_FUT_3'], 1):\n    plt.subplot(2, 3, i)\n    plt.hist(submission[col], bins=50, alpha=0.7, color='skyblue')\n    plt.title(f'{col} Distribution')\n    plt.xlabel(col)\n    plt.ylabel('Frequency')\n\nplt.tight_layout()\nplt.savefig('prediction_distributions.png')\nplt.show()\n\nprint(\"\\n=== 解决方案完成 ===\")\nprint(\"主要创新点:\")\nprint(\"1. 高级特征工程: 时域+频域+非线性+统计特征\")\nprint(\"2. 多任务学习: 同时预测TTF和多个位移值\")\nprint(\"3. 注意力机制: 增强模型对关键特征的关注\")\nprint(\"4. 时间序列交叉验证: 确保模型泛化能力\")\nprint(\"5. 模型集成: 提升预测稳定性\")\nprint(\"6. 物理约束: 确保预测值符合物理现实\")\n\n# 12. 模型性能分析（可选）\ndef analyze_feature_importance(models, feature_names):\n    \"\"\"分析特征重要性（简化版）\"\"\"\n    if hasattr(models[0].layers[0], 'get_weights'):\n        weights = models[0].layers[0].get_weights()[0]\n        importance = np.mean(np.abs(weights), axis=1)\n        \n        # 创建特征重要性DataFrame\n        feature_importance = pd.DataFrame({\n            'feature': feature_names,\n            'importance': importance\n        }).sort_values('importance', ascending=False)\n        \n        return feature_importance.head(10)\n    return None\n\n# 特征名称（根据实际特征工程调整）\nfeature_names = [\n    'RMS', 'Peak', 'Peak-to-Peak', 'Waveform Factor', 'Impulse Factor', 'Margin Factor',\n    'Spectral Centroid', 'Spectral Width', 'Dominant Frequency',\n    'Approx Entropy', 'Change Rate',\n    'Mean', 'Std', 'Median', 'Min', 'Max', 'Q25', 'Q75', 'Abs Mean', 'Delta'\n]\n\nimportance_df = analyze_feature_importance(models, feature_names)\nif importance_df is not None:\n    print(\"\\nTop 10 重要特征:\")\n    print(importance_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-04T06:35:58.622267Z","iopub.execute_input":"2025-12-04T06:35:58.622699Z","iopub.status.idle":"2025-12-04T07:50:37.513705Z","shell.execute_reply.started":"2025-12-04T06:35:58.622659Z","shell.execute_reply":"2025-12-04T07:50:37.511605Z"}},"outputs":[],"execution_count":null}]}