五日目 iQSPR

iQSPR

今回の検討で分かったことがある。

「SMILESの微小改変に追従するぐらい予測モデルの値が変化しないとXenonPyによる構造改変が機能しない。」

少ない記述子で十分に予想できるモデルを作成したとしても、iQSPRを実行するとターゲット領域の構造を探索しきれないことに繋がってしまう。これは記述子だけではなく、学習モデル自体に対しても要求されるはずだ。

今回は比較的手軽な機械学習モデルを用いているが、Fingerprint・機械学習モデルの組み合わせによっては、iQSPRの探索範囲が限定される結果をもたらす可能性がある。そのため、化学構造探索数・精度の観点で、Fingerprintと機械学習モデルの最適化を行うのが、実用上望まれるであろう。

16から256にn_Bitsを変更することで、標的領域に化学構造が現れるようになった。

まず最初に、作成したN-gramを読み込む。

prd_mdls=like_mdl

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# load a pre-trained n-gram from the pickle file
n_gram = joblib.load('ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'rb') as f:
#     n_gram = pk.load(f)
    
# load a pre-trained n-gram from the pickle file
n_gram2 = joblib.load('ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'rb') as f:
#     n_gram2 = pk.load(f)

_ = n_gram.merge_table(n_gram2)

次に、初期の探索点を選定する。

# set up initial molecules for iQSPR
np.random.seed(201906) # fix the random seed
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi)) for i, smi in enumerate(data['SMILES'])
       if (data['HOMO-LUMO gap'].iloc[i] > 3.5)]
init_samples = np.random.choice(cans, 25)
print(init_samples)

また、アニーリング（書籍参照のこと）パラメーターを設定する。

# set up annealing schedule in iQSPR
beta = np.hstack([np.linspace(0.01,0.2,20),np.linspace(0.21,0.4,10),np.linspace(0.4,1,10),np.linspace(1,1,10)])
print('Number of steps: %i' % len(beta))
print(beta)

ここから、N-gramとLogLikelihood関数をIQSPRクラスに設定して、探索を行う。

# library for running iQSPR in XenonPy-iQSPR
from xenonpy.inverse.iqspr import IQSPR

# update NGram parameters for this exampleHOMO-LUMO gap
n_gram.set_params(del_range=[1,20],max_len=500, reorder_prob=0.5, sample_order=(1,20))

# set up likelihood and n-gram models in iQSPR
iqspr_reorder = IQSPR(estimator=prd_mdls, modifier=n_gram)
    
# main loop of iQSPR
iqspr_samples1, iqspr_loglike1, iqspr_prob1, iqspr_freq1 = [], [], [], []
for s, ll, p, freq in iqspr_reorder(init_samples, beta, yield_lpf=True):
    iqspr_samples1.append(s)
    iqspr_loglike1.append(ll)
    iqspr_prob1.append(p)
    iqspr_freq1.append(freq)
# record all outputs
iqspr_results_reorder = {
    "samples": iqspr_samples1,
    "loglike": iqspr_loglike1,
    "prob": iqspr_prob1,
    "freq": iqspr_freq1,
    "beta": beta
}
# save results
with open('iQSPR_results_reorder.obj', 'wb') as f:
    pk.dump(iqspr_results_reorder, f)

探索する際に、likelihood関数がどのような構造が得られたのか、履歴を確認する。Target領域近傍の点がサンプルできた。

# plot the likelihood evolution

# set up the min and max boundary for the plots
tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]
flat_list = np.asarray([item for sublist in tmp_list for item in sublist])
y_max, y_min = max(flat_list), min(flat_list)

_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Log-likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), ll ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_loglike_reorder.png',dpi = 500)
#plt.close()

y_max, y_min = np.exp(y_max), np.exp(y_min)
_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), np.exp(ll) ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_like_reorder.png',dpi = 500)
#plt.close()

探索された点に対して、改めて予測モデルを適用してどの程度の予測ブレがあるのか調べる。結構大きい結果となったので、境界にいる構造は外れている可能性もある。

# re-calculate the property values for the proposed molecules
x_mean, x_std, y_mean, y_std = [], [], [], []
r_std = []
FPs_samples = []
for i, smis in enumerate(iqspr_results_reorder["samples"]):
    tmp_fps = FP_E.transform(smis)
    FPs_samples.append(tmp_fps)
    
    tmp1, tmp2 = custom_mdls1["E"].predict(tmp_fps)
    x_mean.append(tmp1)
    x_std.append(tmp2)
    
    tmp_fps1 = FP_H.transform(smis)
    tmp1, tmp2 = custom_mdls2["HOMO-LUMO gap"].predict(tmp_fps1)
    y_mean.append(tmp1)
    y_std.append(tmp2)
    
    r_std.append([np.sqrt(x_std[-1][i]**2 + y_std[-1][i]**2) for i in range(len(x_std[-1]))])

# flatten the list for max/min calculation
flat_list = [item for sublist in r_std for item in sublist]
print('Range of std. dev.: (%.4f,%.4f)' % (min(flat_list),max(flat_list)))

Range of std. dev.: (0.7748,0.8449)

探索履歴を、二つの物性値のグラフ上で動かした履歴を確認できるようにする。

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_prd/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

flat_list = np.asarray([item for sublist in r_std for item in sublist])
s_max, s_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["E"],
    np.asarray([item for sublist in x_mean for item in sublist])))
x_max, x_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["HOMO-LUMO gap"],
    np.asarray([item for sublist in y_mean for item in sublist])))
y_max, y_min = max(flat_list), min(flat_list)
tmp_beta = iqspr_results_reorder["beta"]

for i in range(len(r_std)):
    dot_size = 45*((np.asarray(r_std[i])-s_min)/(s_max-s_min)) + 5
    
    _ = plt.figure(figsize=(5,5))
    rectangle = plt.Rectangle((0,0),2.5,3,fc='y',alpha=0.1)
    _ = plt.gca().add_patch(rectangle)
    _ = plt.scatter(data["E"], data["HOMO-LUMO gap"],s=3, c='b', alpha=0.2)
    _ = plt.scatter(x_mean[i], y_mean[i],s=dot_size, c='r', alpha=0.5)
    _ = plt.title('Step: %i (beta = %.3f)' % (i,tmp_beta[i]))
    _ = plt.xlim(x_min,x_max)
    _ = plt.ylim(y_min,y_max)
    _ = plt.xlabel('Internal energy')
    _ = plt.ylabel('HOMO-LUMO gap')
    #plt.show()
    _ = plt.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

化学構造の変遷を見てみる。

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_smiles/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]    

n_S = 25
for i, smis in enumerate(iqspr_results_reorder['samples']):
    tmp_smis = iqspr_results_reorder['samples'][i][
        np.argsort(tmp_list[i])[::-1]]
    fig, ax = plt.subplots(5, 5)
    _ = fig.set_size_inches(20, 20)
    _ = fig.set_tight_layout(True)
    for j in range(n_S):
        xaxis = j // 5
        yaxis = j % 5
        try:
            img = Draw.MolToImage(Chem.MolFromSmiles(tmp_smis[j]))
            _ = ax[xaxis, yaxis].clear()
            _ = ax[xaxis, yaxis].set_frame_on(False)
            _ = ax[xaxis, yaxis].imshow(img)
        except:
            pass
        _ = ax[xaxis, yaxis].set_axis_off()
    _ = fig.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()
    

ターゲット領域の設定からして、分子量が小さくて、電子共鳴しやすい構造であるなというのが想像つくのだが、おおよそそうなっていたので、うまく行ったということだろう。

まとめ

iQSPRについては、ほとんど寄り道がなかったので、まとめは実施しない。コード自体はチュートリアルの域を出なかったが、パラメーターを色々変えてみた場合にどうなるか、実施してみることで、次の応用への理解が深まった。

github.com

"Thank you for using XenonPy-iQSPR. We would appreciate any feedback and code contribution to this open-source project.The sample code can be downloaded at: https://github.com/yoshida-lab/XenonPy/blob/master/samples/iQSPR.ipynb"

結局のところ、フィードバックを与えるような試みもしてはいなかった。

likelihoood関数とか記述子くらいか。一次元探索の問題にしたいときは、ダミーの物性予測（RDKitで算出できる分子量とか）をもう一次元に設定して、予測モデルもRDKitの関数にするなりして、範囲も(0.0,np.inf)とかにすれば良さそう。

興味を持った読者は、自身で実行しながら、githubサイト、書籍を読んではいかがだろうか。サンプルコードも寄り道をしながらわかりやすいように説明してくれているので、XenonPyで用意されたフレームワークが柔軟に設計されていることは大変ありがたい。記述子周りも、すっきりさせることができるので、表現方法探索をするのにもこのコードを応用してやれば、泥臭いコードでエラーが頻発するなんてことも提言できそうなので、仕事にも役立つかもしれない。

書籍版のiQSPR用のコードでは、学習モデルのハイパーパラメーター探索もしていたようなので、後で確認してみようと思う。

最後にコード全体がどうなっていたのかを再表示してみる。

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.5), '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=256,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',"ECFP","MACCS",'LayeredFP'],radius=10,n_bits=256,counting=True, input_type='smiles')
model=RandomForestRegressor(100, n_jobs=-1, max_features=None)
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.3, 'HOMO-LUMO gap': 0.3}

# 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.5,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()

prd_mdls=like_mdl

# N-gram library in XenonPy-iQSPR
from xenonpy.inverse.iqspr import NGram

# load a pre-trained n-gram from the pickle file
n_gram = joblib.load('ngram_reO15_O10.xz')
# with open('ngram_pubchem_ikebata_reO15_O10.obj', 'rb') as f:
#     n_gram = pk.load(f)
    
# load a pre-trained n-gram from the pickle file
n_gram2 = joblib.load('ngram_reO15_O11to20.xz')
# with open('ngram_pubchem_ikebata_reO15_O11to20.obj', 'rb') as f:
#     n_gram2 = pk.load(f)

_ = n_gram.merge_table(n_gram2)

# set up initial molecules for iQSPR
np.random.seed(201906) # fix the random seed
cans = [Chem.MolToSmiles(Chem.MolFromSmiles(smi)) for i, smi in enumerate(data['SMILES'])
       if (data['HOMO-LUMO gap'].iloc[i] > 3.5)]
init_samples = np.random.choice(cans, 25)
print(init_samples)

# set up annealing schedule in iQSPR
beta = np.hstack([np.linspace(0.01,0.2,20),np.linspace(0.21,0.4,10),np.linspace(0.4,1,10),np.linspace(1,1,10)])
print('Number of steps: %i' % len(beta))
print(beta)

# library for running iQSPR in XenonPy-iQSPR
from xenonpy.inverse.iqspr import IQSPR

# update NGram parameters for this exampleHOMO-LUMO gap
n_gram.set_params(del_range=[1,20],max_len=500, reorder_prob=0.5, sample_order=(1,20))

# set up likelihood and n-gram models in iQSPR
iqspr_reorder = IQSPR(estimator=prd_mdls, modifier=n_gram)
    
# main loop of iQSPR
iqspr_samples1, iqspr_loglike1, iqspr_prob1, iqspr_freq1 = [], [], [], []
for s, ll, p, freq in iqspr_reorder(init_samples, beta, yield_lpf=True):
    iqspr_samples1.append(s)
    iqspr_loglike1.append(ll)
    iqspr_prob1.append(p)
    iqspr_freq1.append(freq)
# record all outputs
iqspr_results_reorder = {
    "samples": iqspr_samples1,
    "loglike": iqspr_loglike1,
    "prob": iqspr_prob1,
    "freq": iqspr_freq1,
    "beta": beta
}
# save results
with open('iQSPR_results_reorder.obj', 'wb') as f:
    pk.dump(iqspr_results_reorder, f)


with open('iQSPR_results_reorder.obj', 'rb') as f:
    iqspr_results_reorder = pk.load(f)

# plot the likelihood evolution

# set up the min and max boundary for the plots
tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]
flat_list = np.asarray([item for sublist in tmp_list for item in sublist])
y_max, y_min = max(flat_list), min(flat_list)

_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Log-likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), ll ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_loglike_reorder.png',dpi = 500)
#plt.close()

y_max, y_min = np.exp(y_max), np.exp(y_min)
_ = plt.figure(figsize=(10,5))
_ = plt.xlim(0,len(iqspr_results_reorder["loglike"]))
_ = plt.ylim(y_min,y_max)
_ = plt.xlabel('Step')
_ = plt.ylabel('Likelihood')
for i, ll in enumerate(tmp_list):
    _ = plt.scatter([i]*len(ll), np.exp(ll) ,s=10, c='b', alpha=0.2)
_ = plt.show()
#plt.savefig('iqspr_like_reorder.png',dpi = 500)
#plt.close()



# re-calculate the property values for the proposed molecules
x_mean, x_std, y_mean, y_std = [], [], [], []
r_std = []
FPs_samples = []
for i, smis in enumerate(iqspr_results_reorder["samples"]):
    tmp_fps = FP_E.transform(smis)
    FPs_samples.append(tmp_fps)
    
    tmp1, tmp2 = custom_mdls1["E"].predict(tmp_fps)
    x_mean.append(tmp1)
    x_std.append(tmp2)
    
    tmp_fps1 = FP_H.transform(smis)
    tmp1, tmp2 = custom_mdls2["HOMO-LUMO gap"].predict(tmp_fps1)
    y_mean.append(tmp1)
    y_std.append(tmp2)
    
    r_std.append([np.sqrt(x_std[-1][i]**2 + y_std[-1][i]**2) for i in range(len(x_std[-1]))])

# flatten the list for max/min calculation
flat_list = [item for sublist in r_std for item in sublist]
print('Range of std. dev.: (%.4f,%.4f)' % (min(flat_list),max(flat_list)))


import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_prd/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

flat_list = np.asarray([item for sublist in r_std for item in sublist])
s_max, s_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["E"],
    np.asarray([item for sublist in x_mean for item in sublist])))
x_max, x_min = max(flat_list), min(flat_list)
flat_list = np.concatenate((data["HOMO-LUMO gap"],
    np.asarray([item for sublist in y_mean for item in sublist])))
y_max, y_min = max(flat_list), min(flat_list)
tmp_beta = iqspr_results_reorder["beta"]

for i in range(len(r_std)):
    dot_size = 45*((np.asarray(r_std[i])-s_min)/(s_max-s_min)) + 5
    
    _ = plt.figure(figsize=(5,5))
    rectangle = plt.Rectangle((0,0),2.5,3,fc='y',alpha=0.1)
    _ = plt.gca().add_patch(rectangle)
    _ = plt.scatter(data["E"], data["HOMO-LUMO gap"],s=3, c='b', alpha=0.2)
    _ = plt.scatter(x_mean[i], y_mean[i],s=dot_size, c='r', alpha=0.5)
    _ = plt.title('Step: %i (beta = %.3f)' % (i,tmp_beta[i]))
    _ = plt.xlim(x_min,x_max)
    _ = plt.ylim(y_min,y_max)
    _ = plt.xlabel('Internal energy')
    _ = plt.ylabel('HOMO-LUMO gap')
    #plt.show()
    _ = plt.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()

import os

# prepare a folder to save all the figures
ini_dir = './iQSPR_tutorial_smiles/'
if not os.path.exists(ini_dir):
    os.makedirs(ini_dir)

tmp_list = [x.sum(axis = 1, skipna = True).values for x in iqspr_results_reorder["loglike"]]    

n_S = 25
for i, smis in enumerate(iqspr_results_reorder['samples']):
    tmp_smis = iqspr_results_reorder['samples'][i][
        np.argsort(tmp_list[i])[::-1]]
    fig, ax = plt.subplots(5, 5)
    _ = fig.set_size_inches(20, 20)
    _ = fig.set_tight_layout(True)
    for j in range(n_S):
        xaxis = j // 5
        yaxis = j % 5
        try:
            img = Draw.MolToImage(Chem.MolFromSmiles(tmp_smis[j]))
            _ = ax[xaxis, yaxis].clear()
            _ = ax[xaxis, yaxis].set_frame_on(False)
            _ = ax[xaxis, yaxis].imshow(img)
        except:
            pass
        _ = ax[xaxis, yaxis].set_axis_off()
    _ = fig.savefig(ini_dir+'Step_%02i.png' % i,dpi = 500)
    _ = plt.close()
    

次はiQSMDとnn_modelで学習モデルを作って、transfer_learningを試してみるかな。