stacking-1004.py

import os
import time
import warnings
import traceback
import logging

import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, clone
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.model_selection import TimeSeriesSplit
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import MinMaxScaler
from xgboost import XGBRegressor
from typing import Dict

warnings.filterwarnings('ignore')



def calculate_rmse(y_true, y_pred):
    return np.sqrt(mean_squared_error(y_true, y_pred))


def calculate_nse(y_true, y_pred):
    mean_y_true = np.mean(y_true)
    numerator = np.sum((y_true - y_pred) ** 2)
    denominator = np.sum((y_true - mean_y_true) ** 2)
    return 1 - (numerator / denominator) if denominator != 0 else float('inf')


def calculate_r2(y_true, y_pred):
    mean_y_true = np.mean(y_true)
    mean_y_pred = np.mean(y_pred)
    numerator = np.sum((y_true - mean_y_true) * (y_pred - mean_y_pred)) ** 2
    denominator = np.sum((y_true - mean_y_true) ** 2) * np.sum((y_pred - mean_y_pred) ** 2)
    return numerator / denominator if denominator != 0 else float('inf')


def calculate_bias(y_true, y_pred):
    return np.mean(y_pred - y_true)


def blocked_bootstrap_ci(y_true, y_pred, metric_func, n_bootstrap=1000,
                         alpha=0.05, block_size=60, seed=42):
    rng = np.random.default_rng(seed)
    stats = []
    n = len(y_true)
    y_true, y_pred = np.asarray(y_true), np.asarray(y_pred)

    if n < block_size:
        return np.nan, np.nan, np.nan

    num_blocks = n - block_size + 1

    for _ in range(n_bootstrap):
        num_samples = (n + block_size - 1) // block_size
        start_indices = rng.choice(num_blocks, size=num_samples, replace=True)

        idx = [i for start in start_indices for i in range(start, start + block_size)]
        idx = idx[:n]

        if len(idx) != n:
            continue

        stat_val = metric_func(y_true[idx], y_pred[idx])
        if np.isfinite(stat_val):
            stats.append(stat_val)

    if len(stats) < 100:
        return np.nan, np.nan, np.nan

    lower = np.percentile(stats, 100 * alpha / 2)
    upper = np.percentile(stats, 100 * (1 - alpha / 2))
    width = upper - lower
    return lower, upper, width


def assign_season(month):
    if month in [12, 1, 2]: return "winter"
    if month in [3, 4, 5]: return "spring"
    if month in [6, 7, 8]: return "summer"
    return "autumn"

class TimeSeriesStacking:
    def __init__(self, base_models: Dict[str, BaseEstimator], meta_model: BaseEstimator, cv_splits=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.cv_splits = cv_splits
        self.fitted_base_models = {}
        self.fitted_meta_model = None

    def fit(self, X, y):
        tscv = TimeSeriesSplit(n_splits=self.cv_splits)
        meta_features_list, meta_target_list = [], []

        for train_idx, val_idx in tscv.split(X):
            X_train, X_val = X[train_idx], X[val_idx]
            y_train, y_val = y[train_idx], y[val_idx]

            fold_meta_features = []
            for name, model in self.base_models.items():
                model_clone = clone(model)
                model_clone.fit(X_train, y_train)
                y_pred_val = model_clone.predict(X_val)
                fold_meta_features.append(y_pred_val)

            meta_features_list.append(np.column_stack(fold_meta_features))
            meta_target_list.append(y_val)

        meta_features = np.vstack(meta_features_list)
        meta_target = np.concatenate(meta_target_list)
        self.fitted_meta_model = clone(self.meta_model).fit(meta_features, meta_target)

        for name, model in self.base_models.items():
            self.fitted_base_models[name] = clone(model).fit(X, y)

        return self

    def predict(self, X):
        base_preds = [model.predict(X) for model in self.fitted_base_models.values()]
        meta_features = np.column_stack(base_preds)
        return self.fitted_meta_model.predict(meta_features)


def evaluate_models_for_site(file_path, site_hyperparams_df, features, input_set_name, seed):
    file_name = os.path.splitext(os.path.basename(file_path))[0]
    logging.info(f"    -> Processing site: {file_name:15s} | Input: {input_set_name}")

    data = pd.read_excel(file_path)
    data['date'] = pd.to_datetime(data['date'])
    data['year'] = data['date'].dt.year
    data['month'] = data['date'].dt.month
    data['season'] = data['month'].apply(assign_season)

    years = sorted(data['year'].unique())
    split_idx = int(len(years) * 0.8)
    train_years, test_years = years[:split_idx], years[split_idx:]
    train_data = data[data['year'].isin(train_years)]
    test_data = data[data['year'].isin(test_years)]

    X_train_df, y_train = train_data[features], train_data['ET'].values
    X_test_df, y_test = test_data[features], test_data['ET'].values

    scaler = MinMaxScaler()
    X_train = scaler.fit_transform(X_train_df)
    X_test = scaler.transform(X_test_df)

    try:
        rf_hyperparams = site_hyperparams_df[site_hyperparams_df['模型'] == 'RF'].iloc[0]
        xgb_hyperparams = site_hyperparams_df[site_hyperparams_df['模型'] == 'XGB'].iloc[0]
        mlp_hyperparams = site_hyperparams_df[site_hyperparams_df['模型'] == 'MLP'].iloc[0]
    except IndexError:
        logging.error(f"      [ERROR] Missing one or more model hyperparameters for site {file_name}. Skipping.")
        return pd.DataFrame(), pd.DataFrame()

    rf_params = {'n_estimators': int(rf_hyperparams['n_estimators']), 'max_depth': int(rf_hyperparams['max_depth']),
                 'min_samples_split': int(rf_hyperparams['min_samples_split']),
                 'min_samples_leaf': int(rf_hyperparams['min_samples_leaf']),
                 'random_state': seed}

    xgb_params = {'n_estimators': int(xgb_hyperparams['n_estimators']), 'max_depth': int(xgb_hyperparams['max_depth']),
                  'learning_rate': xgb_hyperparams['learning_rate'], 'subsample': xgb_hyperparams['subsample'],
                  'colsample_bytree': xgb_hyperparams['colsample_bytree'],
                  'random_state': seed, 'seed': seed}

    mlp_params = {'hidden_layer_sizes': eval(str(mlp_hyperparams['hidden_layer_sizes'])),
                  'learning_rate_init': mlp_hyperparams['learning_rate_init'],
                  'alpha': mlp_hyperparams['alpha'], 'random_state': seed, 'max_iter': 1000}

    base_models = {'RF': RandomForestRegressor(**rf_params), 'XGB': XGBRegressor(**xgb_params),
                   'MLP': MLPRegressor(**mlp_params)}

    predictions = {name: model.fit(X_train, y_train).predict(X_test) for name, model in base_models.items()}

    ts_stacking = TimeSeriesStacking(base_models=base_models, meta_model=LinearRegression(), cv_splits=5)
    ts_stacking.fit(X_train, y_train)
    predictions['Stacking'] = ts_stacking.predict(X_test)

    site_results = []
    test_seasons = test_data['season'].values

    for model_name, y_pred in predictions.items():
        for season_name in ["Overall", "spring", "summer", "autumn", "winter"]:
            mask = (np.ones_like(y_test, dtype=bool) if season_name == "Overall" else (test_seasons == season_name))
            sample_size = np.sum(mask)

            if sample_size < 30:
                continue

            y_true_s, y_pred_s = y_test[mask], y_pred[mask]

            metrics = {
                'RMSE': calculate_rmse(y_true_s, y_pred_s),
                'R2': calculate_r2(y_true_s, y_pred_s),
                'NSE': calculate_nse(y_true_s, y_pred_s),
                'Bias': calculate_bias(y_true_s, y_pred_s),
                'MAE': mean_absolute_error(y_true_s, y_pred_s)
            }

            ci_results = {}
            ci_funcs = {'RMSE': calculate_rmse, 'R2': calculate_r2, 'NSE': calculate_nse, 'Bias': calculate_bias,
                        'MAE': mean_absolute_error}
            for name, func in ci_funcs.items():
                low, high, width = blocked_bootstrap_ci(y_true_s, y_pred_s, func, seed=seed)
                ci_results[f'{name}_CI_low'] = low
                ci_results[f'{name}_CI_high'] = high
                ci_results[f'{name}_CI_width'] = width

            site_results.append({
                'Site': file_name, 'InputSet': input_set_name, 'Model': model_name, 'Season': season_name,
                'Sample_Size': sample_size, **metrics, **ci_results
            })

    results_df = pd.DataFrame(site_results)
    predictions_df = pd.DataFrame(
        {'Date': test_data['date'], 'True_ET': y_test, **{f'{k}_Pred': v for k, v in predictions.items()}})

    return results_df, predictions_df


if __name__ == "__main__":
    start_time = time.time()

    FOLDER_PATH = r'X:\X'
    OUTPUT_FOLDER = r'X:\X'
    SEED = 42

    if not os.path.exists(OUTPUT_FOLDER):
        os.makedirs(OUTPUT_FOLDER)

    log_file_path = os.path.join(OUTPUT_FOLDER, 'evaluation_log.log')

    logging.basicConfig(
        level=logging.INFO,
        format='%(asctime)s - %(levelname)s - %(message)s',
        handlers=[
            logging.FileHandler(log_file_path, mode='w', encoding='utf-8'),
            logging.StreamHandler()
        ]
    )

    logging.info("=" * 80)
    logging.info(" TIME-SERIES SAFE STACKING EVALUATION SYSTEM (Publication-Ready Version)")
    logging.info("=" * 80 + "\n")

    INPUT_SETS = {"X": ["X"], "X": ["X"],
                  "X": ["X"]}
    HYPERPARAM_FILES = {"X": r'X:\XXXXXXX\XXXXX',
                        "X": r'X:\XXXXXXX\XXXXX',
                        "X": r'X:\XXXXXXX\XXXXX'}

    logging.info("Validating hyperparameter files...")
    for input_set, path in HYPERPARAM_FILES.items():
        if os.path.exists(path):
            logging.info(f"  [OK] {input_set}: {path}")
        else:
            logging.warning(f"  [NOT FOUND] {input_set}: {path}")
    logging.info("")

    all_results_list = []
    files = [f for f in os.listdir(FOLDER_PATH) if f.endswith('.xlsx')]
    logging.info(f"Found {len(files)} station files to process.\n")

    for input_set_name, features in INPUT_SETS.items():
        logging.info("=" * 80)
        logging.info(f" Starting INPUT SET: {input_set_name} | Features: {', '.join(features)} | Seed: {SEED}")
        logging.info("=" * 80)

        current_hyperparam_file = HYPERPARAM_FILES[input_set_name]
        if not os.path.exists(current_hyperparam_file):
            logging.error(f"Hyperparameter file not found: {current_hyperparam_file}. Skipping {input_set_name}.\n")
            continue
        try:
            hyperparams_df = pd.read_excel(current_hyperparam_file)
            logging.info(f"[OK] Successfully loaded {len(hyperparams_df)} hyperparameter configurations.\n")
        except Exception as e:
            logging.error(f"Failed to load hyperparameter file: {e}. Skipping {input_set_name}.\n")
            continue

        for file in files:
            file_path = os.path.join(FOLDER_PATH, file)
            file_name = os.path.splitext(file)[0]

            if file_name not in hyperparams_df['X'].values:
                logging.warning(f"    [WARNING] No hyperparameters found for site {file_name}. Skipping.")
                continue

            hyperparams_for_site = hyperparams_df[hyperparams_df['文件名'] == file_name]


            logging.info(f"      > Verifying Hyperparameters for Site: {file_name}\n{hyperparams_for_site.to_string()}")

            try:
                results_df, predictions_df = evaluate_models_for_site(
                    file_path, hyperparams_for_site, features, input_set_name, SEED
                )

                if not results_df.empty:
                    all_results_list.append(results_df)
                    per_site_path = os.path.join(OUTPUT_FOLDER, f"site_{file_name}_{input_set_name}_results.csv")
                    results_df.to_csv(per_site_path, index=False)
                    logging.info(f"      [OK] Saved per-site results to: {os.path.basename(per_site_path)}")

                if not predictions_df.empty:
                    per_site_pred_path = os.path.join(OUTPUT_FOLDER,
                                                      f"site_{file_name}_{input_set_name}_predictions.csv")
                    predictions_df.to_csv(per_site_pred_path, index=False)
                    logging.info(f"      [OK] Saved per-site predictions to: {os.path.basename(per_site_pred_path)}")

            except Exception as e:
                logging.error(f"      [FATAL ERROR] processing {file_name}: {e}")
                logging.error(traceback.format_exc())
                continue

    if not all_results_list:
        logging.warning("\n[ERROR] No results were generated. Please check logs for errors.")
    else:
        logging.info("\n" + "=" * 80)
        logging.info(" AGGREGATING AND SAVING FINAL RESULTS")
        logging.info("=" * 80 + "\n")

        try:
            final_results_df = pd.concat(all_results_list, ignore_index=True)
            agg_path = os.path.join(OUTPUT_FOLDER, 'aggregated_metrics.csv')
            final_results_df.to_csv(agg_path, index=False)
            logging.info(f"[OK] Successfully aggregated {len(final_results_df)} records.")
            logging.info(f"[OK] Final aggregated metrics saved to: {agg_path}")
        except Exception as e:
            logging.error(f"\n[ERROR] Aggregation failed: {e}")
            logging.error(traceback.format_exc())

    end_time = time.time()
    logging.info("\n" + "=" * 80)
    logging.info(" EXECUTION SUMMARY")
    logging.info("=" * 80)
    logging.info(f"  Total execution time: {(end_time - start_time) / 60:.2f} minutes.")
    logging.info(f"  Total sites processed: {len(files)}")
    logging.info(f"  Input sets evaluated: {', '.join(INPUT_SETS.keys())}")
    logging.info(f"  Output directory: {OUTPUT_FOLDER}")
    logging.info(f"  Detailed log saved to: {log_file_path}")
    logging.info("=" * 80)