Datawhale AI夏令营（AI+物质）Task2：RNN建模SMILES进行反应产率预测

Task2：RNN建模SMILES进行反应产率预测部分

摘要的学习手册知识点：

一、基础知识点
（简单说明,详细见手册，链接见本节末） 

1、Ai4Chemistry发展历史
在前期阶段，人们主要尝试使用不同的方法，尽可能地将化学知识和信息以计算机的形式进行存储，并以此为基础开始构建数据库；中期阶段，使用一些手动的特征工程对已有数据进行编码、特征提取等操作。再辅以传统的机器学习的方法，做一些预测；后期阶段，随着深度学习网络的广泛使用，这也导致后来又非常多的新型的分子指纹出现。

2、SMILES —— 最流行的将分子表示为序列类型数据的方法
SMILES，提出者Weininger et al[1]，全称是Simplified Molecular Input Line Entry System，是一种将化学分子用ASCII字符表示的方法，在化学信息学领域有着举足轻重的作用。

SMILES将化学分子中涉及的原子、键、电荷等信息，用对应的ASCII字符表示；环、侧链等化学结构信息，用特定的书写规范表达。

在SMILES中，原子由他们的化学符号表示，=表示双键、#表示三键、[]里面的内容表示侧基或者特殊原子（例如[Cu+2]表示带电+2电荷的Cu离子）。 

3、 分子指纹 —— 分子向量化
分子的指纹就像人的指纹一样，用于表示特定的分子。分子指纹是一个具有固定长度的位向量（即由0，1组成），其中，每个为1的值表示这个分子具有某些特定的化学结构。例如，对于一个只有长度为2的分子指纹，我们可以设定第一个值表示分子是否有甲基，第二个位置的值表示分子是都有苯环，那么[0,1]的分子指纹表示的就是一个有苯环而没有甲基的分子。

4、RDkit

RDkit是化学信息学中主要的工具，是开源的。网址：http://www.rdkit.org，支持WIN\MAC\Linux，可以被python、Java、C调用。

备注：在手册中也充分说明介绍了，在次就不进行阐述介绍了。手册链接：Datawhale (linklearner.com)

二、机器学习
机器学习按照目标可以分为分类任务（classification）和回归（regression）任务两大类。所谓分类任务，就是模型预测的结果是离散的值，例如类别；那么，回归任务中，模型预测的结果就是连续的值，例如房价等等。

就比说在task1的学习中，就使用了随机森林模型的机器学习方法。将多个决策树结合在一起，训练每个决策树的数据集都是随机有放回地从原数据中选出。预测的时候，输入会通过每个决策树进行预测，然后考虑每个树地输出结果，得到最终的预测值。

三、深度学习
深度学习可以归为机器学习的一个子集，主要通过神经网络学习数据的特征和分布。深度学习的一个重要进化是不再需要繁琐的特征工程，让神经网络自己从里面学习特征。

SMILES是一种以ASCII组成的序列，可以被理解为一种“化学语言”。既然是一种语言，那么很自然地想到了可以使用NLP中的方法对SMILES进行建模。

使用RNN对SMILES建模是早期的一个主要方法。RNN（Recurrent Neural Network）是处理序列数据的一把好手。RNN的网络每层除了会有自己的输出以外，还会输出一个隐向量到下一层。

其中，每一层相当于做了一次线性变换：

h_n = \sigma(W_{hh}h_{n-1} + W_{hx}x_n + b_n)

每层的输出： y_n = Softmax(Vh_n + c)

通过隐向量的不断传递，序列后面的部分就通过“阅读”隐向量，获取前面序列的信息，从而提升学习能力。

但是RNN也有缺点：如果序列太长，那么两个相距比较远的字符之间的联系需要通过多个隐藏向量。这就像人和人之间传话一样，传递的人多了，很容易导致信息的损失或者扭曲。因此，它对长序列的记忆能力较弱。

同时，RNN需要一层一层地传递，所以并行能力差，比较容易出现梯度消失或梯度爆炸问题。

下面是代码每一行的分析。

 一、导入库和包

import re  # 导入正则表达式库，用于字符串处理
import time  # 导入时间库，用于计时
import pandas as pd  # 导入Pandas库，用于数据处理
from typing import List, Tuple  # 导入类型提示工具
import torch  # 导入PyTorch库
import torch.nn as nn  # 导入PyTorch的神经网络模块
import torch.optim as optim  # 导入PyTorch的优化器模块
from torch.utils.data import Dataset, DataLoader, Subset  # 导入数据集和数据加载相关工具

 二、定义RNN模型

class RNNModel(nn.Module):  # 定义一个RNN模型类，继承自PyTorch的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)  # 定义双向RNN层
        self.fc = nn.Sequential(  # 定义全连接层，用于最终的预测输出
            nn.Linear(2 * num_layers * hidden_size, output_size),  # 全连接层，将RNN输出转换为特定维度
            nn.Sigmoid(),  # Sigmoid激活函数
            nn.Linear(output_size, 1),  # 再次全连接层，将输出转换为单一值
            nn.Sigmoid())  # 再次应用Sigmoid激活函数

    def forward(self, x):  # 定义前向传播函数
        x = self.embed(x)  # 将输入数据通过嵌入层进行转换
        _, hn = self.rnn(x)  # 通过RNN层，返回隐藏状态
        hn = hn.transpose(0,1)  # 转置隐藏状态的维度
        z = hn.reshape(hn.shape[0], -1)  # 将隐藏状态展平，准备输入全连接层
        output = self.fc(z).squeeze(-1)  # 通过全连接层计算输出，并压缩最后一个维度
        return output  # 返回模型的输出

 

 三、定义SMILES Tokenizer

class Smiles_tokenizer():  # 定义一个SMILES Tokenizer类
    def __init__(self, pad_token, regex, vocab_file, max_length):
        self.pad_token = pad_token  # 初始化填充符号
        self.regex = regex  # 初始化正则表达式，用于匹配SMILES中的元素
        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:  # 遍历所有SMILES字符串
            tokens = prog.findall(smi)  # 使用正则表达式匹配SMILES中的元素
            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)  # 调用内部的正则表达式匹配方法
        tokens = [["<CLS>"] + token + ["<SEP>"] for token in tokens]  # 添加特殊标记<CLS>和<SEP>
        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)  # 使用Pandas读取CSV文件
    
    reactant1 = df["Reactant1"].tolist()  # 提取Reactant1列并转换为列表
    reactant2 = df["Reactant2"].tolist()  # 提取Reactant2列并转换为列表
    product = df["Product"].tolist()  # 提取Product列并转换为列表
    additive = df["Additive"].tolist()  # 提取Additive列并转换为列表
    solvent = df["Solvent"].tolist()  # 提取Solvent列并转换为列表
    
    if train:  # 如果是训练数据
        react_yield = df["Yield"].tolist()  # 提取Yield列并转换为列表
    else:  # 如果是测试数据
        react_yield = [0 for i in range(len(reactant1))]  # 用0填充产率列表
    
    input_data_list = []  # 初始化输入数据列表
    for react1, react2, prod, addi, sol in zip(reactant1, reactant2, product, additive, solvent):  # 遍历所有反应信息
        input_info = ".".join([react1, react2])  # 将Reactant1和Reactant2用"."连接
        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)]  # 将输入数据和产率组合成元组
    return output  # 返回处理后的数据

 五、定义自定义数据集

class ReactionDataset(Dataset):  # 定义自定义数据集类，继承自PyTorch的Dataset
    def __init__(self, data: List[Tuple[List[str], float]]):  # 初始化方法，接收数据
        self.data = data  # 将数据赋值给实例变量self.data
        
    def __len__(self):  # 返回数据集的长度
        return len(self.data)  # 返回self.data的长度
 
    def __getitem__(self, idx):  # 根据索引获取数据集中的元素
        return self.data[idx]  # 返回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)  # 创建SMILES Tokenizer对象
    
    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])  # 使用Tokenizer将化学分子式转为张量
    yield_list = torch.tensor(yield_list)  # 将产量转换为张量
    
    return tokenizer_batch, yield_list  # 返回标记化后的张量和产量张量

七、训练模型

def train():  # 定义训练函数
    N = 10  # 设置数据集大小的一部分，或数据集大小的一定比例
    NUM_EMBED = 294  # 嵌入层的维度
    INPUT_SIZE = 300  # 输入序列长度
    HIDDEN_SIZE = 512  # RNN隐藏层大小
    OUTPUT_SIZE = 512  # 输出层大小
    NUM_LAYERS = 10  # RNN层数
    DROPOUT = 0.2  # Dropout概率
    CLIP = 1  # 梯度裁剪阈值
    N_EPOCHS = 100  # 训练轮数
    LR = 0.0001  # 学习率
 
    start_time = time.time()  # 记录开始时间
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')  # 检查是否有GPU可用，设置设备
    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)  # 初始化RNN模型并移动到设备上
    model.train()  # 将模型设置为训练模式
 
    optimizer = optim.Adam(model.parameters(), lr=LR)  # 定义Adam优化器
    criterion = nn.L1Loss()  # 使用L1损失函数
 
    best_loss = 10  # 初始化最佳损失值
    for epoch in range(N_EPOCHS):  # 开始训练，迭代每个epoch
        epoch_loss = 0  # 初始化当前epoch的损失值
        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)  # 计算当前epoch的平均损失值
        print(f'Epoch: {epoch+1:02} | Train Loss: {loss_in_a_epoch:.3f}')  # 打印当前epoch的损失值
        if loss_in_a_epoch < best_loss:  # 如果当前epoch的损失值小于最佳损失值
            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()  # 调用训练函数

八、预测并生成提交文件

import torch  # 导入PyTorch
from torch.utils.data import DataLoader  # 导入数据加载器
from dataset import ReactionDataset, collate_fn  # 导入自定义数据集和数据整理函数
from model import RNNModel  # 导入RNN模型
from data_utils import read_data  # 导入数据读取函数

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  # RNN层数
    DROPOUT = 0.2  # Dropout概率
    
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')  # 检查是否有GPU可用，设置设备
    
    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)  # 初始化RNN模型并加载到设备上
    model.load_state_dict(torch.load(model_file))  # 加载已训练好的模型参数
    
    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_submit.txt")  # 调用预测函数并生成提交文件

我实际是在自己的本地服务器上运行的，初次运行就花了两个小时左右，但是成绩和随机森林的结果完全不能比。。。

测试结果如下

（不知道是不是这个模型并不适合这个项目