import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.sandbox.stats.runs import runstest_1samp
import seaborn as sns
import numpy as np
from matplotlib import rcParams
from scipy.stats import skew, kurtosis
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error
from statsmodels.graphics.gofplots import qqplot

# 配置 Matplotlib 中文字体支持，避免中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 替换为本地的中文字体，如 "SimSun" 或 "SimHei"
plt.rcParams['axes.unicode_minus'] = False   # 避免负号显示问题

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.sandbox.stats.runs import runstest_1samp
import seaborn as sns
import numpy as np
from matplotlib import rcParams

# 配置 Matplotlib 中文字体支持，避免中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 替换为本地的中文字体，如 "SimSun" 或 "SimHei"
plt.rcParams['axes.unicode_minus'] = False   # 避免负号显示问题

# 1. 数据导入
file_path = 'flu_positive.xlsx'  
sheet_name = 'Positive Rate'  

# 加载数据
df = pd.read_excel(file_path, sheet_name=sheet_name)

# 检查数据
print("数据预览:")
display(df.head())

# 确保目标列存在，提取目标列数据
if 'ILI+' not in df.columns:
    raise ValueError("目标列 'ILI+' 不存在，请检查数据文件！")

# 提取目标列数据
data_values = df['ILI+'].dropna().values  # 删除缺失值，确保数据完整



# 绘制时间序列图
plt.figure(figsize=(12, 6))
plt.plot(data_values,alpha=0.7)
plt.title("时间序列数据可视化")
plt.xlabel("时间步")
plt.ylabel("ILI+ 值")
plt.grid()
plt.show()


#  平稳性检验（ADF检验）
adf_result = adfuller(data_values)
print("\nADF（Augmented Dickey-Fuller）单位根检验结果:")
print(f"ADF 统计量: {adf_result[0]:.3f}")
print(f"p 值: {adf_result[1]:.3f}")
print("临界值:")
for key, value in adf_result[4].items():
    print(f"  {key}: {value:.3f}")
if adf_result[1] < 0.05:
    print("结论: 数据是平稳的，可以直接用于时间序列建模。")
else:
    print("结论: 数据是非平稳的，建议进行差分或对数变换。")

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import MinMaxScaler
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error
from statsmodels.graphics.gofplots import qqplot
from statsmodels.graphics.tsaplots import plot_acf

# 数据导入
file_path = 'flu_positive.xlsx'
data = pd.ExcelFile(file_path)
df = data.parse('Positive Rate')

# 提取目标列 "ILI+" 并标准化数据
scaler = MinMaxScaler()
df['num_scaled'] = scaler.fit_transform(df[['ILI+']])

# 创建时间序列数据 (设定序列长度)
sequence_length = 30

def create_sequences(data, sequence_length):
    X, y = [], []
    for i in range(sequence_length, len(data)):
        X.append(data[i - sequence_length:i])
        y.append(data[i])
    return np.array(X), np.array(y)

# 准备序列数据
data_scaled = df['num_scaled'].values
X, y = create_sequences(data_scaled, sequence_length)

# 创建对应的日期
dates = df['date'].iloc[sequence_length:].reset_index(drop=True)

# 检查 GPU 是否可用
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(f"Using device: {device}")

# 定义 GRU 模型
class GRUModel(nn.Module):
    def __init__(self, input_size, hidden_size, output_size, num_layers):
        super(GRUModel, self).__init__()
        self.gru = nn.GRU(input_size, hidden_size, num_layers, batch_first=True)
        self.fc = nn.Linear(hidden_size, output_size)
        self.activation = nn.Softplus()  # 添加 Softplus 激活函数

    def forward(self, x):
        h0 = torch.zeros(self.gru.num_layers, x.size(0), self.gru.hidden_size).to(device)
        out, _ = self.gru(x, h0)
        out = self.fc(out[:, -1, :])  # 取最后一个时间步的输出
        out = self.activation(out)  # 确保输出非负
        return out

# 自定义损失函数，加入对负值的惩罚
def custom_loss_function(outputs, targets):
    mse_loss = nn.MSELoss()(outputs, targets)
    penalty = torch.mean(torch.relu(-outputs))  # 对负值部分进行惩罚
    return mse_loss + penalty * 10  # 10 是惩罚系数，可以调整

# 定义 GRU 模型参数
input_size = 1
hidden_size = 128
output_size = 1
num_layers = 2

# 定义 Bootstrap 训练函数
def bootstrap_train_model(base_model, train_loader, criterion, optimizer_fn, n_bootstrap=100, epochs_per_bootstrap=50, noise_scale=0.4):
    all_results = {
        "train_predictions": [],
        "test_predictions": [],
        "future_predictions": []
    }
    for i in range(n_bootstrap):
        print(f"Bootstrap iteration {i+1}/{n_bootstrap}")
        # 初始化模型和优化器
        model = GRUModel(input_size, hidden_size, output_size, num_layers).to(device)
        optimizer = optimizer_fn(model.parameters())

        model.train()
        for epoch in range(epochs_per_bootstrap):
            epoch_loss = 0
            # 动态调整学习率
            if epoch == 350:
                for param_group in optimizer.param_groups:
                    param_group['lr'] = 0.001

            for inputs, targets in train_loader:
                inputs = inputs.unsqueeze(-1).to(device)
                targets = targets.unsqueeze(-1).to(device)

                # 添加噪声到输入
                noisy_inputs = inputs + torch.randn_like(inputs) * noise_scale

                # 前向传播
                outputs = model(noisy_inputs)
                loss = criterion(outputs, targets)

                # 反向传播和优化
                optimizer.zero_grad()
                loss.backward()
                optimizer.step()

                epoch_loss += loss.item()

            if (epoch + 1) % 100 == 0:
                print(f"Epoch [{epoch + 1}/{epochs_per_bootstrap}], Loss: {epoch_loss / len(train_loader):.4f}")

        # 保存训练集预测
        model.eval()
        with torch.no_grad():
            train_preds = model(torch.tensor(X_train, dtype=torch.float32).unsqueeze(-1).to(device)).cpu().numpy()
            test_preds = model(torch.tensor(X_test, dtype=torch.float32).unsqueeze(-1).to(device)).cpu().numpy()

            # 获取未来预测
            future_preds = []
            current_sequence = data_scaled[-sequence_length:]
            for _ in range(21):
                input_tensor = torch.tensor(current_sequence, dtype=torch.float32).unsqueeze(0).unsqueeze(-1).to(device)
                next_pred = model(input_tensor).cpu().numpy().flatten()
                future_preds.append(next_pred[0])
                current_sequence = np.append(current_sequence[1:], next_pred)

        all_results["train_predictions"].append(train_preds)
        all_results["test_predictions"].append(test_preds)
        all_results["future_predictions"].append(future_preds)

    return all_results

# 计算均值和置信区间
def calculate_mean_and_ci(results):
    results = np.array(results)  # 转为 NumPy 数组 (n_bootstrap, n_samples)
    mean_preds = results.mean(axis=0)
    lower_bound = np.percentile(results, 2.5, axis=0)
    upper_bound = np.percentile(results, 97.5, axis=0)
    return mean_preds, lower_bound, upper_bound

# 创建存储指标的 DataFrame
metrics_df = pd.DataFrame(columns=['Fold', 'MAE', 'RMSE', 'MAPE'])

# 创建一个列表来存储所有残差、预测值和真实值
all_residuals = []
all_predictions = []
all_true_values = []

# 定义交叉验证 (3 折)
ts_split = TimeSeriesSplit(n_splits=4)
fold = 1
for train_index, test_index in ts_split.split(X):
    print(f"Fold {fold}...")

    # 划分训练集和测试集
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]

    # 转换为 PyTorch 张量
    X_train_tensor = torch.tensor(X_train, dtype=torch.float32)
    X_test_tensor = torch.tensor(X_test, dtype=torch.float32)
    y_train_tensor = torch.tensor(y_train, dtype=torch.float32)
    y_test_tensor = torch.tensor(y_test, dtype=torch.float32)

    # 创建数据加载器
    train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
    test_dataset = TensorDataset(X_test_tensor, y_test_tensor)
    train_loader = DataLoader(train_dataset, batch_size=32, shuffle=False)
    test_loader = DataLoader(test_dataset, batch_size=32, shuffle=False)

    # 训练 Bootstrap 模型
    bootstrap_results = bootstrap_train_model(
        GRUModel, train_loader, custom_loss_function, 
        optimizer_fn=lambda params: optim.Adam(params, lr=0.005), 
        n_bootstrap=10, epochs_per_bootstrap=500, noise_scale=0.2
    )

    # 计算训练集、测试集、未来预测的均值和置信区间
    train_mean = np.mean(bootstrap_results["train_predictions"], axis=0)
    test_mean, test_lower, test_upper = calculate_mean_and_ci(bootstrap_results["test_predictions"])
    future_mean, future_lower, future_upper = calculate_mean_and_ci(bootstrap_results["future_predictions"])

    # 将训练集和测试集预测值反归一化
    train_mean_rescaled = scaler.inverse_transform(train_mean.reshape(-1, 1))
    test_mean_rescaled = scaler.inverse_transform(test_mean.reshape(-1, 1))
    test_lower_rescaled = scaler.inverse_transform(test_lower.reshape(-1, 1))
    test_upper_rescaled = scaler.inverse_transform(test_upper.reshape(-1, 1))
    future_mean_rescaled = scaler.inverse_transform(future_mean.reshape(-1, 1))
    future_lower_rescaled = scaler.inverse_transform(future_lower.reshape(-1, 1))
    future_upper_rescaled = scaler.inverse_transform(future_upper.reshape(-1, 1))

# 计算误差指标
    epsilon = 1e-6  # 防止分母为零的平滑常数
    y_test_rescaled = scaler.inverse_transform(y_test.reshape(-1, 1))
    mae = mean_absolute_error(y_test_rescaled, test_mean_rescaled)
    rmse = mean_squared_error(y_test_rescaled, test_mean_rescaled, squared=False)
    mape = np.mean(np.abs((y_test_rescaled - test_mean_rescaled) / (y_test_rescaled + epsilon))) * 100
    smape = np.mean(2 * np.abs(y_test_rescaled - test_mean_rescaled) / 
                (np.abs(y_test_rescaled) + np.abs(test_mean_rescaled) + epsilon)) * 100

    print(f"Fold {fold} Metrics:")
    print(f"MAE: {mae:.4f}")
    print(f"RMSE: {rmse:.4f}")
    print(f"MAPE: {mape:.2f}%")
    print(f"SMAPE: {smape:.2f}%")

# 将指标写入 DataFrame
    new_row = pd.DataFrame({
        'Fold': [fold],
        'MAE': [mae],
        'RMSE': [rmse],
        'MAPE': [mape],
        'SMAPE': [smape]
    })
    metrics_df = pd.concat([metrics_df, new_row], ignore_index=True)

    # 计算残差
    residuals = y_test_rescaled.flatten() - test_mean_rescaled.flatten()

    # 将残差添加到 all_residuals 列表中
    all_residuals.extend(residuals)

    # 存储预测值和真实值
    all_predictions.extend(test_mean_rescaled.flatten())
    all_true_values.extend(y_test_rescaled.flatten())

    # 创建未来日期范围
    last_date = pd.to_datetime(dates.iloc[-1])
    future_dates = [last_date + pd.Timedelta(days=i+1) for i in range(21)]

    # 可视化结果 (包含训练集、测试集和未来预测)
    plt.figure(figsize=(15, 8))

    # 绘制所有真实值
    plt.plot(df['date'], df['ILI+'], label='True Values (All)', color='blue')

    # 绘制训练集预测值
    plt.plot(dates.iloc[train_index], train_mean_rescaled, label='Train Predictions', color='green')

    # 绘制测试集预测值及置信区间
    plt.plot(dates.iloc[test_index], test_mean_rescaled, label='Test Predictions', color='red')
    plt.fill_between(dates.iloc[test_index], test_lower_rescaled.flatten(), test_upper_rescaled.flatten(), color='gray', alpha=0.3, label='95% Confidence Interval (Test)')

    # 绘制未来预测值及置信区间
    plt.plot(future_dates, future_mean_rescaled, label='Future Predictions (Next 21 Days)', color='orange', linestyle='--')
    plt.fill_between(future_dates, future_lower_rescaled.flatten(), future_upper_rescaled.flatten(), color='yellow', alpha=0.3, label='95% Confidence Interval (Future)')

    plt.title(f'Bootstrap GRU Model Predictions for Fold {fold}')
    plt.xlabel('Date')
    plt.ylabel('ILI+')
    plt.legend()
    plt.xticks(rotation=45)
    plt.tight_layout()
    plt.show()

    fold += 1

# 可视化各折叠的 MAE、RMSE、MAPE 和 SMAPE
fig, axes = plt.subplots(1, 4, figsize=(40, 5))

# MAE 可视化
axes[0].plot(metrics_df['Fold'], metrics_df['MAE'], marker='o', label='MAE')
axes[0].set_title('MAE Across Folds')
axes[0].set_xlabel('Fold')
axes[0].set_ylabel('MAE')
axes[0].legend()

# RMSE 可视化
axes[1].plot(metrics_df['Fold'], metrics_df['RMSE'], marker='o', label='RMSE', color='orange')
axes[1].set_title('RMSE Across Folds')
axes[1].set_xlabel('Fold')
axes[1].set_ylabel('RMSE')
axes[1].legend()

# MAPE 可视化
axes[2].plot(metrics_df['Fold'], metrics_df['MAPE'], marker='o', label='MAPE', color='green')
axes[2].set_title('MAPE Across Folds')
axes[2].set_xlabel('Fold')
axes[2].set_ylabel('MAPE (%)')
axes[2].legend()

# SMAPE 可视化
axes[3].plot(metrics_df['Fold'], metrics_df['SMAPE'], marker='o', label='SMAPE', color='purple')
axes[3].set_title('SMAPE Across Folds')
axes[3].set_xlabel('Fold')
axes[3].set_ylabel('SMAPE (%)')
axes[3].legend()

plt.tight_layout()
plt.show()

# 残差分析和可视化
all_residuals = np.array(all_residuals)
all_predictions = np.array(all_predictions)
all_true_values = np.array(all_true_values)

# 计算基本统计量
residual_mean = np.mean(all_residuals)
residual_std = np.std(all_residuals)
print(f"Residual Mean: {residual_mean:.4f}")
print(f"Residual Std Dev: {residual_std:.4f}")

# 残差分布图
plt.figure(figsize=(10, 6))
sns.histplot(all_residuals, bins=50, kde=True, color='purple')
plt.title('Residuals Distribution')
plt.xlabel('Residual')
plt.ylabel('Frequency')
plt.show()



# Q-Q 图
plt.figure(figsize=(10, 6))
qqplot(all_residuals, line='s')
plt.title('Q-Q Plot of Residuals')
plt.show()

# 自相关函数 (ACF) 图
plt.figure(figsize=(10, 6))
plot_acf(all_residuals, lags=100)
plt.title('Autocorrelation of Residuals')
plt.show()

# 残差时间序列图
plt.figure(figsize=(15, 6))
plt.plot(all_true_values, label='True Values', color='blue')
plt.plot(all_predictions, label='Predicted Values', color='red')
plt.scatter(range(len(all_residuals)), all_residuals, label='Residuals', color='green', alpha=0.5)
plt.axhline(0, color='black', linestyle='--')
plt.title('Residuals Over Time')
plt.xlabel('Time')
plt.ylabel('Values')
plt.legend()
plt.show()