【AI for Science】量子化学：分子属性预测-第2次打卡-特征工程baseline上分

第一次打卡跑通了 baseline 【AI for Science】量子化学：分子属性预测-第1次打卡-机器学习baseline

baseline 采用了 CatBoost 模型。
🔗 集成学习之xgboost、lightGBM、CatBoost、NGBoost对比

跑通 baseline 还是比较顺利的。

接下来是继续学习 baseline 代码，并且想办法降低 score。

群里众多大佬提出了非常有用的方法，鱼佬也在b站分享了上分技巧。

【首届世界智能科学大赛】量子化学赛道上分思路与技巧

1. 特征选择

👉 群里的鱼佬 (“鱼遇雨欲语与余”) 提出：

• 增加 connectivity_li 和 edge_li 特征；
• 把 edge_attr 特征加在 edge_list 后面；
• 分析同个分子内相同原子电荷是否相同：formal_charge 特征。

还有朋友提出可以考虑用 coordinates 距离信息。

但是如果把这所有特征都加上，内存会爆掉。后面就是删删补补。

还可以根据后面的特征分数进行筛选。

2. 压缩空间

2.1 内存压缩

内存使用率在慢慢提升

由于云环境内存受限，如果不及时压缩和删除变量，内核可能就会炸了：

👉 群里的 “sheng” 大佬分享了 内存压缩程序

# reduce memory
def reduce_mem_usage(df, verbose=True):
    start_mem = df.memory_usage().sum() / 1024**2
    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']

    for col in df.columns:
        col_type = df[col].dtypes
        if col_type in numerics:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
               
    end_mem = df.memory_usage().sum() / 1024**2
    print('Start Memory usage is: {:.2f} MB'.format(start_mem))
    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))
    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))
    return df

另外，还可以添加下列代码：

gc.collect()

gc.collect()会扫描Python的全局变量和所有创建的对象，并通过引用计数或垃圾回收器判断对象是否仍在使用。如果一个对象不再被使用，gc.collect()就会将其释放掉。

2.2 修改特征提取参数

一般在构建TfidfVectorizer和CountVectorizer特征的时候，内存会增大很多，这时候可以调整TfidfVectorizer和CountVectorizer的参数。

例如：

tfidf = TfidfVectorizer(max_df = 0.95, min_df = 1)

对于出现次数小于1 (min_df) 的数会被过滤掉，可以设置大一点以防过拟合
对于出现频率大于0.95 (max_df) 的数会被过滤掉，可以设置小一点以保留一些信息

另外，Tfidf 和 CVec 可以只保留一个。

3. 模型融合

可以分别用 lightgbm、xgboost 和 catboost 建立预测模型，对三者的预测结果进行融合。融合的方式有很多种，最简单的就是均值融合。

参考：
AI量化模型挑战赛【学习笔记2】
百行代码入手数据挖掘竞赛~

4. 报错解决

问题：
训练 CatBoost 模型时遇到报错：

CatBoostError: Bad value for num_feature[non_default_doc_idx=0,feature_idx=19]="600600 600600 600100 600100 600600 600700 600100 600100 700600 700600 700100 600700 600600 600600 600100 600600 600700 600100 600100 700600 700600 700600 600700 600600 600100 600100 600600 600600 600100 600100 600600 600600 600100 600100 600700 600600 600100 600100 600600 6001600 600100 600100 1600600 1600600 600600 6001600 600100 600100 100600 100600 100600 100600 100700 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600": Cannot convert 'b'600600 600600 600100 600100 600600 600700 600100 600100 700600 700600 700100 600700 600600 600600 600100 600600 600700 600100 600100 700600 700600 700600 600700 600600 600100 600100 600600 600600 600100 600100 600600 600600 600100 600100 600700 600600 600100 600100 600600 6001600 600100 600100 1600600 1600600 600600 6001600 600100 600100 100600 100600 100600 100600 100700 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600 100600'' to float

在这个问题上卡了很久，后来在网上搜到可能是特征顺序问题：CatBoostClassifier功能的静态顺序

如果把新增的 connectivity_li 和 edge_li 特征的引入以及后面的词频统计去掉就没有报错了。

后来发现后面的入模特征没有加上’edge_li’，加上去之后就没有报错了。

cols = [f for f in test_df.columns if f not in ['elements','energy','mol_name','elements','edge_li','istest']]
cat_oof, cat_test = catboost_model(train_df[cols], train_df['energy'], test_df[cols]) 

5. 我的代码

5.1 代码1

加入 edge_li 特征，不使用 CountVectorizer ，微调 TfidfVectorizer 参数

# 导入numpy库，用于进行数值计算
import numpy as np
# 导入pandas库，用于数据处理和分析
import pandas as pd
# 导入polars库，用于处理大规模数据集
import polars as pl
# 导入collections库中的defaultdict和Counter，用于统计
from collections import defaultdict, Counter
# 导入CatBoostRegressor库，用于梯度提升树模型
from catboost import CatBoostRegressor
# 导入StratifiedKFold、KFold和GroupKFold，用于交叉验证
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold
# 使用Joblib中的Parallel和delayed实现并行化处理
from joblib import Parallel, delayed
# 导入sys、os、gc、argparse和warnings库，用于处理命令行参数和警告信息
import sys, os, gc, argparse, warnings, tqdm
# 均绝对误差）是一个用于回归问题评估模型性能的指标之一，它衡量了预测值与实际观测值之间的平均绝对差异。
from sklearn.metrics import mean_absolute_error
import stats
# 忽略警告信息
warnings.filterwarnings('ignore')

# reduce memory
def reduce_mem_usage(df, verbose=True):
    start_mem = df.memory_usage().sum() / 1024**2
    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']

    for col in df.columns:
        col_type = df[col].dtypes
        if col_type in numerics:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
               
    end_mem = df.memory_usage().sum() / 1024**2
    print('Start Memory usage is: {:.2f} MB'.format(start_mem))
    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))
    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))
    return df

path = 'data'
test0 = np.load(f'{path}/QMB_round1_test_230725_0.npy', allow_pickle=True).tolist()
test1 = np.load(f'{path}/QMB_round1_test_230725_1.npy', allow_pickle=True).tolist()
test = test0 + test1
del test0, test1

#一个训练集加载后大概20G，根据自己的算力情况选择性加载，baseline采用A10，30G内存的环境，只加载一个训练集
train0 = np.load(f'{path}/QMB_round1_train_230725_0.npy', allow_pickle=True).tolist()

train = train0
del train0

def get_parallel_feature(data, IS_TRAIN=False):
    # 相连原子组成的列表的最长和最大统计
    max_len = len(max(data['connectivity'], key=len))
    min_len = len(min(data['connectivity'], key=len))
    
    # 提取最大出度和入度，以及边的数量
    # max_out_degree = stats.mode(data['edge_list'][:,0])[1][0]
    # max_in_degree = stats.mode(data['edge_list'][:,1])[1][0]
    edge_list_len = len(data['edge_list'])
    
    # 坐标位置的均值、最大值、最小值
    coordinates = data['coordinates'].mean(axis=0).tolist() + \
                  data['coordinates'].max(axis=0).tolist() + \
                  data['coordinates'].min(axis=0).tolist()
    
    # elements的不同元素数
    elements_nunique = len(set(data['elements']))
    elements = ' '.join([str(i) for i in data['elements']])
    
    # formal_charge最大值和均值
    formal_charge = [data['formal_charge'].max(), data['formal_charge'].mean()]
    
    # edge_attr键类型占比数
    edge_attr_1_ratio = len(np.where(np.array(data['edge_attr'])=='1')[0]) / edge_list_len
    edge_attr_2_ratio = len(np.where(np.array(data['edge_attr'])=='2')[0]) / edge_list_len
    edge_attr_3_ratio = len(np.where(np.array(data['edge_attr'])=='3')[0]) / edge_list_len
    edge_attr_nunique = len(set(data['edge_attr']))
    
    idx2element = dict(zip([i for i in range(data['atom_count'])], data['elements']))
    # 原子边
    edge_li = ' '.join([''.join([str(idx2element[i]*100) for i in li]) for li in data['edge_list']])
                                              
    # 合并到一个list中
    res = [data['mol_name'], data['atom_count'], data['bond_count'], max_len, min_len, edge_list_len] + \
           coordinates + [elements_nunique, elements] + formal_charge + \
          [edge_attr_1_ratio, edge_attr_2_ratio, edge_attr_3_ratio, edge_attr_nunique] + [edge_li]
    
    # 返回结果
    if IS_TRAIN:
        return res + [data['energy']]
    else:
        return res

### 测试数据       
test_samples = Parallel(n_jobs=40)(
    delayed(get_parallel_feature)(data, False)
      for data in tqdm.tqdm(test)
)

test_df = pd.DataFrame(test_samples, columns=['mol_name','atom_count','bond_count','maxlen','maxin','edgelen',\
                      'mean1','mean2','mean3','max1','max2','max3','min1','min2','min3','elements_nunique','elements',\
                      'formal_charge_max','formal_charge_min','edge_attr_1_ratio','edge_attr_2_ratio','edge_attr_3_ratio',\
                      'edge_attr_nunique','edge_li'])

### 训练数据
train_samples = Parallel(n_jobs=40)(
    delayed(get_parallel_feature)(data, True)
      for data in tqdm.tqdm(train)
)

train_df = pd.DataFrame(train_samples, columns=['mol_name','atom_count','bond_count','maxlen','minlen','edgelen',\
                      'mean1','mean2','mean3','max1','max2','max3','min1','min2','min3','elements_nunique','elements',\
                      'formal_charge_max','formal_charge_min','edge_attr_1_ratio','edge_attr_2_ratio','edge_attr_3_ratio',\
                      'edge_attr_nunique','edge_li','energy'])
del train_samples
del test_samples
reduce_mem_usage(train_df)
reduce_mem_usage(test_df)

from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer
def tfidf(data, seqs):
    tfidf = TfidfVectorizer(max_df = 0.92, min_df = 5)
    res = tfidf.fit_transform(data[seqs])
    res = res.toarray()
    for i in range(len(res[0])):
        data['{}_tfidf_{}'.format(seqs,str(i))] = res[:,i]
    gc.collect()
    return data

### 合并训练数据和测试数据
test_df['istest'] = 1
train_df['istest'] = 0
df = pd.concat([test_df, train_df], axis=0, ignore_index=True)

### 进行Tfidf
# elements
df = tfidf(df,'elements')

# edge_li
df = tfidf(df,'edge_li')

### 切分训练数据和测试数据
test_df = df[df.istest==1].reset_index(drop=True)
train_df = df[df.istest==0].reset_index(drop=True)
del df
gc.collect()

def catboost_model(train_x, train_y, test_x, seed = 2023):
    folds = 5
    kf = KFold(n_splits=folds, shuffle=True, random_state=seed)
    oof = np.zeros(train_x.shape[0])
    test_predict = np.zeros(test_x.shape[0])
    cv_scores = []
    
    for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
        print('************************************ {} ************************************'.format(str(i+1)))
        trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
        
        params = {'learning_rate': 0.1, 
                  'depth': 4, 
                  'bootstrap_type':'Bernoulli',
                  'random_seed':2023,
                  'od_type': 'Iter', 
                  'od_wait': 200, 
                  'random_seed': 11, 
                  'allow_writing_files': False,       
                  'task_type':"GPU",  # 任务类型，表示模型运行在GPU还是CPU上。设置为"GPU"表示模型运行在GPU上，如果计算机没有GPU，可以设置为"CPU"。
                  'devices':'0:1'}
        
        #iterations是迭代次数，可以根据自己的算力配置与精力调整
        model = CatBoostRegressor(iterations=10000, **params)
        model.fit(trn_x, trn_y, eval_set=(val_x, val_y),
                  metric_period=400,
                  use_best_model=True, 
                  cat_features=[],
                  verbose=1)

        val_pred  = model.predict(val_x)
        test_pred = model.predict(test_x)
        
        oof[valid_index] = val_pred
        test_predict += test_pred / kf.n_splits
        
        score = mean_absolute_error(val_y, val_pred)
        cv_scores.append(score)
        print(cv_scores)
        
        # 获取特征重要性打分，便于评估特征
        if i == 0:
            fea_ = model.feature_importances_
            fea_name = model.feature_names_
            fea_score = pd.DataFrame({'fea_name':fea_name, 'score':fea_})
            fea_score = fea_score.sort_values('score', ascending=False)
            fea_score.to_csv('feature_importances8.csv', index=False)
        
    return oof, test_predict
cols = [f for f in test_df.columns if f not in ['elements','energy','mol_name','elements','edge_li','istest']]
cat_oof, cat_test = catboost_model(train_df[cols], train_df['energy'], test_df[cols]) 

# 输出赛题提交格式的结果
test_df['energy'] = cat_test
test_df['force'] = test_df['atom_count'].apply(lambda x: ','.join(['0.0' for _ in range(x*3)]))

test_df[['energy','force']].to_csv("submission8.csv", index=True)

跑完内存使用率90+%，还是挺危险的。
分数：17715.0737
相比之前降低了将近4000 👍

5.2 代码2

删去 max_len 、min_len、edge_attr_3_ratio 、formal_charge_max 、edge_attr_nunique 特征，加入 connectivity_li 特征，最终还是爆内存了 😞

5.3 代码3

融合 lightgbm 、catboost 和 xgboost 模型，加入 edge_list 特征。
先做 lightgbm + catboost 的。

其余代码不变，加入lighgbm：

# LGBMRegressor: 调用Scikit-learn接口的回归
def lgb_model(train_x, train_y, test_x, seed = 2023):
    folds = 5
    kf = KFold(n_splits=folds, shuffle=True, random_state=seed)
    oof = np.zeros(train_x.shape[0])
    test_predict = np.zeros(test_x.shape[0])
    cv_scores = []
    
    for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
        print('************************************ {} ************************************'.format(str(i+1)))
        trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
        
        model = lgb.LGBMRegressor(objective='regression', 
                                  num_leaves=500, 
                                  learning_rate=0.1, 
                                  n_estimators=500)
        
        model.fit(trn_x, trn_y, eval_set=(val_x, val_y))

        val_pred  = model.predict(val_x)
        test_pred = model.predict(test_x)
        
        oof[valid_index] = val_pred
        test_predict += test_pred / kf.n_splits
        
        score = mean_absolute_error(val_y, val_pred)
        cv_scores.append(score)
        print(cv_scores)
        
    return oof, test_predict
cols = [f for f in test_df.columns if f not in ['elements','energy','mol_name','elements','edge_li','istest']]
lgb_oof, lgb_test = lgb_model(train_df[cols], train_df['energy'], test_df[cols]) 

然后取均值：

# 输出赛题提交格式的结果
# test_df['energy'] = (cat_test + lgb_test + xgb_test)/3
test_df['energy'] = (cat_test + lgb_test)/2
test_df['force'] = test_df['atom_count'].apply(lambda x: ','.join(['0.0' for _ in range(x*3)]))

test_df[['energy','force']].to_csv("submission11.csv", index=True)

分数：17895.7788

再次调整：调整 CatBoostRegressor 的 ‘depth’: 6，CatBoostRegressor 和 LGBMRegressor 的 folds = 6，删去 max_len 和 edge_attr_3_ratio 特征。
分数：15898.3176
又有了一点点提升 👍