Datawhale AI夏令营 ---催化反应产率预测

比赛链接:上海科学智能研究院 (sais.com.cn)

---

先对比赛的数据进行分析

官方发布的数据是化学分子的SMILES表达式。

SMILES

SMILES(Simplified Molecular Input Line Entry System)是一种将化学分子表达式转换为ASCII码的形式的方法，是化学信息学领域很重要的工具。

通过分析数据，我们很容易就可以知道：

Reactant1、Reactant2表示反应物

Product表示产物

Additive表示添加剂(不对反应贡献原子)

Solvent表示溶剂

Yield表示产率

以上数据除rxnid列之外其余列的数据都用SMILES表达。因此可以通过SMILES提取分子指纹(向量)来作为特征以便后续的训练

分子指纹

分子指纹(Molecular Fingerprint)是化学信息学中一种用于表示分子结构的简洁而有效的方式。它将分子结构转换成固定长度的位向量(bit vector)，每一位（bit）表示分子中某种特定的结构特征是否存在。

比如乙醇分子用SMILES表示为(CCO)，使用Morgan指纹(分子指纹的一种，广泛用于回归)转化后可以表示为以下16位的位向量:

[0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1]

接下来就上代码，我们先导入必要的库

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
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
RDLogger.DisableLog('rdApp.*')

然后写两个函数用于特征提取

def mfgen(mol, nBits=2048, radius=2):
    '''
    Parameters
    ----------
    mol : mol
        RDKit mol object.
    nBits : int
        Number of bits for the fingerprint.
    radius : int
        Radius of the Morgan fingerprint.
    Returns
    -------
    mf_desc_map : ndarray
        ndarray of molecular fingerprint descriptors.
    '''
    # 返回分子的位向量形式的Morgan fingerprint
    fp = rdMolDescriptors.GetMorganFingerprintAsBitVect(mol, radius=radius,nBits=nBits)
    return np.array(list(map(eval, list(fp.ToBitString()))))


# 加载数据
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)

接下来正式开始读入数据和预处理数据(即将数据转化为分子指纹)

# 先读入数据
train_df = pd.read_csv('round1_train_data.csv')
test_df = pd.read_csv('round1_test_data.csv')

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

# 从csv中读取数据
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)

接下来便可以开始进行模型训练和预测了

目前使用的是随机森林预测法

随机森林

随机森林(Random Forest)，一种流行且强大的机器学习方法，适用于分类和回归任务。它通过结合多个决策树的预测结果来提高模型的准确性和鲁棒性。

随机森林的优点

        高准确性：通过集成多个决策树，随机森林通常比单个决策树具有更高的预测准确性。

        抗过拟合：由于引入了数据和特征的随机性，随机森林在处理高维数据和复杂模型时，具有较好的抗过拟合能力。

        处理缺失值：随机森林能够处理数据中的缺失值。

        特征重要性：随机森林可以计算特征的重要性，有助于特征选择和理解模型。

随机森林的缺点

        计算复杂度高：由于需要训练大量的决策树，随机森林的训练和预测时间较长。

        内存消耗大：随机森林需要存储多个决策树，可能会占用较多内存。

上代码

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

解释以下各个参数：

n_estimators:决策树的个数(这里设置为10颗)，越多越好，但是越多计算开销越大

max_depth:默认为None,设置树的最大深度

min_samples_split:根据属性划分节点时最少样本数(这里设置为5)

min_samples_leaf:叶子节点的最少样本数(这里设置为1)

n_jobs:并行job个数-1表示使用所有CPU进行并行计算，一般设置为-1或4(这里设置为-1)

模拟拟合完毕后，就可以正常进行预测啦

# 保存模型
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)

# 生成赛题的submit.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))

但是提交数据之后发现分数其实并不理想，只有0.18分

原因之一可能是模型的参数不太理想，我们可以尝试使用超参数调优的方法。

超参数调优

超参数指的是模型在训练之前需要设置的参数，它们对模型性能有显著影响。为了找到最佳的超参数组合，我们可以使用网格搜索(Grid Search)和随机搜索(Random Search)两种方法。

网格搜索

网格搜索(Grid Search)是一种系统的超参数调优方法，通过穷举所有可能的超参数组合来找到最佳的超参数集。

优点

使用简单方便，易于理解

遍历了所有的超参数组合，确保找到全局最优解

缺点

计算成本高，时间消耗大

在大数据集下，计算时间不可接受

以下是针对这次赛题的Grid Search代码:

# 导入Grid Search需要的库
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV
# 数据标准化和模型训练的流水线
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('model', RandomForestRegressor(n_jobs=-1))
])

# 定义参数网格
param_grid = {
    'model__n_estimators': [1000, 3000, 4500, 5000],
    'model__max_depth': [None, 10, 20, 30],
    'model__min_samples_split': [2, 5, 10],
    'model__min_samples_leaf': [1, 2, 4],
    'model__max_features': ['auto', 'sqrt', 'log2']
}

# 网格搜索
grid_search = GridSearchCV(estimator=pipeline, param_grid=param_grid, cv=5, n_jobs=4, scoring='r2')
grid_search.fit(train_x, train_y)

# 打印最佳参数和得分
print("Best parameters found: ", grid_search.best_params_)
print("Best cross-validation score: ", grid_search.best_score_)

在经过长时间(一般是几个小时，和自己给定的搜索范围大大相关)的搜索后，可以得到模型的最佳参数，即grid_search.best_params_，我们可以直接提取里面的数据然后进行训练。

# 使用最佳参数训练模型, 然后拟合
best_params = grid_search.best_params_
model = RandomForestRegressor(
    n_estimators=best_params['model__n_estimators'],
    max_depth=best_params['model__max_depth'],
    min_samples_split=best_params['model__min_samples_split'],
    min_samples_leaf=best_params['model__min_samples_leaf'],
    max_features=best_params['model__max_features'],
    n_jobs=-1
)
model.fit(train_x, train_y)

随机搜索

随机搜索(Random Search)是一种通过随机采样超参数组合来进行超参数调优的方法。

优点

相比网格搜索，随机搜索在高维参数空间中的效率更高。

可以通过设定搜索次数，在更大的参数空间中进行搜索，减少计算时间。

缺点

由于是随机采样，可能会错过某些潜在的最佳超参数组合。

以下是针对本次赛题的 Random Search代码:

# 导入Random Search需要的库
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import RandomizedSearchCV

# 数据标准化和模型训练的流水线
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('model', RandomForestRegressor(n_jobs=-1))
])

# 定义参数分布
param_distributions = {
    'model__n_estimators': [100, 500, 1000, 2000, 3000],
    'model__max_depth': [None, 10, 20, 30, 40],
    'model__min_samples_split': [2, 5, 10, 15],
    'model__min_samples_leaf': [1, 2, 4, 6],
    'model__max_features': ['auto', 'sqrt', 'log2']
}

# 随机搜索
random_search = RandomizedSearchCV(estimator=pipeline, param_distributions=param_distributions, 
                                   n_iter=50, cv=5, n_jobs=4, scoring='r2', random_state=42)

# 训练模型
random_search.fit(train_x, train_y)

# 打印最佳参数和得分
print("Best parameters found: ", random_search.best_params_)
print("Best cross-validation score: ", random_search.best_score_)

同样地，在经过一定时间的搜索后，可以得到模型的最佳参数，即random_search.best_params_，我们可以直接提取里面的数据然后进行训练。

# 使用最佳参数训练模型, 然后拟合
best_params = random_search.best_params_
model = RandomForestRegressor(
    n_estimators=best_params['model__n_estimators'],
    max_depth=best_params['model__max_depth'],
    min_samples_split=best_params['model__min_samples_split'],
    min_samples_leaf=best_params['model__min_samples_leaf'],
    max_features=best_params['model__max_features'],
    n_jobs=-1
)
model.fit(train_x, train_y)

经过超参数调优的方法之后，我们便可以选取出给定范围的最佳参数啦，通过提交的submit.txt文件可知，分数虽然不算太高，但是提升了近2倍

其它的方法有待探索。

补充

使用以上工具搜索的时候可能会遇到TerminatedWorkerError，具体描述可能如下:

“A worker process managed by the executor was unexpectedly terminated. This could be caused by a segmentation fault while calling the function or by an excessive memory usage causing the Operating System to kill the worker.”

解决方法是将GridSearchCV或者RandomizedSearchCV函数里面的n_jobs参数由-1改为4或者其它合适的数字即可(如果一开始设置为-1的话)。