AcidificationModel/api/model_optimize/model_analyzer.py

# -*- coding: utf-8 -*-
"""
模型分析工具

此模块提供了综合的模型分析功能，包括：
- 模型性能详细分析
- 特征重要性分析
- 预测行为分析
- 参数敏感性分析
- 模型诊断和可视化

"""

import os
import sys
import pickle
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.model_selection import cross_val_score
from sklearn.inspection import permutation_importance
import sqlite3
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class ModelAnalyzer:
    """
    模型分析器类
    
    提供全面的模型分析功能，支持模型性能评估、特征分析、预测分析等
    """
    
    def __init__(self, db_path='SoilAcidification.db'):
        """
        初始化模型分析器
        
        @param {str} db_path - 数据库文件路径
        """
        self.db_path = db_path
        self.model_info = None
        self.model = None
        self.test_data = None
        self.predictions = None
        self.analysis_results = {}
        
        # 创建输出目录
        self.output_dir = 'model_optimize/analysis_results'
        if not os.path.exists(self.output_dir):
            os.makedirs(self.output_dir)
    
    def load_model_by_id(self, model_id):
        """
        根据模型ID加载模型信息和模型对象
        
        @param {int} model_id - 模型ID
        @return {bool} - 是否成功加载
        """
        try:
            # 连接数据库
            conn = sqlite3.connect(self.db_path)
            cursor = conn.cursor()
            
            # 查询模型信息
            cursor.execute("""
                SELECT ModelID, Model_name, Model_type, ModelFilePath, 
                       Data_type, Performance_score, MAE, RMSE, CV_score,
                       Description, Created_at
                FROM Models 
                WHERE ModelID = ?
            """, (model_id,))
            
            result = cursor.fetchone()
            if not result:
                print(f"未找到ID为 {model_id} 的模型")
                return False
            
            # 存储模型信息
            self.model_info = {
                'ModelID': result[0],
                'Model_name': result[1],
                'Model_type': result[2],
                'ModelFilePath': result[3],
                'Data_type': result[4],
                'Performance_score': result[5],
                'MAE': result[6],
                'RMSE': result[7],
                'CV_score': result[8],
                'Description': result[9],
                'Created_at': result[10]
            }
            
            conn.close()
            
            # 加载模型文件
            if os.path.exists(self.model_info['ModelFilePath']):
                with open(self.model_info['ModelFilePath'], 'rb') as f:
                    self.model = pickle.load(f)
                print(f"成功加载模型: {self.model_info['Model_name']}")
                return True
            else:
                print(f"模型文件不存在: {self.model_info['ModelFilePath']}")
                return False
                
        except Exception as e:
            print(f"加载模型时出错: {str(e)}")
            return False
    
    def load_test_data(self):
        """
        根据模型数据类型加载相应的测试数据
        
        @return {bool} - 是否成功加载测试数据
        """
        try:
            data_type = self.model_info['Data_type']
            
            if data_type == 'reflux':
                X_test_path = 'uploads/data/X_test_reflux.csv'
                Y_test_path = 'uploads/data/Y_test_reflux.csv'
            elif data_type == 'reduce':
                X_test_path = 'uploads/data/X_test_reduce.csv'
                Y_test_path = 'uploads/data/Y_test_reduce.csv'
            else:
                print(f"不支持的数据类型: {data_type}")
                return False
            
            if os.path.exists(X_test_path) and os.path.exists(Y_test_path):
                self.X_test = pd.read_csv(X_test_path)
                self.Y_test = pd.read_csv(Y_test_path)
                print(f"成功加载 {data_type} 类型的测试数据")
                print(f"测试集大小: X_test {self.X_test.shape}, Y_test {self.Y_test.shape}")
                return True
            else:
                print(f"测试数据文件不存在")
                return False
                
        except Exception as e:
            print(f"加载测试数据时出错: {str(e)}")
            return False
    
    def load_training_data(self):
        """
        加载训练数据以获取参数的真实分布范围
        
        @return {bool} - 是否成功加载训练数据
        """
        try:
            data_type = self.model_info['Data_type']
            
            if data_type == 'reflux':
                # 加载current_reflux表的数据
                conn = sqlite3.connect(self.db_path)
                training_data = pd.read_sql_query("SELECT * FROM current_reflux", conn)
                conn.close()
                
                # 分离特征和目标变量
                if 'Delta_pH' in training_data.columns:
                    self.X_train = training_data.drop(['id', 'Delta_pH'], axis=1, errors='ignore')
                    self.y_train = training_data['Delta_pH']
                else:
                    print("未找到目标变量 Delta_pH")
                    return False
                    
            elif data_type == 'reduce':
                # 加载current_reduce表的数据
                conn = sqlite3.connect(self.db_path)
                training_data = pd.read_sql_query("SELECT * FROM current_reduce", conn)
                conn.close()
                
                # 分离特征和目标变量
                if 'Q_over_b' in training_data.columns:
                    self.X_train = training_data.drop(['id', 'Q_over_b'], axis=1, errors='ignore')
                    self.y_train = training_data['Q_over_b']
                else:
                    print("未找到目标变量 Q_over_b")
                    return False
            else:
                print(f"不支持的数据类型: {data_type}")
                return False
            
            print(f"成功加载训练数据: {self.X_train.shape}")
            print(f"训练数据特征统计:")
            print(self.X_train.describe())
            return True
            
        except Exception as e:
            print(f"加载训练数据时出错: {str(e)}")
            return False
    
    def calculate_detailed_metrics(self):
        """
        计算详细的模型性能指标
        
        @return {dict} - 包含各种性能指标的字典
        """
        try:
            if self.model is None or self.X_test is None:
                print("模型或测试数据未加载")
                return {}
            
            # 获取预测值
            y_pred = self.model.predict(self.X_test)
            y_true = self.Y_test.iloc[:, 0] if len(self.Y_test.columns) == 1 else self.Y_test
            
            # 计算各种指标
            metrics = {
                'R2_Score': r2_score(y_true, y_pred),
                'RMSE': mean_squared_error(y_true, y_pred, squared=False),
                'MAE': mean_absolute_error(y_true, y_pred),
                'MSE': mean_squared_error(y_true, y_pred),
                'MAPE': np.mean(np.abs((y_true - y_pred) / y_true)) * 100,
                'Max_Error': np.max(np.abs(y_true - y_pred)),
                'Mean_Residual': np.mean(y_true - y_pred),
                'Std_Residual': np.std(y_true - y_pred)
            }
            
            self.predictions = {'y_true': y_true, 'y_pred': y_pred}
            self.analysis_results['metrics'] = metrics
            
            # 中文指标名称映射
            metric_names_zh = {
                'R2_Score': 'R²得分',
                'RMSE': '均方根误差',
                'MAE': '平均绝对误差',
                'MSE': '均方误差',
                'MAPE': '平均绝对百分比误差',
                'Max_Error': '最大误差',
                'Mean_Residual': '平均残差',
                'Std_Residual': '残差标准差'
            }
            
            print("=== 模型性能指标 ===")
            for key, value in metrics.items():
                zh_name = metric_names_zh.get(key, key)
                print(f"{key} ({zh_name}): {value:.4f}")
            
            return metrics
            
        except Exception as e:
            print(f"计算性能指标时出错: {str(e)}")
            return {}
    
    def analyze_feature_importance(self):
        """
        分析特征重要性
        
        @return {dict} - 特征重要性分析结果
        """
        try:
            if self.model is None:
                print("模型未加载")
                return {}
            
            importance_data = {}
            
            # 1. 模型内置特征重要性（如果支持）
            if hasattr(self.model, 'feature_importances_'):
                importance_data['model_importance'] = {
                    'features': self.X_test.columns.tolist(),
                    'importance': self.model.feature_importances_.tolist()
                }
            
            # 2. 排列特征重要性
            if self.X_test is not None and self.Y_test is not None:
                y_true = self.Y_test.iloc[:, 0] if len(self.Y_test.columns) == 1 else self.Y_test
                perm_importance = permutation_importance(
                    self.model, self.X_test, y_true, 
                    n_repeats=10, random_state=42
                )
                importance_data['permutation_importance'] = {
                    'features': self.X_test.columns.tolist(),
                    'importance_mean': perm_importance.importances_mean.tolist(),
                    'importance_std': perm_importance.importances_std.tolist()
                }
            
            self.analysis_results['feature_importance'] = importance_data
            
            print("=== 特征重要性分析完成 ===")
            return importance_data
            
        except Exception as e:
            print(f"分析特征重要性时出错: {str(e)}")
            return {}
    
    def analyze_prediction_behavior(self, sample_params=None):
        """
        分析模型的预测行为
        
        @param {dict} sample_params - 示例参数，如用户提供的参数
        @return {dict} - 预测行为分析结果
        """
        try:
            if self.model is None:
                print("模型未加载")
                return {}
            
            behavior_analysis = {}
            
            # 1. 分析用户提供的示例参数
            if sample_params:
                # 转换参数为DataFrame
                if self.model_info['Data_type'] == 'reflux':
                    # 为reflux类型，移除target_pH
                    sample_params_clean = sample_params.copy()
                    sample_params_clean.pop('target_pH', None)
                    sample_df = pd.DataFrame([sample_params_clean])
                else:
                    sample_df = pd.DataFrame([sample_params])
                
                # 确保列顺序与训练数据一致
                if hasattr(self, 'X_test'):
                    sample_df = sample_df.reindex(columns=self.X_test.columns, fill_value=0)
                
                prediction = self.model.predict(sample_df)[0]
                
                behavior_analysis['sample_prediction'] = {
                    'input_params': sample_params,
                    'prediction': float(prediction),
                    'model_type': self.model_info['Data_type']
                }
                
                print(f"=== 示例参数预测结果 ===")
                print(f"输入参数: {sample_params}")
                print(f"预测结果: {prediction:.4f}")
            
            # 2. 预测值分布分析
            if self.predictions:
                y_pred = self.predictions['y_pred']
                behavior_analysis['prediction_distribution'] = {
                    'mean': float(np.mean(y_pred)),
                    'std': float(np.std(y_pred)),
                    'min': float(np.min(y_pred)),
                    'max': float(np.max(y_pred)),
                    'percentiles': {
                        '25%': float(np.percentile(y_pred, 25)),
                        '50%': float(np.percentile(y_pred, 50)),
                        '75%': float(np.percentile(y_pred, 75))
                    }
                }
            
            self.analysis_results['prediction_behavior'] = behavior_analysis
            return behavior_analysis
            
        except Exception as e:
            print(f"分析预测行为时出错: {str(e)}")
            return {}
    
    def sensitivity_analysis(self, sample_params=None):
        """
        基于训练数据分布的参数敏感性分析
        
        @param {dict} sample_params - 基准参数
        @return {dict} - 敏感性分析结果
        """
        try:
            if self.model is None or sample_params is None:
                print("模型未加载或未提供基准参数")
                return {}
            
            # 确保训练数据已加载
            if not hasattr(self, 'X_train') or self.X_train is None:
                print("训练数据未加载，尝试加载...")
                if not self.load_training_data():
                    print("无法加载训练数据，使用默认参数范围")
                    return self.sensitivity_analysis_fallback(sample_params)
            
            sensitivity_results = {}
            
            # 准备基准数据
            base_params = sample_params.copy()
            if self.model_info['Data_type'] == 'reflux':
                base_params.pop('target_pH', None)
            
            base_df = pd.DataFrame([base_params])
            if hasattr(self, 'X_test'):
                base_df = base_df.reindex(columns=self.X_test.columns, fill_value=0)
            
            base_prediction = self.model.predict(base_df)[0]
            print(f"基准预测值: {base_prediction:.6f}")
            
            # 对每个参数进行基于训练数据的敏感性分析
            for param_name, base_value in base_params.items():
                if param_name not in base_df.columns:
                    continue
                
                # 检查参数是否在训练数据中存在
                if param_name not in self.X_train.columns:
                    print(f"警告: 参数 {param_name} 不在训练数据中，跳过")
                    continue
                
                param_sensitivity = {
                    'base_value': base_value,
                    'variations': [],
                    'predictions': [],
                    'sensitivity_score': 0,
                    'training_stats': {}
                }
                
                # 获取训练数据中该参数的统计信息
                param_series = self.X_train[param_name]
                param_stats = {
                    'min': float(param_series.min()),
                    'max': float(param_series.max()),
                    'mean': float(param_series.mean()),
                    'std': float(param_series.std()),
                    'q25': float(param_series.quantile(0.25)),
                    'q50': float(param_series.quantile(0.50)),
                    'q75': float(param_series.quantile(0.75)),
                    'q05': float(param_series.quantile(0.05)),
                    'q95': float(param_series.quantile(0.95))
                }
                param_sensitivity['training_stats'] = param_stats
                
                print(f"\n分析参数: {param_name} (基准值: {base_value})")
                print(f"  训练数据范围: [{param_stats['min']:.3f}, {param_stats['max']:.3f}]")
                print(f"  训练数据均值±标准差: {param_stats['mean']:.3f} ± {param_stats['std']:.3f}")
                print(f"  训练数据分位数 [5%, 25%, 50%, 75%, 95%]: [{param_stats['q05']:.3f}, {param_stats['q25']:.3f}, {param_stats['q50']:.3f}, {param_stats['q75']:.3f}, {param_stats['q95']:.3f}]")
                
                # 基于训练数据分布设计参数变化策略
                # 策略1: 使用训练数据的实际范围
                min_val = param_stats['min']
                max_val = param_stats['max']
                mean_val = param_stats['mean']
                std_val = param_stats['std']
                
                # 创建合理的参数变化点
                variations = []
                
                # 1. 训练数据的关键分位数点
                variations.extend([
                    param_stats['min'],     # 最小值
                    param_stats['q05'],     # 5%分位数
                    param_stats['q25'],     # 25%分位数
                    param_stats['q50'],     # 中位数
                    param_stats['q75'],     # 75%分位数
                    param_stats['q95'],     # 95%分位数
                    param_stats['max'],     # 最大值
                ])
                
                # 2. 围绕均值的标准差变化
                variations.extend([
                    max(min_val, mean_val - 2*std_val),  # 均值-2σ
                    max(min_val, mean_val - std_val),    # 均值-σ
                    mean_val,                            # 均值
                    min(max_val, mean_val + std_val),    # 均值+σ
                    min(max_val, mean_val + 2*std_val),  # 均值+2σ
                ])
                
                # 3. 基准值本身
                variations.append(base_value)
                
                # 4. 如果基准值在合理范围内，添加围绕基准值的小幅变化
                if min_val <= base_value <= max_val:
                    # 计算基准值的相对变化（不超出训练数据范围）
                    range_size = max_val - min_val
                    small_change = range_size * 0.05  # 5%的范围变化
                    
                    variations.extend([
                        max(min_val, base_value - small_change),
                        min(max_val, base_value + small_change),
                    ])
                
                # 去重并排序
                variations = sorted(list(set(variations)))
                
                print(f"  使用参数变化点: {len(variations)} 个")
                
                for i, variation in enumerate(variations):
                    temp_params = base_params.copy()
                    temp_params[param_name] = variation
                    temp_df = pd.DataFrame([temp_params])
                    temp_df = temp_df.reindex(columns=self.X_test.columns, fill_value=0)
                    
                    try:
                        pred = self.model.predict(temp_df)[0]
                        param_sensitivity['variations'].append(float(variation))
                        param_sensitivity['predictions'].append(float(pred))
                        
                        # 计算相对于训练数据范围的位置
                        if max_val > min_val:
                            relative_pos = (variation - min_val) / (max_val - min_val) * 100
                        else:
                            relative_pos = 50.0
                        
                        pred_change = pred - base_prediction
                        print(f"  {param_name}={variation:.3f} (训练数据{relative_pos:.0f}位) → 预测={pred:.6f} (变化={pred_change:+.6f})")
                        
                    except Exception as e:
                        print(f"  警告: {param_name}={variation} 预测失败: {str(e)}")
                        continue
                
                # 计算敏感性指标
                if len(param_sensitivity['predictions']) > 1:
                    predictions = np.array(param_sensitivity['predictions'])
                    
                    # 指标1: 预测值变化范围
                    pred_range = np.max(predictions) - np.min(predictions)
                    
                    # 指标2: 预测值标准差
                    pred_std = np.std(predictions)
                    
                    # 指标3: 与基准值的最大偏差
                    max_deviation = np.max(np.abs(predictions - base_prediction))
                    
                    # 指标4: 标准化敏感性（考虑参数的变化范围）
                    param_range = max_val - min_val if max_val > min_val else 1.0
                    normalized_sensitivity = pred_range / param_range if param_range > 0 else 0
                    
                    # 使用最大偏差作为主要敏感性指标
                    param_sensitivity['sensitivity_score'] = float(max_deviation)
                    param_sensitivity['pred_range'] = float(pred_range)
                    param_sensitivity['pred_std'] = float(pred_std)
                    param_sensitivity['normalized_sensitivity'] = float(normalized_sensitivity)
                    
                    print(f"  {param_name} 敏感性得分: {max_deviation:.6f}")
                    print(f"    预测范围: {pred_range:.6f}, 标准差: {pred_std:.6f}")
                    print(f"    标准化敏感性: {normalized_sensitivity:.6f}")
                else:
                    param_sensitivity['sensitivity_score'] = 0.0
                
                sensitivity_results[param_name] = param_sensitivity
            
            # 按敏感性得分排序
            sorted_sensitivity = dict(sorted(
                sensitivity_results.items(), 
                key=lambda x: x[1]['sensitivity_score'], 
                reverse=True
            ))
            
            self.analysis_results['sensitivity_analysis'] = sorted_sensitivity
            
            print("\n=== 基于训练数据的参数敏感性分析结果 ===")
            for param, result in sorted_sensitivity.items():
                score = result['sensitivity_score']
                normalized_score = result.get('normalized_sensitivity', 0)
                if score > 0.001:
                    sensitivity_level = "高" if score > 0.1 else "中" if score > 0.01 else "低"
                else:
                    sensitivity_level = "极低"
                print(f"{param:12}: 敏感性得分 = {score:.6f} ({sensitivity_level})")
                print(f"             标准化敏感性 = {normalized_score:.6f}")
            
            return sorted_sensitivity
            
        except Exception as e:
            print(f"敏感性分析时出错: {str(e)}")
            import traceback
            print(traceback.format_exc())
            # 如果分析失败，回退到原来的方法
            return self.sensitivity_analysis_fallback(sample_params)
    
    def generate_visualizations(self):
        """
        生成分析可视化图表
        """
        try:
            timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
            model_id = self.model_info['ModelID']
            
            # 1. 预测值 vs 真实值散点图
            if self.predictions:
                plt.figure(figsize=(10, 8))
                
                y_true = self.predictions['y_true']
                y_pred = self.predictions['y_pred']
                
                plt.scatter(y_true, y_pred, alpha=0.6, s=50)
                
                # 绘制理想线 (y=x)
                min_val = min(min(y_true), min(y_pred))
                max_val = max(max(y_true), max(y_pred))
                plt.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Prediction')
                
                plt.xlabel('真实值 (True Values)')
                plt.ylabel('预测值 (Predicted Values)')
                plt.title(f'模型 {model_id} 预测性能\nR^2 = {self.analysis_results["metrics"]["R2_Score"]:.4f}')
                plt.legend()
                plt.grid(True, alpha=0.3)
                
                plt.tight_layout()
                plt.savefig(f'{self.output_dir}/model_{model_id}_scatter_{timestamp}.png', dpi=300)
                plt.show()
                
                # 2. 残差图
                plt.figure(figsize=(12, 5))
                
                # 确保数据类型一致，转换为numpy数组
                y_true_array = np.array(y_true).flatten()
                y_pred_array = np.array(y_pred).flatten()
                residuals = y_true_array - y_pred_array
                
                # 移除NaN值
                valid_indices = ~(np.isnan(residuals) | np.isnan(y_pred_array))
                residuals_clean = residuals[valid_indices]
                y_pred_clean = y_pred_array[valid_indices]
                
                # 残差散点图
                plt.subplot(1, 2, 1)
                plt.scatter(y_pred_clean, residuals_clean, alpha=0.6)
                plt.axhline(y=0, color='r', linestyle='--', linewidth=2)
                plt.xlabel('预测值')
                plt.ylabel('残差 (真实值 - 预测值)')
                plt.title('残差散点图')
                plt.grid(True, alpha=0.3)
                
                # 残差直方图
                plt.subplot(1, 2, 2)
                # 根据数据量自动调整bins数量
                n_samples = len(residuals_clean)
                n_bins = min(30, max(5, int(np.sqrt(n_samples))))  # Sturges' rule的变体
                
                n, bins, patches = plt.hist(residuals_clean, bins=n_bins, alpha=0.7, edgecolor='black', density=False)
                plt.xlabel('残差')
                plt.ylabel('频数')
                plt.title(f'残差分布 (n={n_samples}, bins={n_bins})')
                plt.grid(True, alpha=0.3)
                
                # 添加统计信息
                mean_residual = np.mean(residuals_clean)
                std_residual = np.std(residuals_clean)
                plt.axvline(mean_residual, color='red', linestyle='--', alpha=0.8, 
                           label=f'均值: {mean_residual:.4f}')
                plt.axvline(mean_residual + std_residual, color='orange', linestyle=':', alpha=0.8,
                           label=f'±1σ: ±{std_residual:.4f}')
                plt.axvline(mean_residual - std_residual, color='orange', linestyle=':', alpha=0.8)
                plt.legend(fontsize=8)
                
                plt.tight_layout()
                plt.savefig(f'{self.output_dir}/model_{model_id}_residuals_{timestamp}.png', dpi=300)
                plt.show()
            
            # 3. 特征重要性图
            if 'feature_importance' in self.analysis_results:
                importance_data = self.analysis_results['feature_importance']
                
                if 'model_importance' in importance_data:
                    plt.figure(figsize=(10, 6))
                    features = importance_data['model_importance']['features']
                    importance = importance_data['model_importance']['importance']
                    
                    indices = np.argsort(importance)[::-1]
                    
                    plt.bar(range(len(features)), [importance[i] for i in indices])
                    plt.xticks(range(len(features)), [features[i] for i in indices], rotation=45)
                    plt.xlabel('特征')
                    plt.ylabel('重要性')
                    plt.title(f'模型 {model_id} 特征重要性')
                    plt.tight_layout()
                    plt.savefig(f'{self.output_dir}/model_{model_id}_feature_importance_{timestamp}.png', dpi=300)
                    plt.show()
            
            # 4. 敏感性分析图
            if 'sensitivity_analysis' in self.analysis_results:
                sensitivity_data = self.analysis_results['sensitivity_analysis']
                
                # 敏感性得分条形图
                plt.figure(figsize=(10, 6))
                params = list(sensitivity_data.keys())
                scores = [sensitivity_data[param]['sensitivity_score'] for param in params]
                
                plt.bar(params, scores)
                plt.xlabel('参数')
                plt.ylabel('敏感性得分')
                plt.title(f'模型 {model_id} 参数敏感性分析')
                plt.xticks(rotation=45)
                plt.tight_layout()
                plt.savefig(f'{self.output_dir}/model_{model_id}_sensitivity_{timestamp}.png', dpi=300)
                plt.show()
                
                # 详细敏感性曲线
                n_params = len(params)
                fig, axes = plt.subplots(2, 3, figsize=(15, 10))
                axes = axes.flatten()
                
                for i, param in enumerate(params[:6]):  # 显示前6个参数
                    if i < len(axes):
                        data = sensitivity_data[param]
                        axes[i].plot(data['variations'], data['predictions'], 'b-o', markersize=4)
                        axes[i].set_xlabel(param)
                        axes[i].set_ylabel('预测值')
                        axes[i].set_title(f'{param} 敏感性')
                        axes[i].grid(True, alpha=0.3)
                
                # 隐藏多余的子图
                for i in range(len(params), len(axes)):
                    axes[i].set_visible(False)
                
                plt.tight_layout()
                plt.savefig(f'{self.output_dir}/model_{model_id}_sensitivity_curves_{timestamp}.png', dpi=300)
                plt.show()
            
            print(f"=== 可视化图表已保存到 {self.output_dir} ===")
            
        except Exception as e:
            print(f"生成可视化时出错: {str(e)}")
    
    def generate_report(self):
        """
        生成综合分析报告
        
        @return {str} - 报告内容
        """
        try:
            timestamp = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
            model_id = self.model_info['ModelID']
            
            report = f"""
模型分析报告
{'='*60}

生成时间: {timestamp}
分析模型: {self.model_info['Model_name']} (ID: {model_id})
模型类型: {self.model_info['Model_type']}
数据类型: {self.model_info['Data_type']}
创建时间: {self.model_info['Created_at']}

模型描述:
{self.model_info['Description'] or '无描述'}

{'='*60}
性能指标分析
{'='*60}
"""
            
            if 'metrics' in self.analysis_results:
                metrics = self.analysis_results['metrics']
                for key, value in metrics.items():
                    report += f"{key:15}: {value:8.4f}\n"
            
            report += f"\n{'='*60}\n特征重要性分析\n{'='*60}\n"
            
            if 'feature_importance' in self.analysis_results:
                importance_data = self.analysis_results['feature_importance']
                if 'model_importance' in importance_data:
                    features = importance_data['model_importance']['features']
                    importance = importance_data['model_importance']['importance']
                    
                    # 按重要性排序
                    sorted_indices = np.argsort(importance)[::-1]
                    
                    for i in sorted_indices:
                        report += f"{features[i]:15}: {importance[i]:8.4f}\n"
            
            report += f"\n{'='*60}\n敏感性分析\n{'='*60}\n"
            
            if 'sensitivity_analysis' in self.analysis_results:
                sensitivity_data = self.analysis_results['sensitivity_analysis']
                for param, data in sensitivity_data.items():
                    report += f"{param:15}: {data['sensitivity_score']:.8f}\n"
            
            if 'prediction_behavior' in self.analysis_results:
                behavior = self.analysis_results['prediction_behavior']
                if 'sample_prediction' in behavior:
                    sample = behavior['sample_prediction']
                    report += f"\n{'='*60}\n示例预测分析\n{'='*60}\n"
                    report += f"输入参数: {sample['input_params']}\n"
                    report += f"预测结果: {sample['prediction']:.4f}\n"
            
            report += f"\n{'='*60}\n分析建议\n{'='*60}\n"
            
            # 生成分析建议
            if 'metrics' in self.analysis_results:
                r2 = self.analysis_results['metrics']['R2_Score']
                if r2 > 0.9:
                    report += "• 模型表现优秀，R²得分大于0.9\n"
                elif r2 > 0.8:
                    report += "• 模型表现良好，R²得分在0.8-0.9之间\n"
                elif r2 > 0.7:
                    report += "• 模型表现中等，R²得分在0.7-0.8之间，建议优化\n"
                else:
                    report += "• 模型表现较差，R²得分低于0.7，需要重新训练\n"
            
            if 'sensitivity_analysis' in self.analysis_results:
                sensitivity_data = self.analysis_results['sensitivity_analysis']
                most_sensitive = max(sensitivity_data.items(), key=lambda x: x[1]['sensitivity_score'])
                report += f"• 最敏感参数: {most_sensitive[0]}，在调参时需要特别注意\n"
                
                least_sensitive = min(sensitivity_data.items(), key=lambda x: x[1]['sensitivity_score'])
                report += f"• 最不敏感参数: {least_sensitive[0]}，对预测结果影响较小\n"
            
            report += f"\n{'='*60}\n报告结束\n{'='*60}\n"
            
            # 保存报告
            report_path = f'{self.output_dir}/model_{model_id}_analysis_report_{datetime.now().strftime("%Y%m%d_%H%M%S")}.txt'
            with open(report_path, 'w', encoding='utf-8') as f:
                f.write(report)
            
            print(f"分析报告已保存到: {report_path}")
            return report
            
        except Exception as e:
            print(f"生成报告时出错: {str(e)}")
            return ""
    
    def run_full_analysis(self, model_id, sample_params=None):
        """
        运行完整的模型分析流程
        
        @param {int} model_id - 模型ID
        @param {dict} sample_params - 示例参数
        @return {dict} - 完整的分析结果
        """
        print(f"开始分析模型 {model_id}...")
        
        # 1. 加载模型
        if not self.load_model_by_id(model_id):
            return {}
        
        # 2. 加载测试数据
        if not self.load_test_data():
            print("警告: 无法加载测试数据，部分分析功能将不可用")
        
        # 3. 加载训练数据（用于敏感性分析）
        if not self.load_training_data():
            print("警告: 无法加载训练数据，敏感性分析将使用默认方法")
        
        # 4. 计算性能指标
        self.calculate_detailed_metrics()
        
        # 5. 特征重要性分析
        self.analyze_feature_importance()
        
        # 6. 预测行为分析
        self.analyze_prediction_behavior(sample_params)
        
        # 7. 敏感性分析
        if sample_params:
            self.sensitivity_analysis(sample_params)
        
        # 8. 生成可视化
        self.generate_visualizations()
        
        # 9. 生成报告
        self.generate_report()
        
        print("分析完成!")
        return self.analysis_results


def main():
    """
    主函数 - 分析模型24的示例
    """
    # 创建分析器实例
    analyzer = ModelAnalyzer()
    
    # 用户提供的示例参数
    sample_params = {
        "OM": 19.913,
        "CL": 373.600,
        "CEC": 7.958,
        "H_plus": 0.774,
        "N": 0.068,
        "Al3_plus": 3.611,
        "target_pH": 7.0  # 这个参数在reflux模型中会被移除
    }
    
    # 运行完整分析
    results = analyzer.run_full_analysis(model_id=24, sample_params=sample_params)
    
    if results:
        print("\n=== 分析完成 ===")
        print("所有分析结果和可视化图表已保存到 model_optimize/analysis_results 目录")
    else:
        print("分析失败，请检查模型ID和数据文件")


if __name__ == "__main__":
    main()