Function_BO_2.py

import numpy as np
import matplotlib.pyplot as plt
from skopt.learning import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
from skopt.acquisition import gaussian_ei
from scipy.optimize import minimize 
import pandas as pd

def runBO(df, input_cols, output_col, bounds=None):
    
    X_in = df[input_cols].values
    Y_out = df[output_col].values
    
    # 建议加上 log，因为熵产通常不能为负，且变化范围大
    # 如果数据里有负数，请注释掉这一行
    # Y_out = np.log(Y_out) 

    # =======================================================
    # 🛠️ 2. 边界处理逻辑
    # =======================================================
    if bounds is None:
        print("⚠️ 未提供边界，正在根据当前数据自动生成...")
        # 自动获取每一列的 min-10% 和 max+10%
        bounds = [ (df[col].min()*0.9, df[col].max()*1.1) for col in input_cols ]
    
    if len(bounds) != len(input_cols):
        raise ValueError(f"错误：变量数 {len(input_cols)} 与边界数 {len(bounds)} 不一致")

    print(f"Running BO with inputs: {input_cols}")
    print(f"Search Bounds: {bounds}")

    # 拟合 GP 模型
    kernel = Matern(nu=2.5)
    gpr_model = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, random_state=0)
    gpr_model.fit(X_in, Y_out)
    y_opt = np.min(Y_out) # 当前观测到的最小值
    
    next_point = []

    # ================= CASE A: 2D 可视化模式 =================
    if len(input_cols) == 2:
        print(">>> 2D 模式：正在绘图...")
        
        # 1. 创建网格 (由 bounds 决定范围)
        res = 50 # 网格密度 (太大画图会卡)
        x1 = np.linspace(bounds[0][0], bounds[0][1], res) 
        x2 = np.linspace(bounds[1][0], bounds[1][1], res)
        X1, X2 = np.meshgrid(x1, x2)
        X_grid = np.vstack((X1.ravel(), X2.ravel())).T
        
        # 2. 【关键】计算模型预测值 (您之前漏了这一步)
        pred_mean, pred_std = gpr_model.predict(X_grid, return_std=True)
        
        # 3. 计算 EI (采集函数)
        ei_values = gaussian_ei(X=X_grid, model=gpr_model, y_opt=y_opt)
        
        # 4. 找到 EI 最大的点
        best_idx = np.argmax(ei_values)
        best_x = X_grid[best_idx]
        next_point = best_x
        
        # 5. predict next Obj value
        pred_val = gpr_model.predict(np.array([best_x]))


        # ==========================================
        # 🎨 【核心修复】这里是画图代码
        # ==========================================
        fig = plt.figure(figsize=(16, 7))

        # --- 子图 1: GP 模型预测均值 ---
        ax1 = fig.add_subplot(121, projection='3d')
        surf1 = ax1.plot_surface(X1, X2, pred_mean.reshape(X1.shape), 
                                                                    cmap='viridis',
                                                                    alpha=0.7, 
                                                                    edgecolor='none', 
                                                                    linewidth=0)
        ax1.scatter(X_in[:,0], X_in[:,1], Y_out, 
                                                c='red', 
                                                s=50, 
                                                edgecolors='k', 
                                                label='Observed Data')
        
        # 💡 强制锁定坐标轴范围
        ax1.set_xlim(bounds[0][0], bounds[0][1])
        ax1.set_ylim(bounds[1][0], bounds[1][1])
        
        ax1.set_title(f'GP Model Prediction (Mean)\nTarget: {output_col}')
        ax1.set_xlabel(input_cols[0])
        ax1.set_ylabel(input_cols[1])
        fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=10)

        # --- 子图 2: EI 采集函数 ---
        ax2 = fig.add_subplot(122, projection='3d')
        # 画 EI 曲面
        surf2 = ax2.plot_surface(X1, X2, ei_values.reshape(X1.shape), 
                                                                    cmap='plasma', 
                                                                    alpha=0.8, 
                                                                    edgecolor='none')
        # 标记推荐的下一个点 (绿星)
        ax2.scatter(best_x[0], best_x[1], np.max(ei_values), 
                                                             c='lime', 
                                                            s=200, 
                                                            marker='*', 
                                                            edgecolors='k', 
                                                            label='Next Point', 
                                                            zorder=10)
        
        ax2.set_title('Acquisition Function (Expected Improvement)')
        ax2.set_xlabel(input_cols[0])
        ax2.set_ylabel(input_cols[1])
        ax2.legend()


        # for ax in [ax1, ax2]:
        #     # 1. 关闭网格线
        #     ax.grid(False)
            
        #     # 2. 隐藏 Z 轴（既然是俯视，Z 轴会重叠，隐藏后更美观）
        #     ax.set_zticks([])
        #     ax.zaxis.line.set_lw(0.) # 隐藏 Z 轴的那条线
            
        #     # 3. 设置初始镜头：正上方俯视
        #     # elev=90 表示仰角 90 度（正上方），azim=-90 是常用的方位角设置
        #     ax.view_init(elev=90, azim=-90)

        plt.tight_layout()
        plt.show() # 这一句必须要有，图才会弹出来

    # ================= CASE B: 高维数值优化模式 =================
    else:
        def min_func(x):
            return -gaussian_ei(X=x.reshape(1, -1), model=gpr_model, y_opt=y_opt)
        
        best_ei_val = np.inf
        for i in range(10): 
            x0 = [np.random.uniform(b[0], b[1]) for b in bounds]
            res = minimize(min_func, x0=x0, bounds=bounds, method='L-BFGS-B')
            if res.fun < best_ei_val:
                best_ei_val = res.fun
                next_point = res.x

    # 返回结果
    next_point_dict = {col: val for col, val in zip(input_cols, next_point)}
    next_point_dict["pred_obj"] = pred_val[0]
    return next_point_dict


if __name__ =="__main__":
    # 您的调试代码保持不变
    # num_samples = 10
    # width_data = np.round(np.random.uniform(0.2, 2.0, num_samples), 4)
    # flow_data  = np.round(np.random.uniform(0.1, 1.0, num_samples), 4)
    # entropy_data = np.round(np.random.uniform(0.5, 100.0, num_samples), 4)

    # df_debug = pd.DataFrame({
    #                         "Channel_Width": width_data,
    #                         "Flow_Rate": flow_data,
    #                         "Entropy": entropy_data
    #                                                 })

    # input_cols = ["Channel_Width", "Flow_Rate"]
    # output_col = "Entropy"
    # debug_bounds = [(0.2, 2.0), (0.1, 1.0)]

    # print("🚀 开始测试 runBO ...")
    # next_point = runBO(df_debug, input_cols, output_col, bounds=debug_bounds)
    # print("推荐点:", next_point)





    import pandas as pd

    df = pd.read_csv(r'C:\Users\lichengd\OneDrive - Oregon State University\1. PhD\Chiller\automation\Chiller_auto_Python\ALL_complete_code_here\CFD_Final_Result_run8_fixiter50.csv',encoding='gbk')

    # 1. 定义你需要的列名列表
    input_cols = ['width', 'flowrate']
    output_cols = 'entropy'
    all_cols = input_cols+ [output_cols]

    top_10 = df[all_cols].head(15)

    # 2. 仅筛选这些列，并转换为字典
    # df[target_columns] 会创建一个只包含这些列的新 DataFrame
    # data_dict = df[target_columns].to_dict(orient='records')

    print(top_10)

    debug_bounds = [(0.2, 2.0), (0.001, 1.0)]

    print("🚀 开始测试 runBO ...")
    next_point = runBO(top_10, input_cols, output_cols, bounds=debug_bounds)
    print(next_point)