【AI+物质科学】催化反应产率预测#AI夏令营#Datawhale AI夏令营

南墨琛齐

已于 2024-07-31 17:33:47 修改

阅读量956

点赞数 23

文章标签：人工智能

于 2024-07-28 13:35:39 首次发布

本文链接：https://blog.csdn.net/qq_62841331/article/details/140749298

版权

Task1

数据准备

训练数据与测试数据下载连接：上海科学智能研究院 (sais.com.cn)

练集中包含23538条反应数据，测试集中包含2616条反应数据。训练集与测试集的比例接近9：1。每条训练数据包含 rxnid, Reactant1, Reactant2 , Product , Additive, Solvent, Yield 字段。其中 Reactant1 , Reactant2, Product, Additive, Solvent字段中为对应物质的SMILES字符串，Yield字段为目标字段，是经过归一化的浮点数。

即 rxnid, Reactant1, Reactant2, Product, Additive, Solvent为我们进行预测时用到的特征，Yield为我们的输出值

导入库

编写代码时，首先导入我们过程中所需的库

import pickle
import pandas as pd
from tqdm import tqdm
from sklearn.ensemble import RandomForestRegressor
from rdkit.Chem import rdMolDescriptors
from rdkit import RDLogger,Chem
import numpy as np

特征处理

首先对于数据中特征进行了解

rxnid 对数据的id标识，无实际意义
Reactant1 反应物1
Reactant2 反应物2
Product 产物
Additive 添加剂（包括催化剂catalyst等辅助反应物合成但是不对产物贡献原子的部分）
Solvent 溶剂
Yield 产率

由此可以发现rxnid对于我们训练模型的意义不大，过程中可忽略或减弱这一特征的影响

利用rdkit库中函数编写函数将数据特征转换为由0，1组成的位向量，便于我们将数据输入到模型当中进行训练

def mfgen(mol,nBits=2048, radius=2):

    # 返回分子的位向量形式的Morgan fingerprint
    fp = rdMolDescriptors.GetMorganFingerprintAsBitVect(mol,radius=radius,nBits=nBits)
    return np.array(list(map(eval,list(fp.ToBitString()))))

编写加载数据函数，并在过程中使用tqdm库显示进度条，这一操作便于我们了解到加载数据进度，将进度可视化

# 加载数据
def vec_cpd_lst(smi_lst):
    smi_set = list(set(smi_lst))
    smi_vec_map = {}
    for smi in tqdm(smi_set): # tqdm：显示进度条
        mol = Chem.MolFromSmiles(smi)
        smi_vec_map[smi] = mfgen(mol)
    smi_vec_map[''] = np.zeros(2048)
    
    vec_lst = [smi_vec_map[smi] for smi in smi_lst]
    return np.array(vec_lst)

加载数据集

读取训练数据集和测试数据集，并查看训练集和测试集尺寸。

dataset_dir = '../dataset'

train_df = pd.read_csv(f'{dataset_dir}/round1_train_data.csv')
test_df = pd.read_csv(f'{dataset_dir}/round1_test_data.csv')

print(f'Training set size: {len(train_df)}, test set size: {len(test_df)}')

对数据集进行特征处理

读取训练集和测试集中的每个特征为单独的变量，此处我们忽略rxnid这一毫无意义的特征，对每个单独变量使用前面我们编写的特征提取的函数进行处理，将每个特征都转换为位向量，便于后续输入模型，而后将每个特征位向量进行拼接，成为我们输入到模型的特征，而后将Yield转换为数组，是为了使其成为模型输入的标准形式。

train_rct1_smi = train_df['Reactant1'].to_list()
train_rct2_smi = train_df['Reactant2'].to_list()
train_add_smi = train_df['Additive'].to_list()
train_sol_smi = train_df['Solvent'].to_list()

# 将SMILES转化为分子指纹
train_rct1_fp = vec_cpd_lst(train_rct1_smi)
train_rct2_fp = vec_cpd_lst(train_rct2_smi)
train_add_fp = vec_cpd_lst(train_add_smi)
train_sol_fp = vec_cpd_lst(train_sol_smi)
# 在dim=1维度进行拼接。即：将一条数据的Reactant1,Reactant2,Product,Additive,Solvent字段的morgan fingerprint拼接为一个向量。
train_x = np.concatenate([train_rct1_fp,train_rct2_fp,train_add_fp,train_sol_fp],axis=1)
train_y = train_df['Yield'].to_numpy()

# 测试集也进行同样的操作
test_rct1_smi = test_df['Reactant1'].to_list()
test_rct2_smi = test_df['Reactant2'].to_list()
test_add_smi = test_df['Additive'].to_list()
test_sol_smi = test_df['Solvent'].to_list()

test_rct1_fp = vec_cpd_lst(test_rct1_smi)
test_rct2_fp = vec_cpd_lst(test_rct2_smi)
test_add_fp = vec_cpd_lst(test_add_smi)
test_sol_fp = vec_cpd_lst(test_sol_smi)
test_x = np.concatenate([test_rct1_fp,test_rct2_fp,test_add_fp,test_sol_fp],axis=1)

随机森林模型

使用随机森林模型进行训练与预测，随机森林模型可以使用sklearn库中api直接进行调用，并使用fit在训练集上进行训练

model = RandomForestRegressor(n_estimators=10,max_depth=10,min_samples_split=2,min_samples_leaf=1,n_jobs=-1) # 实例化模型，并指定重要参数
model.fit(train_x,train_y)

将模型进行保存

# 保存模型
with open('./random_forest_model.pkl', 'wb') as file:
    pickle.dump(model, file)

加载模型

# 加载模型
with open('random_forest_model.pkl', 'rb') as file:
    loaded_model = pickle.load(file)

使用加载的模型进行预测

# 预测\推理
test_pred = loaded_model.predict(test_x)

若不加载模型直接使用训练好的模型，可在fit训练后直接进行预测

test_pred = model.predict(test_x)

调参调优

调整随机森林模型中的参数，尝试进行优化

model = RandomForestRegressor(n_estimators=20,max_depth=10,min_samples_split=2,min_samples_leaf=1,n_jobs=-1)
model.fit(train_x,train_y)
test_pred = model.predict(test_x)

对参数进行调整，得到更加优化的模型

生成提交所需txt文本文件

根据赛事所需要求，生成相关文本文件

ans_str_lst = ['rxnid,Yield']
for idx,y in enumerate(test_pred):
    ans_str_lst.append(f'test{idx+1},{y:.4f}')
with open('./submit.txt','w') as fw:
    fw.writelines('\n'.join(ans_str_lst))

Task2

经过前一部分的尝试，对于机器学习有了一定的了解，接下来我们选择使用神经网络模型进行训练与预测。对于神经网络模型来说，他的自由性相对于前面随机森林模型这种现成的api来说要更高，我们不仅需要对于模型中每层参数进行调整，还需要对于模型结构本身进行不断优化以寻找最佳模型。

在此次任务中我们选择使用RNN循环神经网络进行尝试，并使用pytorch作为神经网络架构。

导入库

import re
import time
import pandas as pd
from typing import List, Tuple
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader, Subset

数据处理

在使用神经网络模型时，我们选择对于数据进行重新处理，对数据进行分词处理，以便数据能够被模型更好地学习特征。

分词处理在人工智能领域中十分常见，它可以根据我们的需求将数据转换为机器可以识别的类型。在日常使用，例如对于中英文分词，可以使用jieba库，transformer中的官方分词器进行分词较为简便，当我们有自己的需求时，可以自行编写分词器对数据进行处理。

collate_fn函数是在生成dataloader时调用这一函数对数据进行统一处理。

class Smiles_tokenizer():
    def __init__(self, pad_token, regex, vocab_file, max_length):
        self.pad_token = pad_token
        self.regex = regex
        self.vocab_file = vocab_file
        self.max_length = max_length

        with open(self.vocab_file, "r") as f:
            lines = f.readlines()
        lines = [line.strip("\n") for line in lines]
        vocab_dic = {}
        for index, token in enumerate(lines):
            vocab_dic[token] = index
        self.vocab_dic = vocab_dic

    def _regex_match(self, smiles):
        regex_string = r"(" + self.regex + r"|"
        regex_string += r".)"
        prog = re.compile(regex_string)

        tokenised = []
        for smi in smiles:
            tokens = prog.findall(smi)
            if len(tokens) > self.max_length:
                tokens = tokens[:self.max_length]
            tokenised.append(tokens) # 返回一个所有的字符串列表
        return tokenised
    
    def tokenize(self, smiles):
        tokens = self._regex_match(smiles)
        # 添加上表示开始和结束的token：<cls>, <end>
        tokens = [["<CLS>"] + token + ["<SEP>"] for token in tokens]
        tokens = self._pad_seqs(tokens, self.pad_token)
        token_idx = self._pad_token_to_idx(tokens)
        return tokens, token_idx

    def _pad_seqs(self, seqs, pad_token):
        pad_length = max([len(seq) for seq in seqs])
        padded = [seq + ([pad_token] * (pad_length - len(seq))) for seq in seqs]
        return padded

    def _pad_token_to_idx(self, tokens):
        idx_list = []
        for token in tokens:
            tokens_idx = []
            for i in token:
                if i in self.vocab_dic.keys():
                    tokens_idx.append(self.vocab_dic[i])
                else:
                    self.vocab_dic[i] = max(self.vocab_dic.values()) + 1
                    tokens_idx.append(self.vocab_dic[i])
            idx_list.append(tokens_idx)
        
        return idx_list

# 读数据并处理
def read_data(file_path, train=True):
    df = pd.read_csv(file_path)
    reactant1 = df["Reactant1"].tolist()
    reactant2 = df["Reactant2"].tolist()
    product = df["Product"].tolist()
    additive = df["Additive"].tolist()
    solvent = df["Solvent"].tolist()
    if train:
        react_yield = df["Yield"].tolist()
    else:
        react_yield = [0 for i in range(len(reactant1))]
    
    # 将reactant拼到一起，之间用.分开。product也拼到一起，用>分开
    input_data_list = []
    for react1, react2, prod, addi, sol in zip(reactant1, reactant2, product, additive, solvent):
        input_info = ".".join([react1, react2])
        input_info = ">".join([input_info, prod])
        input_data_list.append(input_info)
    output = [(react, y) for react, y in zip(input_data_list, react_yield)]

    # 下面的代码将reactant\additive\solvent拼到一起，之间用.分开。product也拼到一起，用>分开
    '''
    input_data_list = []
    for react1, react2, prod, addi, sol in zip(reactant1, reactant2, product, additive, solvent):
        input_info = ".".join([react1, react2, addi, sol])
        input_info = ">".join([input_info, prod])
        input_data_list.append(input_info)
    output = [(react, y) for react, y in zip(input_data_list, react_yield)]
    '''

    # # 统计seq length，序列的长度是一个重要的参考，可以使用下面的代码统计查看以下序列长度的分布
    # seq_length = [len(i[0]) for i in output]
    # seq_length_400 = [len(i[0]) for i in output if len(i[0])>200]
    # print(len(seq_length_400) / len(seq_length))
    # seq_length.sort(reverse=True)
    # plt.plot(range(len(seq_length)), seq_length)
    # plt.title("templates frequence")
    # plt.show()
    return output

class ReactionDataset(Dataset):
    def __init__(self, data: List[Tuple[List[str], float]]):
        self.data = data
        
    def __len__(self):
        return len(self.data)

    def __getitem__(self, idx):
        return self.data[idx]
    
def collate_fn(batch):
    REGEX = r"\[[^\]]+]|Br?|Cl?|N|O|S|P|F|I|b|c|n|o|s|p|\(|\)|\.|=|#|-|\+|\\\\|\/|:|~|@|\?|>|\*|\$|\%[0-9]{2}|[0-9]"
    tokenizer = Smiles_tokenizer("<PAD>", REGEX, "../vocab_full.txt", max_length=300)
    smi_list = []
    yield_list = []
    for i in batch:
        smi_list.append(i[0])
        yield_list.append(i[1])
    tokenizer_batch = torch.tensor(tokenizer.tokenize(smi_list)[1])
    yield_list = torch.tensor(yield_list)
    return tokenizer_batch, yield_list

定义模型结构

之后我们定义模型结构及参数，这里我们的主体是全连接层，使用三层线性连接层进行连接，并使用sigmoid激活函数，在后续调整中也可尝试relu，softmax等激活函数。

class RNNModel(nn.Module):
    def __init__(self, num_embed, input_size, hidden_size, output_size, num_layers, dropout, device):
        super(RNNModel, self).__init__()
        self.embed = nn.Embedding(num_embed, input_size)
        self.rnn = nn.RNN(input_size, hidden_size, num_layers=num_layers, 
                          batch_first=True, dropout=dropout, bidirectional=True)
        self.fc = nn.Sequential(nn.Linear(2 * num_layers * hidden_size, output_size),
                                nn.Sigmoid(),
                                nn.Linear(output_size,1024 ),
                                nn.Sigmoid(),
                                nn.Linear(1024, 1),
                                nn.Sigmoid())

    def forward(self, x):
        # x : [bs, seq_len]
        x = self.embed(x)
        # x : [bs, seq_len, input_size]
        _, hn = self.rnn(x) # hn : [2*num_layers, bs, h_dim]
        hn = hn.transpose(0,1)
        z = hn.reshape(hn.shape[0], -1) # z shape: [bs, 2*num_layers*h_dim]
        output = self.fc(z).squeeze(-1) # output shape: [bs, 1]
        return output

设定参数进行训练

编写训练函数，设定好我们所需要的参数，例如训练的epoch数，隐藏层神经元数，输出层神经元数，学习率等。

同时使用神经网络模型时，我们不再简单的使用数组将数据输入给模型，我们需要将其转换为张量类型，并打包为dataset,并转换为dataloader，在转换为dataloader时，设定shuffle=True对数据进行充分打乱清洗，batch_size即为一次查看多少组数据等等。

选定优化器，损失函数等即可开始训练

之后编写训练过程，在此次训练过程中，为了防止梯度爆炸问题，我们加入了梯度截断步骤

def train():
    ## super param
    N = 10  #int / int(len(dataset) * 1)  # 或者你可以设置为数据集大小的一定比例，如 int(len(dataset) * 0.1)
    NUM_EMBED = 294 # nn.Embedding()
    INPUT_SIZE = 300 # src length
    HIDDEN_SIZE = 512
    OUTPUT_SIZE = 512
    NUM_LAYERS = 10
    DROPOUT = 0.2
    CLIP = 1 # CLIP value
    N_EPOCHS = 50
    LR = 0.005
    
    start_time = time.time()  # 开始计时
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    # device = 'cpu'
    data = read_data("../dataset/round1_train_data.csv")
    dataset = ReactionDataset(data)
    subset_indices = list(range(N))
    subset_dataset = Subset(dataset, subset_indices)
    train_loader = DataLoader(dataset, batch_size=128, shuffle=True, collate_fn=collate_fn)

    model = RNNModel(NUM_EMBED, INPUT_SIZE, HIDDEN_SIZE, OUTPUT_SIZE, NUM_LAYERS, DROPOUT, device).to(device)
    model.train()
    
    optimizer = optim.Adam(model.parameters(), lr=LR)
    # criterion = nn.MSELoss() # MSE
    criterion = nn.L1Loss() # MAE

    best_loss = 10
    for epoch in range(N_EPOCHS):
        epoch_loss = 0
        for i, (src, y) in enumerate(train_loader):
            src, y = src.to(device), y.to(device)
            optimizer.zero_grad()
            output = model(src)
            loss = criterion(output, y)
            loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), CLIP)
            optimizer.step()
            epoch_loss += loss.item()
            loss_in_a_epoch = epoch_loss / len(train_loader)
        print(f'Epoch: {epoch+1:02} | Train Loss: {loss_in_a_epoch:.3f}')
        if loss_in_a_epoch < best_loss:
            # 在训练循环结束后保存模型
            torch.save(model.state_dict(), '../model/RNN.pth')
    end_time = time.time()  # 结束计时
    # 计算并打印运行时间
    elapsed_time_minute = (end_time - start_time)/60
    print(f"Total running time: {elapsed_time_minute:.2f} minutes")

if __name__ == '__main__':
    train()

生成提交文件

对于神经网络模型，我们根据其特性，需要对生成文件的代码进行少许修改与增添

def predicit_and_make_submit_file(model_file, output_file):
    NUM_EMBED = 294
    INPUT_SIZE = 300
    HIDDEN_SIZE = 512
    OUTPUT_SIZE = 512
    NUM_LAYERS = 10
    DROPOUT = 0.2
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    test_data = read_data("../dataset/round1_test_data.csv", train=False)
    test_dataset = ReactionDataset(test_data)
    test_loader = DataLoader(test_dataset, batch_size=64, shuffle=False, collate_fn=collate_fn) 

    model = RNNModel(NUM_EMBED, INPUT_SIZE, HIDDEN_SIZE, OUTPUT_SIZE, NUM_LAYERS, DROPOUT, device).to(device)
    # 加载最佳模型
    model.load_state_dict(torch.load(model_file))
    model.eval()
    output_list = []
    for i, (src, y) in enumerate(test_loader):
        src, y = src.to(device), y.to(device)
        with torch.no_grad():
            output = model(src)
            output_list += output.detach().tolist()
    ans_str_lst = ['rxnid,Yield']
    for idx,y in enumerate(output_list):
        ans_str_lst.append(f'test{idx+1},{y:.4f}')
    with open(output_file,'w') as fw:
        fw.writelines('\n'.join(ans_str_lst))

    print("done!!!")
    
predicit_and_make_submit_file("../model/RNN.pth",
                              "../output/RNN_submit2.txt")

调优

通过简单测试可以发现，RNN网络模型相比于之前使用的随机森林模型效果较差，这就需要我们对于模型进行不断调参调优。

首先是数据处理上，可以尝试对于分词过程进行进一步优化。

而后可以修改模型结构，例如增加不同的连接层，修改隐藏层神经元数量等

还可以修改学习率，训练的epoch数，调整使用的激活函数。

在过程中还可以加入dropout,early stop，批归一化等技术进行优化，这些均可直接通过pytorch中的api直接调用十分方便。

最后还可以直接选择其他神经网络模型进行尝试，例如卷积神经网络，长短时记忆神经网络模型等。

南墨琛齐

关注

23
点赞
踩
21

收藏

觉得还不错? 一键收藏
0
评论
【AI+物质科学】催化反应产率预测#AI夏令营#Datawhale AI夏令营

练集中包含条反应数据，测试集中包含条反应数据。训练集与测试集的比例接近9：1。每条训练数据包含rxnidReactant1,Reactant2ProductAdditiveSolventYield字段。其中Reactant1Reactant2ProductAdditiveSolvent字段中为对应物质的Yield字段为目标字段，是经过归一化的浮点数。即rxnidReactant1,Reactant2ProductAdditive。
复制链接

扫一扫