マテリアルズインフォマティクス に関する企業の技術レポート・解析内容を再現する。
はじめに
技術レポート内容
以下は、産総研の研究成果記事(2018/11/27)1より
昭和電工株式会社（社長：森川 宏平 以下、昭和電工）と国立研究開発法人 産業技術総合研究所（理事長：中鉢 良治 以下、産総研）と先端素材高速開発技術研究組合（理事長：腰塚 國博 以下、ADMAT）は、人工知能（AI）の活用により、要求特性を満たすポリマーを設計する際の試行回数を約1/40に低減できることを見いだしました。
具体的には、構造とガラス転移点が判明しているポリマーデータ 417種を使用し、そのうち、無作為に抽出した10件のデータをAIに学習させ、残りの407件の中から最もガラス転移点が高いであろうポリマーをベイズ最適化を用いて予測する。次に、選ばれたポリマーのガラス転移点情報を取得・学習データに加え、再度ガラス転移点が最も高いであろうポリマーを予測する。このサイクルを繰り返すことで、実際に所望のポリマーを発見するまでの試行回数を調べたとのこと。
この投稿では、成果発表資料2により詳細な内容が記載されているので、こちらを参考に解析内容を再現してみる。
(ポリマーのデータがないので、1290化合物の水溶解度データを使用する）
ベイズ最適化について3
Gaussian Process（ガウス過程）によって得られる予測値とその分散（予測の不確かさ）を用いて、次の実験点を決定することで、効率よく探索（目的変数の最大化）を行う手法。（ガウス過程 × 実験計画法のイメージ）
*データ駆動型の材料開発における課題
材料開発においては、目的物性の最大化（最小化）を目指す場合が多く、特に既存データの範疇よりも大きい値や、複数物性のトレードオフの解消を目的とする場合が多く、この場合はデータ適用範囲外（外挿領域）における物性予測が必要となる。
線形回帰が可能な場合は、目的変数に対する寄与の大きい変数が最大（最小）となる条件（化合物）を次の実験点とすれば良いが、非線形回帰を使用する場合には、データ適用範囲外の予測は基本的にはできない。
（例えば決定木ベースの手法においては、データ適用範囲外のサンプルの予測値は、そのサンプルと類似している既存データサンプルの平均値をとなる）
この問題はデータ駆動型の研究開発（マテリアルズインフォマティクス：MI）における共通の課題であり、その解決方法として、検証（実験）を伴わない強化学習・転移学習や、予測と実験を繰り返すことで未知領域を探索するベイズ最適化が活用されている。
目的
- ベイズ最適化により、実験（検討）回数をどれほど低減できるか検証する。
- 1290化合物のうち最も水溶解度が低い（logSが低い）化合物を予測する。
諸条件
データ
1290分子の水溶解度（LogS）データ（sdfファイル形式）
（金子研究室HPのデータセットから取得）
変数作成
- RDkit 記述子を使用（参考HP）
前処理
- 以下の変数を削除
- 分散が0の変数
- 互いに相関係数が高い（0.95以上）変数ペアの内、片方の変数
（金子研モジュール使用） - 90%以上の分子で同じ値をとる変数
- 変数は全て標準化（平均：0, 分散：1）
変数選択
ガウス過程で解析する前に、Boruta4*により変数をスクリーニングする。
（Boruta ランダムフォレストの変数重要度に基づく変数選択手法）
変数選択を行うのは、変数が多いとガウス過程を活用した際の探索空間が次元の累乗で増加＆スパースになり、実質的にランダムサーチと変わらなくなってしまうため。
解析手法
- 予測器：ガウス過程（Gaussian Process）
- ハイパーパラメータ最適化：11種類のカーネルから最適なものクロスバリデーションに決定
- 獲得関数： EI (Expected Improvement：既存データにおける目的変数の最大値（最初値）を上回る幅（更新幅）の期待値。期待改善量ともいう）
(獲得関数：次のサンプルを選択するための指標。EI 以外にも様々な獲得関数がある）
解析の流れ
logSの値が1290化合物の中で最も低い化合物(logS=-11.62)が、次の実験候補（化合物）として選ばれるまで、以下の 2~4 を繰り返す。
- logSが0以上のサンプル郡の中から初期サンプル(20サンプル）をランダムにピックアップし、既知データとし、それ以外を未知データとした。
- 既知データの前処理と変数選択
1.1 変数の削除（分散0の変数など）
1.2 標準化
1.3 Borutaによる変数選択 - 学習と予測
2.1 クロスバリデーションによるカーネル最適化
2.2 構築したモデルで、既知データを使って学習
2.3 未知データのlogSの予測とその分散を算出 - 次の実験候補の選出
3.1 期待改善量(EI)の計算
3.2 EI が最も高いサンプルの logSを取得
3.3 (取得したlogSが1290化合物の中でもっとも低い値であれば検討終了)
3.4 そうでなければ、既存データにlogSを取得したサンプルを加える。
コード
前準備）モジュールのインストール
conda install -c conda-forge rdkit
pip install boruta
pip install dcekitモジュールインポート
import os
import pandas as pd
import numpy as np
import math
import random
from scipy.stats import norm
from scipy.spatial.distance import cdist
import matplotlib.figure as figure
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.ML.Descriptors import MoleculeDescriptors
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor as rf
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, DotProduct, WhiteKernel, RBF, ConstantKernel
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error
from sklearn.ensemble import RandomForestRegressor
from sklearn.manifold import TSNE
from boruta import BorutaPy
from dcekit.variable_selection import search_high_rate_of_same_values, search_highly_correlated_variables
np.random.seed(0)
データセットの準備
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols = [mol for mol in suppl if mol is not None]
#SDFファイル中のプロパティの取得
property_list = list(mols[0].GetPropNames())
print(property_list)
#['CAS_Number', 'logS']
#RDkit記述子の計算
descriptor_names = [descriptor_name[0] for descriptor_name in Chem.Descriptors._descList]
descriptor_calculation = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
RDkit = [descriptor_calculation.CalcDescriptors(mol) for mol in mols]
#データフレーム作成
df_ = pd.DataFrame(RDkit, columns = descriptor_names)
#溶解度(logS)の値を取得
df_['logS'] = [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols]
df_ = df_.astype(float)ちなみに、各サンプルのlogSの分布と、t-SNEで二次元空間にマッピングするとこんな感じ
t-SNEによる可視化
（色はlogSを表している。赤枠のサンプルがlogSが最も低い目標サンプル）
いざ検証
i = 1
while df_known['logS'].min() != df_['logS'].min():
print('RUN',i)
df = df_known.copy()
df_X = df.drop(['logS'],axis=1)
df_y = df.loc[:,['logS']]
#欠損のある変数を削除
df_X.dropna(axis=1, inplace=True)
#分散が０の変数削除
del_num1 = np.where(df_X.var() == 0)
df_X.drop(df_X.columns[del_num1], axis=1,inplace=True)
#変数選択（互いに相関関係にある変数の内一方を削除）
corr_var = search_highly_correlated_variables(df, 0.95)
df_X.drop(df_X.columns[corr_var], axis=1, inplace=True)
#同じ値を多くもつ変数の削除
inner_fold_number = 10 #CVでの分割数（予定）
rate_of_same_value = []
threshold_of_rate_of_same_value = (inner_fold_number-1)/inner_fold_number
for col in df_X.columns:
same_value_number = df_X[col].value_counts()
rate_of_same_value.append(float(same_value_number[same_value_number.index[0]] / df_X.shape[0]))
del_var_num = np.where(np.array(rate_of_same_value) >= threshold_of_rate_of_same_value)
df_X.drop(df_X.columns[del_var_num], axis=1, inplace=True)
print('num_featrues(before boruta):', df_X.shape[1])
#borutaによる変数選択
feat_selector = BorutaPy(RandomForestRegressor(), n_estimators='auto', verbose=0,
alpha=0.05, max_iter=100, perc=80, random_state=0)
feat_selector.fit(df_X.values, df_y.values)
selected_features = df_X.iloc[:,feat_selector.support_].columns
print('num_featrues(after_boruta):',len(selected_features))
df_X = df_X.loc[:, selected_features]
X = df_X.values
y = df_y.values
#標準化（オートスケーリング）
X_scaled = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
y_scaled = (y - y.mean(axis=0)) / y.std(axis=0, ddof=1)
#ガウシアンプロセスで解析
kernels = [ConstantKernel() * DotProduct() + WhiteKernel(),
ConstantKernel() * RBF() + WhiteKernel(),
ConstantKernel() * RBF() + WhiteKernel() + DotProduct(),
ConstantKernel() * RBF(np.ones(X.shape[1])) + WhiteKernel(),
ConstantKernel() * RBF(np.ones(X.shape[1])) + WhiteKernel() + DotProduct(),
ConstantKernel() * Matern(nu=1.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=1.5) + WhiteKernel() + DotProduct(),
ConstantKernel() * Matern(nu=0.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=0.5) + WhiteKernel() + DotProduct(),
ConstantKernel() * Matern(nu=2.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=2.5) + WhiteKernel() + DotProduct()]
# CV による ε の最適化
model_in_cv = GridSearchCV(GaussianProcessRegressor(),
{'kernel': kernels}, cv=10, verbose=0, n_jobs=-1)
model_in_cv.fit(X_scaled, y_scaled)
print('score:',round(model_in_cv.best_score_,3))
model = GaussianProcessRegressor(kernel=model_in_cv.best_params_['kernel'])
model.fit(X_scaled, y_scaled)
X_unknown = df_unknown.loc[:, df_X.columns]
X_unknown_scaled = (X_unknown - df_X.mean(axis=0)) / df_X.std(axis=0, ddof=1)
pred_y_mean_, pred_y_std_ = model.predict(X_unknown_scaled, return_std=True)
pred_y_mean = pred_y_mean_.flatten()
pred_y_std = pred_y_std_.flatten()
#EI を基準に次のサンプルを選ぶ場合
pred_y_mean_i = pred_y_mean * -1
y_scaled_i = y_scaled * -1
EI_values = (pred_y_mean_i - max(y_scaled_i) - 0.01) * \
norm.cdf((pred_y_mean_i - max(y_scaled_i) - 0.01) /pred_y_std) + \
pred_y_std * norm.pdf((pred_y_mean_i - max(y_scaled_i) - 0.01) / pred_y_std)
next_sample = df_unknown.iloc[np.where(EI_values == max(EI_values))[0],:]
#単純に予測の一番高いサンプルを選ぶ場合
#next_sample = df_unknown.iloc[np.where(pred_y_mean == min(pred_y_mean))[0],:]
#ランダムに次のサンプルを選ぶ場合
#next_sample = df_unknown.loc[[random.choice(df_unknown.index.tolist())],:]
#ピックアップした化合物のlogS
print('結果:logS =',next_sample['logS'].values[0])
#調査したサンプルを加えて再度学習
df_known = pd.concat([df_known, next_sample], axis=0)
df_unknown.drop(next_sample.index, axis=0, inplace=True)
i += 1
print('====================')結果の可視化
描画t-SNE を使ってプロットを二次元に写す
df_tsne = df_.copy()
df_tsne_X = df_tsne.drop(['logS'],axis=1)
df_tsne_y = df_tsne.loc[:,['logS']]
df_tsne_X.dropna(axis=1, inplace=True)
del_num1 = np.where(df_tsne_X.var() == 0)
df_tsne_X.drop(df_tsne_X.columns[del_num1], axis=1, inplace=True)
corr_var = search_highly_correlated_variables(df_tsne_X, 0.95)
df_tsne_X.drop(df_tsne_X.columns[corr_var], axis=1, inplace=True)
X_tsne_scaled = (df_tsne_X - df_tsne_X.mean(axis=0)) / df_tsne_X.std(axis=0, ddof=1)
y_tsne_scaled = (df_tsne_y - df_tsne_y.mean(axis=0)) / df_tsne_y.std(axis=0, ddof=1)
Z = TSNE(perplexity=30, n_components=2,
init='pca',
random_state=0, n_jobs=-1).fit_transform(X_tsne_scaled)
df_plot = pd.DataFrame(Z, columns=['z1','z2'], index=df_tsne.index)予測と検証のサイクルを可視化してみる
known_samples = df_known.index[:20].tolist()
trial_samples = df_known.index[20:].tolist()
# plot after tSNE
fig = plt.figure(figsize=(13, 6))
plt.subplots_adjust(wspace=0.4, hspace=0.6)
ax1 = fig.add_subplot(1, 2, 1)
ax1.scatter(df_plot.loc[:,'z1'], df_plot.loc[:,'z2'], alpha=0.05, s = 50,
c=df_tsne_y.values, vmin=df_tsne_y.min(), vmax=df_tsne_y.max())
ax1.scatter(df_plot.loc[known_samples,'z1'], df_plot.loc[known_samples,'z2'], alpha=0.95, s = 50,
c=df_tsne_y.loc[known_samples,'logS'].values, vmin=df_tsne_y.min(), vmax=df_tsne_y.max())
ax1.set_xlabel("z1")
ax1.set_ylabel("z2")
ax1.set_title('first 20 samples', fontsize=20)
ax2 = fig.add_subplot(1, 2, 2)
ax2.plot(0, df_known.loc[known_samples,'logS'].min(),marker='o')
ax2.set_xlim(-1,700)
ax2.set_ylim(2,-12)
ax2.set_xlabel('number of trials', fontsize=20)
ax2.set_ylabel('logS(MIN) in dataset', fontsize=20)
fig.savefig('pic_randomsample/trial_0.png')
fig.show()
min_logS_list = [df_known.loc[known_samples,'logS'].min()]
for i in range(len(trial_samples)):
known_samples.append(trial_samples[i])
min_logS_list.append(df_known.loc[known_samples,'logS'].min())
fig = plt.figure(figsize=(13, 6))
plt.subplots_adjust(wspace=0.4, hspace=0.6)
ax1 = fig.add_subplot(1, 2, 1)
ax1.scatter(df_plot.loc[:,'z1'], df_plot.loc[:,'z2'], alpha=0.05,
c=df_tsne_y.values, vmin=df_tsne_y.min(), vmax=df_tsne_y.max())
ax1.scatter(df_plot.loc[known_samples,'z1'], df_plot.loc[known_samples,'z2'], alpha=0.95, s = 50,
c=df_tsne_y.loc[known_samples,'logS'].values, vmin=df_tsne_y.min(), vmax=df_tsne_y.max())
ax1.scatter(df_plot.loc[trial_samples[i],'z1'], df_plot.loc[trial_samples[i],'z2'], alpha=0.95, s = 70,
edgecolors = 'red',
c=df_tsne_y.loc[trial_samples[i],'logS'], vmin=df_tsne_y.min(), vmax=df_tsne_y.max())
ax1.set_xlabel("z1")
ax1.set_ylabel("z2")
ax1.set_title('trial {}'.format(i+1), fontsize=20)
ax2 = fig.add_subplot(1, 2, 2)
ax2.plot([j for j in range(0, i+2)], min_logS_list, marker='o')
ax2.set_xlim(-1,700)
ax2.set_ylim(2,-12)
ax2.set_xlabel('number of trials', fontsize=20)
ax2.set_ylabel('logS(MIN) in dataset', fontsize=20)
fig.savefig('pic_randomsample/trial_{}.png'.format(i+1))
fig.show()結果
ランダムに次の化合物を指定した場合
約700サイクル回して、やっと目的の化合物が選ばれる。
単純に予測値（の平均）が高い化合物を指定した場合
20サイクルまで行かずとも、目的化合物にたどり着ける。わりと早いが、一度logSが低いサンプルが選ばれると、そのサンプルの周辺を検討してしまっており、その範囲に解がなければ無駄な検討をしてしまっている。
（EIが最も高い化合物を指定した場合（ベイズ最適化）
たった5サイクルで目的化合物にたどり着けた。単純に予測値を参考にするよりも圧倒的に早い。
追加検証
初期サンプルの選び方によって、目的の化合物にたどり着くまでの試行回数は変わってくると思うので、100パターンの初期サンプルを用意して同様の検証を行って、それぞれの条件での平均試行回数を算出した。
結果（各条件での試行回数）
| 平均回数 | 最多回数 | 最小回数 | |
|---|---|---|---|
| EIが高い化合物を選出 | 7.75 | 2 | 21 |
| 予測値が高い化合物を選出 | 14.81 | 2 | 57 |
| ランダムに化合物を選出 | 658.55 | 27 | 1267 |
最小の試行回数には差がないものの、EIを指標とした方が、単純に予測値を指標とするよりも、目標の化合物を 2倍早く見つけられるという結果となった。
ちなみに箱ひげ図を書くと下記のようになる。
箱ひげ図
クローズアップ版
参考
- 人工知能（AI）を用いてポリマー設計・検証サイクルの試行回数を大幅低減（産総研記事）
- AI解析による熱硬化性ポリマー合成の思考回数削減(技術資料）
- ベイズ最適化(Bayesian Optimization, BO)～実験計画法で使ったり、ハイパーパラメータを最適化したり～（金子研究室HP)
- Boruta、ランダムフォレストの変数重要度に基づく変数選択手法
関連書籍
リンク
リンク
リンク