GPyOptで化合物探索プログラムを実装

#データセット約2000件を学習用データ、予測用データに分割して保存。
#seed値を0〜90として10セット生成
import pandas as pd
import random

# CSVファイルの読み込み
data_x = pd.read_csv('data_x.csv', index_col='cid')
data_y = pd.read_csv('data_y.csv', index_col='cid')

# 抽出する行数が元の行数より多い場合は、元の行数に合わせる
num_rows = min(30, len(data_x))

#Seed値を設定
seed = [x for x in range(0,91) if x%10==0]

for seed in seed:
    # ランダムに行を抽出
    random_rows = data_x.sample(n=num_rows, random_state=random.seed(seed))

    # 抽出した行を別のファイルに保存する場合
    random_rows.to_csv(f'train_x_seed_{seed}.csv', index='cid')
    data_y.loc[random_rows.index].to_csv(f'train_y_seed_{seed}.csv', index='cid')

    # 元のデータから抽出した行を削除
    data_x = data_x.drop(random_rows.index)
    data_y = data_y.drop(random_rows.index)

    # 更新されたデータを別のファイルに保存する場合
    data_x.to_csv(f'prediction_data_x_seed_{seed}.csv', index='cid')
    data_y.to_csv(f'prediction_data_y_seed_{seed}.csv', index='cid')

import pickle
import GPy
import GPyOpt
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.preprocessing import StandardScaler
from scipy.stats import norm



#目的関数の設定。今回は用意した候補から選択するためboundsは使用しないが、最適化が実行できないのでダミーとして設定
def f(x):
    global trial
    
    #print(f'trial:{trial}')
    X_1st = x[0]
    y = float(input('実験値は？'))
    return y

#探索範囲の設定。今回は用意した候補から選択するためboundsは使用しないが、最適化が実行できないのでダミーとして設定
bounds = []
x.columns.tolist()

for criptor in x.columns.tolist():
    bound = {}
    bound['name'] = criptor
    bound['type'] = 'discrete'
    bound['domain'] = x_prediction[criptor].tolist()
    bounds.append(bound)

#10個のデータセットについてベイズ最適化を実行
for seed in [seed for seed in range(0,91,10)]:
    
    #試行回数を初期化
    trial=1

    # 予測用データをCSVファイルからデータを読み込む
    x_prediction = pd.read_csv(f'prediction_data_x_seed_{seed}.csv', index_col=0, header=0)
    y_answer = pd.read_csv(f'prediction_data_y_seed_{seed}.csv', index_col=0, header=0)

    # 学習用の目的変数と説明変数をCSVファイルから読み込む
    initial_y = pd.read_csv(f'train_y_seed_{seed}.csv', index_col=0, header=0)  # 目的変数
    initial_y = initial_y.iloc[:,0]
    initial_y = np.array(initial_y)

    initial_x = pd.read_csv(f'train_x_seed_{seed}.csv', index_col=0, header=0)  # 説明変数

    max_initial_y = np.max(initial_y)
    suggested_y = 0

    suggested_data = []

    #学習データの最大値を更新するか、試行回数が200回を超えるまで繰り返しベイズ最適化を行う
    while max_initial_y > float(suggested_y) and trial < 200:

        print(f'試行回数：{trial}')

        # 標準化のためのスケーラーを作成
        scaler = StandardScaler()

        #予測用データの標準化を行う
        autoscaled_x_prediction = pd.DataFrame(scaler.fit_transform(x_prediction), columns=x_prediction.columns)
        autoscaled_x_prediction = np.array(autoscaled_x_prediction)
        x_prediction_array = np.array(x_prediction)
        
        #目的関数の標準化を行う
        initial_y = initial_y.reshape(-1,1)
        mean_of_y = initial_y.mean()
        std_of_y = initial_y.std()
        autoscaled_initial_y = (initial_y - initial_y.mean()) / initial_y.std()
        autoscaled_initial_y = np.array(autoscaled_initial_y)
        autoscaled_initial_y = autoscaled_initial_y.reshape(-1,1)

        # 標準化のためのスケーラーを作成
        scaler = StandardScaler()

        #説明変数の標準化を行う
        autoscaled_initial_x = pd.DataFrame(scaler.fit_transform(initial_x), columns=initial_x.columns)
        autoscaled_initial_x = np.array(autoscaled_initial_x)

        # カーネルの選択
        kernel_choices = [GPy.kern.RBF(input_dim=initial_x.shape[1]), GPy.kern.Matern52(input_dim=initial_x.shape[1])]
        selected_kernel = kernel_choices[0]  # カーネルを選択
        
        #ベイズ最適化の条件を設定
        opt = GPyOpt.methods.BayesianOptimization(f=f, domain=bounds,
                                                 model_type='GP',
                                                 #kernel=GPy.kern.RBF(input_dim=initial_x.shape[1],ARD=True),
                                                 X=autoscaled_initial_x,
                                                 Y=autoscaled_initial_y,
                                                 initial_design_numdata=len(initial_y),
                                                  maximize = True,
                                                  normalize_Y = False,
                                                 )  # acquisition_weightを調整
        
        #既にサンプリングした点を除くためにignored_Xで入力点を指定
        x_suggest = opt.suggest_next_locations(ignored_X=autoscaled_initial_x) 
        y_predict = opt.model.model.predict(x_suggest) 
        
        #学習したGPRモデルを使用して予測
        gp_model = opt.model.model
        gp_mean = gp_model.predict(autoscaled_x_prediction)[0]
        gp_std = gp_model.predict(autoscaled_x_prediction)[1]
        current_optimum = np.max(initial_y)
        
        #予測用データのEIを計算
        EI_test = []

        xi = 0.01
        for i in range(len(gp_mean)):

            mean = gp_mean[i]
            std = gp_std[i]
            f_max = np.max(initial_y)
            z = (mean - f_max - xi) / std
            ei = (mean - f_max - xi) * norm.cdf(z) + std * norm.pdf(z)
            EI_test.append(ei)
            
        #獲得関数をグラフ化
        fig = plt.figure(figsize=(8, 6))
        plt.plot(EI_test)
        plt.xlabel('$x$', fontsize=16)
        plt.ylabel('Expected Improvement', fontsize=16)
        plt.tick_params(labelsize=16)
        plt.show()

        #獲得関数が最も大きかったデータを学習データに追加し、予測用データから削除する
        suggested = x_prediction.iloc[EI_test.index(np.max(EI_test)),:]
        initial_x = initial_x.append(suggested)
        x_prediction = x_prediction.drop(index=x_prediction.index[EI_test.index(np.max(EI_test))])

        suggested_y = y_answer.iloc[EI_test.index(np.max(EI_test)),:]
        initial_y = np.insert(initial_y, len(initial_y), suggested_y)
        initial_y = initial_y.reshape(-1,1)
        y_answer = y_answer.drop(y_answer.index[EI_test.index(np.max(EI_test))])

        print(f'選択されたxlogp : {suggested_y}\n'
        f'トレーニングデータの最大値　：　{max_initial_y}')
        
        #推奨された候補をリストへ追加
        suggested_data.append(float(suggested_y))

        # 保存するリスト
        sava_list = suggested_data

        # リストを保存するファイルパス
        file_path = f'suggested_data_seed_{seed}.pkl'

        # リストを保存
        with open(file_path, 'wb') as f:
            pickle.dump(sava_list, f)

        trial += 1