二日目 学習モデルとLogLikelihood関数について

前回のおさらい

dango0dan5.hatenablog.com

XenonPyを使って記述子を設定する方法を学んだ。記述子にECFP6を選んだ。

run tools.ipynb

from xenonpy.descriptor import Fingerprints

# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
count = [len(x) for x in data['SMILES']]
y,x,_=plt.hist(count, bins=20)
y=list(y)
x=list(x)
print("Peak position of SMILES length histogram",x[y.index(max(y))+1])

# Convert energy unit
data["E"]=data["E"]*0.0104

# check target properties: E & HOMO-LUMO gap
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(data['E'],data['HOMO-LUMO gap'],s=15,c='b',alpha = 0.1)
_ = plt.title('iQSPR sample data')
_ = plt.xlabel('Internal energy')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.show()

# Descriptor
ECFP6_FPs = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=128*4,counting=True, input_type='smiles')
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
tmpFP = ECFP6_FPs.transform(data['SMILES'])

学習モデルの作成

二日目の最終ゴールは、異なる学習方法で作成したモデルから、LogLikelihood関数を作成することである。

まずは、物性名を辞書として使うために、リストを先に生成しておく必要がある。

# property name will be used as a reference for calling models
prop = ['E','HOMO-LUMO gap']

コード使い回し時に、めんどくさければ、dataから"SMILES"を取り除くようにすることでも設定できる。

print(data.columns)
prop=list(data.columns)
prop.remove("SMILES")
print(prop)

Index(['SMILES', 'E', 'HOMO-LUMO gap'], dtype='object')
['E', 'HOMO-LUMO gap']

XenonPyの交差検証用にデータを分割する方法を利用して以下ではモデリングをしていく。一つ一つの物性値に対してまずは実施して、やり方を確かめてみる。k_foldの値は分割数を表す。

from xenonpy.datatools import Splitter
# prepare indices for cross-validation data sets
sp = Splitter(data.shape[0], test_size=0, k_fold=10)

次に、学習データ(train)に対応する教師データ・予測データ、テストデータ(test)に対する教師データ・予測データを格納する配列を定義する。
また、mdlsは辞書ファイルで、mdl["物性名"]とするとモデルオブジェクトが返されるものと定義する。

# initialize output variables
y_trues_test, y_preds_test = [[] for i in range(len(prop))], [[] for i in range(len(prop))]
y_trues_train, y_preds_train = [[] for i in range(len(prop))], [[] for i in range(len(prop))]
mdls = {key: [] for key in prop}

Sampleコードを参照しながら、prop["E"]に対するモデリングを行う。とりあえず、ECFP6_FPsとElasticNetの組み合わせをここでは使ってみた。

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import ElasticNet

Target="E"
Desc=ECFP6_FPs
# cross-validation test
for iTr, iTe in sp.cv():

    x_train = data['SMILES'].iloc[iTr]
    x_test = data['SMILES'].iloc[iTe]
    
    fps_train = Desc.transform(x_train)
    fps_test = Desc.transform(x_test)
    
    y_train = data[prop].iloc[iTr]
    y_test = data[prop].iloc[iTe]

    propid=prop.index(Target)

    #mdl = RandomForestRegressor(500, n_jobs=-1, max_features='sqrt')
    #mdl = ElasticNetCV(cv=5)
    #mdl = GradientBoostingRegressor(loss='absolute_error')
    mdl = ElasticNet(alpha=0.5)
        
    mdl.fit(fps_train, y_train.iloc[:,propid])
    prd_train = mdl.predict(fps_train)
    prd_test = mdl.predict(fps_test)

    y_trues_test[propid].append(y_test.iloc[:,propid].values)
    y_trues_train[propid].append(y_train.iloc[:,propid].values)
    
    y_preds_test[propid].append(prd_test)
    y_preds_train[propid].append(prd_train)
    
    mdls[Target].append(mdl)

from sklearn.metrics import r2_score

y_true_test = np.concatenate(y_trues_test[propid])
y_pred_test = np.concatenate(y_preds_test[propid])
y_true_train = np.concatenate(y_trues_train[propid])
y_pred_train = np.concatenate(y_preds_train[propid])
xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
xy_diff = xy_max - xy_min
r2=r2_score(y_true_test,y_pred_test)
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
_ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
_ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
_ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
_ = plt.title('Property: ' + Target)
_ = plt.xlim(xy_min,xy_max)
_ = plt.ylim(xy_min,xy_max)
_ = plt.legend(loc='upper left')
_ = plt.xlabel('Prediction')
_ = plt.ylabel('Observation')
_ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
_ = plt.show()

sampleの予測結果とあまり変わらないように見えている。ここで、mdlsには、各検証でFittingした予測モデルが入っている。この予測モデルを調べてみる。一つのSMILESに対して、記述子を計算し、入れてみると、それぞれのモデルで値が若干異なっていることがわかる。これは、LogLikelihood関数を作成するときに大事なことになる。

tmpFP = ECFP6_FPs.transform([data["SMILES"][1]])
for x in mdls["E"]:
    print(x.predict(tmpFP))

[1.95511605]
[1.95277628]
[1.95631381]
[1.94501778]
[1.95411992]
[1.94175693]
[1.94580894]
[1.94818088]
[1.95047658]
[1.95249834]

一個一個個別に作って行ってもいいのだが、結局は同じ操作をすることになります。それぞれの物性値の予測モデルの、学習モデルと記述子を個別に与えたいために、以下の関数を導入します。deepcopyを使わなければならないのに注意。

from copy import deepcopy

# initialize model dictionary
mdls = {key: [] for key in prop}

def CrossValidation(sp,propid,Desc,model):
    y_trues_test=[]
    y_preds_test=[]
    y_trues_train=[]
    y_preds_train=[]
    out_mdls=[]
    # cross-validation test
    for iTr, iTe in sp.cv():

        x_train = data['SMILES'].iloc[iTr]
        x_test = data['SMILES'].iloc[iTe]
    
        fps_train = Desc.transform(x_train)
        fps_test = Desc.transform(x_test)
    
        y_train = data[prop].iloc[iTr]
        y_test = data[prop].iloc[iTe]

        mdl = deepcopy( model )
        
        mdl.fit(fps_train, y_train.iloc[:,propid])
        prd_train = mdl.predict(fps_train)
        prd_test = mdl.predict(fps_test)

        y_trues_test.append(y_test.iloc[:,propid].values)
        y_trues_train.append(y_train.iloc[:,propid].values)
    
        y_preds_test.append(prd_test)
        y_preds_train.append(prd_train)
        out_mdls.append(mdl)
    
    y_true_test = np.concatenate(y_trues_test)
    y_pred_test = np.concatenate(y_preds_test)
    y_true_train = np.concatenate(y_trues_train)
    y_pred_train = np.concatenate(y_preds_train)
    xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
    xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
    xy_diff = xy_max - xy_min

    r2=r2_score(y_true_test,y_pred_test)

    _ = plt.figure(figsize=(5,5))
    _ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
    _ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
    _ = plt.title('Property: ' + Target)
    _ = plt.xlim(xy_min,xy_max)
    _ = plt.ylim(xy_min,xy_max)
    _ = plt.legend(loc='upper left')
    _ = plt.xlabel('Prediction')
    _ = plt.ylabel('Observation')
    _ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
    _ = plt.show()
    return out_mdls

n_bitsを小さくして、予測モデル作成を早く終わるようにしてみます。16の値を入れたことはないのだが、やってみました。

Target="E"
ECFP6_FPs = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
Desc=ECFP6_FPs
model=GradientBoostingRegressor(loss='squared_error')
mdls[Target]=CrossValidation(sp,prop.index(Target),Desc,model)

結果として、良いR2の値と、MAEの値となっていましたね。次に、予測モデルを使って、一個のSMILESから予想される値を調べると、バラついていることがわかります。

tmpFP = ECFP6_FPs.transform([data["SMILES"][1]])
for x in mdls["E"]:
    print(x.predict(tmpFP))

[1.55619993]
[1.4899335]
[1.57128927]
[1.58779223]
[1.56506169]
[1.51193849]
[1.61289896]
[1.54215519]
[1.62158862]
[1.52752974]

次に、バンドギャップについても作成してみますが、違う記述子と違うモデルを作成してみました。

from sklearn.ensemble import RandomForestRegressor

Target="HOMO-LUMO gap"
Desc= Fingerprints(featurizers=['PatternFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=RandomForestRegressor(10, n_jobs=-1, max_features='sqrt')
mdls[Target]=CrossValidation(sp,prop.index(Target),Desc,model)

これもまぁまぁなモデルが作れています。Sampleと違って、10eV以上の点で改善されているように見えます。RadomForestがベストなものではないので、いろいろ変えてみると、さらに良くなるかもしれません。

tmpFP = Desc.transform([data["SMILES"][1]])
for x in mdls["HOMO-LUMO gap"]:
    print(x.predict(tmpFP))

LogLikelihood関数

LogLikelihood関数を定義します。ここでは、少し特殊なことをしています。custom_mdlsを定義するところまでは、sampleと同じです。今回は別々の記述子を使うので、それぞれの予測モデルを組み込んだcustom_mdls1とcustom_mdls2を定義しています。それぞれの記述子を指定するために、GaussianLogLikelihood関数をつかって、prd_mdls1とprd_mdls2を定義しています。

# Forward model template in XenonPy-iQSPR 
from xenonpy.inverse.iqspr import GaussianLogLikelihood

class bootstrap_fcn():
    def __init__(self, ms, var):
        self.Ms = ms
        self.Var = var
    def predict(self, x):
        val = np.array([m.predict(x) for m in self.Ms])
        return np.mean(val, axis = 0), np.sqrt(np.var(val, axis = 0) + self.Var)
    
# include basic measurement noise
c_var = {'E': 0.1, 'HOMO-LUMO gap': 0.1}

# generate a dictionary for the final models used to create the likelihood
custom_mdls = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() }
custom_mdls1 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="E"}
custom_mdls2 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="HOMO-LUMO gap"}

# import descriptor calculator and forward model to iQSPR
prd_mdls1 = GaussianLogLikelihood(descriptor=ECFP6_FPs, targets={'E': (0, 2)}, **custom_mdls1)
prd_mdls2 = GaussianLogLikelihood(descriptor=Desc, targets={'HOMO-LUMO gap': (-np.inf, 3)}, **custom_mdls2)

このprd_mdls1が予測モデルとして使えるかどうか、調べてみましょう。

prd_mdls1.predict(data["SMILES"])

次からは難問です。オリジナルのLogLikelihoodを定義するのですが、predictorの設定や、それに伴ってsampleにあったコードを書き換える必要があります。基本的には、GaussianLogLikelihood関数のソースコードを参考にしながら、改変することでうまくいきました。
github.com

from xenonpy.inverse.base import BaseLogLikelihoodSet
from sklearn.linear_model import BayesianRidge

target_range = {'E': (0,2), 'HOMO-LUMO gap': (-np.inf, 3)}
desc={"E":ECFP6_FPs,'HOMO-LUMO gap':Desc}

class MyLogLikelihood(BaseLogLikelihoodSet):
    def __init__(self):
        super().__init__()
        self.loglike = prd_mdls1
        self.loglike = prd_mdls2
        self._mdl=custom_mdls
        self.descriptor=desc
        self.targets=target_range
    def predict(self, smiles, **kwargs):
        tmp = {}
        for k, v in self._mdl.items():
            fps = self.descriptor[k].transform(smiles, return_type='df')
            fps_ = fps.dropna()
            tmp[k + ': mean'], tmp[k + ': std'] = v.predict(fps_)
        tmp = pd.DataFrame(data=tmp, index=fps_.index)
        return pd.DataFrame(data=tmp, index=fps.index)

like_mdl = MyLogLikelihood()

これでpredictorが動くかどうか調べてみましたが、うまく行っていますね。

pred=like_mdl.predict(data['SMILES'])
pred

では、最後にいろいろとプロットさせていきます。

# calculate log-likelihood for a given target property region
tmp_ll = like_mdl(data['SMILES'])
tmp = tmp_ll.sum(axis = 1, skipna = True)

_ = plt.figure(figsize=(8,5))
_ = plt.hist(tmp, bins=50)
_ = plt.title('ln(likelihood) for SMILES in the data subset')
_ = plt.xlabel('ln(likelihood)')
_ = plt.ylabel('Counts')
_ = plt.show()

# plot histogram of likelihood values
_ = plt.figure(figsize=(8,5))
_ = plt.hist(np.exp(tmp), bins=50)
_ = plt.title('Likelihood for SMILES in the data subset')
_ = plt.xlabel('Likelihood')
_ = plt.ylabel('Counts')
_ = plt.show()

# check the predicted likelihood
dot_scale = 0.5
l_std = np.sqrt(pred['E: std']**2+pred['HOMO-LUMO gap: std']**2)

_ = plt.figure(figsize=(6,5))
rectangle = plt.Rectangle((0,0),2,3,fc='y',alpha=0.1)
_ = plt.gca().add_patch(rectangle)
im = plt.scatter(pred['E: mean'], pred['HOMO-LUMO gap: mean'], s=l_std*dot_scale, c=np.exp(tmp),alpha = 0.2,cmap=plt.get_cmap('cool'))
_ = plt.title('Pubchem data')
_ = plt.xlabel('E')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.colorbar(im)
_ = plt.show()

ここで作成したモデルは、高いバンドギャップの領域(10eV以上）での予想の傾向が向上していたので、そちらの探索にも使えるかもしれませんね。次回以降は、領域指定の値を変えて実施してもいいかなと思います。

二日目のまとめ

これらの結果をまとめると、下のようなコードになります。jupyter-notebookでやるときには適当なところで、セルを分割してください。

run tools.ipynb

from xenonpy.descriptor import Fingerprints

# load in-house data
data = pd.read_csv("./iQSPR_sample_data.csv", index_col=0)
data=data[ data["SMILES"] != "FBr(F)(F)(F)F"  ]
# make histogram
count = [len(x) for x in data['SMILES']]
y,x,_=plt.hist(count, bins=20)
y=list(y)
x=list(x)
print("Peak position of SMILES length histogram",x[y.index(max(y))+1])

# Convert energy unit
data["E"]=data["E"]*0.0104

# check target properties: E & HOMO-LUMO gap
_ = plt.figure(figsize=(5,5))
_ = plt.scatter(data['E'],data['HOMO-LUMO gap'],s=15,c='b',alpha = 0.1)
_ = plt.title('iQSPR sample data')
_ = plt.xlabel('Internal energy')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.show()

# set target range
target_range = {'E': (0,2), 'HOMO-LUMO gap': (-np.inf, 3)}

# prepare indices for cross-validation data sets
from xenonpy.datatools import Splitter
sp = Splitter(data.shape[0], test_size=0, k_fold=10)

# set property list
prop=list(data.columns)
prop.remove("SMILES")

from copy import deepcopy
from sklearn.metrics import r2_score

# initialize model dictionary
mdls = {key: [] for key in prop}

def CrossValidation(sp,propid,Desc,model):
    y_trues_test=[]
    y_preds_test=[]
    y_trues_train=[]
    y_preds_train=[]
    out_mdls=[]
    # cross-validation test
    for iTr, iTe in sp.cv():

        x_train = data['SMILES'].iloc[iTr]
        x_test = data['SMILES'].iloc[iTe]
    
        fps_train = Desc.transform(x_train)
        fps_test = Desc.transform(x_test)
    
        y_train = data[prop].iloc[iTr]
        y_test = data[prop].iloc[iTe]

        mdl = deepcopy( model )
        
        mdl.fit(fps_train, y_train.iloc[:,propid])
        prd_train = mdl.predict(fps_train)
        prd_test = mdl.predict(fps_test)

        y_trues_test.append(y_test.iloc[:,propid].values)
        y_trues_train.append(y_train.iloc[:,propid].values)
    
        y_preds_test.append(prd_test)
        y_preds_train.append(prd_train)
        out_mdls.append(mdl)
    
    y_true_test = np.concatenate(y_trues_test)
    y_pred_test = np.concatenate(y_preds_test)
    y_true_train = np.concatenate(y_trues_train)
    y_pred_train = np.concatenate(y_preds_train)
    xy_min = min(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*0.95
    xy_max = max(np.concatenate([y_true_test,y_true_train,y_pred_test,y_pred_train]))*1.05
    xy_diff = xy_max - xy_min

    r2=r2_score(y_true_test,y_pred_test)

    _ = plt.figure(figsize=(5,5))
    _ = plt.scatter(y_pred_train, y_true_train, c='b', alpha=0.1, label='train')
    _ = plt.scatter(y_pred_test, y_true_test, c='r', alpha=0.2, label='test')
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.10,' R2: %5.2f' % r2,fontsize=12)
    _ = plt.text(xy_min+xy_diff*0.7,xy_min+xy_diff*0.05,'MAE: %5.2f' % np.mean(np.abs(y_true_test - y_pred_test)),fontsize=12)
    _ = plt.title('Property: ' + Target)
    _ = plt.xlim(xy_min,xy_max)
    _ = plt.ylim(xy_min,xy_max)
    _ = plt.legend(loc='upper left')
    _ = plt.xlabel('Prediction')
    _ = plt.ylabel('Observation')
    _ = plt.plot([xy_min,xy_max],[xy_min,xy_max],ls="--",c='k')
    _ = plt.show()
    return out_mdls

from sklearn.ensemble import GradientBoostingRegressor
Target="E"
FP_E = Fingerprints(featurizers=['ECFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=GradientBoostingRegressor(loss='squared_error')
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_E,model)

from sklearn.ensemble import RandomForestRegressor
Target="HOMO-LUMO gap"
FP_H= Fingerprints(featurizers=['PatternFP',"MACCS"],radius=3,n_bits=16,counting=True, input_type='smiles')
model=RandomForestRegressor(10, n_jobs=-1, max_features='sqrt')
mdls[Target]=CrossValidation(sp,prop.index(Target),FP_H,model)


# Forward model template in XenonPy-iQSPR 
from xenonpy.inverse.iqspr import GaussianLogLikelihood

class bootstrap_fcn():
    def __init__(self, ms, var):
        self.Ms = ms
        self.Var = var
    def predict(self, x):
        val = np.array([m.predict(x) for m in self.Ms])
        return np.mean(val, axis = 0), np.sqrt(np.var(val, axis = 0) + self.Var)
    
# include basic measurement noise
c_var = {'E': 0.1, 'HOMO-LUMO gap': 0.1}

# generate a dictionary for the final models used to create the likelihood
custom_mdls1 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="E"}
custom_mdls2 = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() if key=="HOMO-LUMO gap"}

# import descriptor calculator and forward model to iQSPR
prd_mdls1 = GaussianLogLikelihood(descriptor=FP_E, targets={'E': target_range["E"]}, **custom_mdls1)
prd_mdls2 = GaussianLogLikelihood(descriptor=FP_H, targets={'HOMO-LUMO gap': target_range['HOMO-LUMO gap']}, **custom_mdls2)


from xenonpy.inverse.base import BaseLogLikelihoodSet
from sklearn.linear_model import BayesianRidge

custom_mdls = {key: bootstrap_fcn(value, c_var[key])  for key, value in mdls.items() }
desc={"E":FP_E,'HOMO-LUMO gap':FP_H}

class MyLogLikelihood(BaseLogLikelihoodSet):
    def __init__(self):
        super().__init__()
        self.loglike = prd_mdls1
        self.loglike = prd_mdls2
        self._mdl=custom_mdls
        self.descriptor=desc
        self.targets=target_range
    def predict(self, smiles, **kwargs):
        tmp = {}
        for k, v in self._mdl.items():
            fps = self.descriptor[k].transform(smiles, return_type='df')
            fps_ = fps.dropna()
            tmp[k + ': mean'], tmp[k + ': std'] = v.predict(fps_)
        tmp = pd.DataFrame(data=tmp, index=fps_.index)
        return pd.DataFrame(data=tmp, index=fps.index)

like_mdl = MyLogLikelihood()

pred=like_mdl.predict(data['SMILES'])

# calculate log-likelihood for a given target property region
tmp_ll = like_mdl(data['SMILES'])
tmp = tmp_ll.sum(axis = 1, skipna = True)

_ = plt.figure(figsize=(8,5))
_ = plt.hist(tmp, bins=50)
_ = plt.title('ln(likelihood) for SMILES in the data subset')
_ = plt.xlabel('ln(likelihood)')
_ = plt.ylabel('Counts')
_ = plt.show()

# plot histogram of likelihood values
_ = plt.figure(figsize=(8,5))
_ = plt.hist(np.exp(tmp), bins=50)
_ = plt.title('Likelihood for SMILES in the data subset')
_ = plt.xlabel('Likelihood')
_ = plt.ylabel('Counts')
_ = plt.show()

# check the predicted likelihood
dot_scale = 0.5
l_std = np.sqrt(pred['E: std']**2+pred['HOMO-LUMO gap: std']**2)

_ = plt.figure(figsize=(6,5))
rectangle = plt.Rectangle((0,0),2,3,fc='y',alpha=0.1)
_ = plt.gca().add_patch(rectangle)
im = plt.scatter(pred['E: mean'], pred['HOMO-LUMO gap: mean'], s=l_std*dot_scale, c=np.exp(tmp),alpha = 0.2,cmap=plt.get_cmap('cool'))
_ = plt.title('Pubchem data')
_ = plt.xlabel('E')
_ = plt.ylabel('HOMO-LUMO gap')
_ = plt.colorbar(im)
_ = plt.show()