智源小分子预测赛进行中：基于CatBoost建模的baseline分享

import pandas as pd

import os

import tqdm

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sns



# 注意改成自己的路径

train_data = pd.read_csv('molecule_open_data/candidate_train.csv')

test_data = pd.read_csv('molecule_open_data/candidate_val.csv')

sample_submission = pd.read_csv('molecule_open_data/sample_submission.csv')

train_answer = pd.read_csv('molecule_open_data/train_answer.csv')



def show_info(df):

    print(df.shape)

    print(df.columns)

    print()



# 查看数据的 shape 和 columns

show_info(train_data)

show_info(test_data)

show_info(sample_submission)

show_info(train_answer)

# 将 p1 ~ p6 拼接到原来的特征

train_data = train_data.merge(train_answer, on='id', how='left')

show_info(train_data)

train_data.head()

# 查看特征的数据类型。有 21 个特征是浮点数，3156 个特征是整数

print(train_data.info())

print()

print(test_data.info())

# 查看是否有缺失值，结果显示没有缺失值。

train_data.isna().sum().sum()

# 查看 p1 ~ p6 的统计指标，如 min, max, mean 等。

train_answer.describe()

# 查看 p1 ~ p6 的分布，可以看到，p5、p6 十分接近

# 所有 property 离群点都不多，p1,2,4,5,6 接近正态分布。

# p1 在右边有个小长尾，也许可以试试对 p1 做一下 log 变换

from scipy.stats import norm

f, ax = plt.subplots(3, 2, figsize=(10, 10))


sns.distplot(train_answer['p1'], kde=True, fit=norm, ax=ax[0, 0])
sns.distplot(train_answer['p2'], kde=True, fit=norm, ax=ax[0, 1])
sns.distplot(train_answer['p3'], kde=True, fit=norm, ax=ax[1, 0])
sns.distplot(train_answer['p4'], kde=True, fit=norm, ax=ax[1, 1])
sns.distplot(train_answer['p5'], kde=True, fit=norm, ax=ax[2, 0])
sns.distplot(train_answer['p6'], kde=True, fit=norm, ax=ax[2, 1])

# 查看 p1 ~ p6 之间的相关性。可以看到，p5 p6 相关性为 1，p3 p4 相关性 0.65。除此之外， property 之间的相关性都比较低。

# 也许可以用某个 property 作为另一个 property 的特征。（这个也许可以尝试）

sns.heatmap(train_answer.corr(), annot=True, linewidths=0.2)

fig=plt.gcf()

fig.set_size_inches(7, 6)

select_feats = ['34', '3157', '3158', '3175']


f, ax = plt.subplots(len(select_feats), 2, figsize=(10, 5 * len(select_feats)))
for i in range(len(select_feats)):

    sns.distplot(train_data[select_feats[i]], kde=True, ax=ax[i, 0])

    sns.distplot(test_data[select_feats[i]], kde=True, ax=ax[i, 1])

    ax[i, 0].set_title('Train data')

    ax[i, 1].set_title('Test data')



# 注：可以看到，在这四个特征上，训练集和测试集的分布还是十分接近的。

# 3157，3158，3175 为数值型特征（连续的）。34 从图上看似乎只有 0,1,2 三个取值。

# 3157 和 3158 不仅取值范围相近，而且密度曲线也很接近。

# （3157,1358） 和 3175 虽然取值范围不同，但密度曲线比较接近。

# 使用 value_counts() 方法， 可以看到，特征 34 的取值实际上有 0,1,2,3,4。只不过特征值为 3、4 的样本数占比太少。

print(train_data['34'].value_counts())

print()

print(test_data['34'].value_counts())

# 特征 34 的分布

f, ax=plt.subplots(1, 2, figsize=(10,5))



explode = [0.03] * train_data['34'].nunique()

train_data['34'].value_counts().plot(kind='bar', ax=ax[0])

train_data['34'].value_counts().plot.pie(explode=explode, autopct='%1.1f%%', ax=ax[1],shadow=True)



# 注：在后面的模型训练中，我将所有特征当成连续数值型特征来处理。

# 但其实也可以尝试把某些特征（比如这里的 34）当成类别型特征处理，做一些特征工程，看效果有没有提升

# 直接选取特征 0 到 20。（可以自己修改 show_num 的值，查看更多特征的分布）

show_num = 20

select_feats = [str(i) for i in range(show_num)]

f, ax = plt.subplots(len(select_feats), 2, figsize=(10, 5 * len(select_feats)))

for i in range(len(select_feats)):

    sns.distplot(train_data[select_feats[i]], kde=True, ax=ax[i, 0])

    sns.distplot(test_data[select_feats[i]], kde=True, ax=ax[i, 1])

    ax[i, 0].set_title('Train data')

# 各个特征在训练集和测试集上的分布依旧一致。
# 一个奇怪的现象是，某些特征（0,1,2,3,4,5,19）似乎只有一个取值。我们在下个 cell 验证一下

single_value_cols = []

columns = [str(i) for i in range(3177)]

for col in columns:

    if train_data[col].nunique() == 1:

        print(col, train_data[col].unique(), test_data[col].unique())

        single_value_cols.append(col)



print('只有一个取值的 col 数：', len(single_value_cols))

# 可以看到，有 51 个特征，只有一个取值 0。可以删掉这些特征。（注：本文并没有删除这些特征）

# 拼接训练集和测试集，做特征的时候就不用分开做了

data = pd.concat([train_data, test_data], axis=0, sort=False).reset_index(drop=True)



from sklearn.preprocessing import PolynomialFeatures



# 不做特征工程的时候跑过一次模型，筛选出来的每个 property 重要性排前 4 的特征

high_imps_list = [

    ['25', '3175', '34', '3158'],

    ['25', '16', '3175', '3157'],

    ['34', '10', '3157', '15'],

    ['3161', '3155', '3173', '3157'],

    ['3155', '3173', '3175', '3151'],

    ['3155', '3173', '3175', '3151']

]



# 每个 property 使用不同的多项式特征

pf_df_list = []

for i in range(len(high_imps_list)):

    high_imps = high_imps_list[i]

    pf = PolynomialFeatures(degree=5)

    pf_array = pf.fit_transform(data[high_imps])

    pf_df = pd.DataFrame(pf_array)

    pf_df.columns = [f'p{i+1}_pf{col}' for col in pf_df.columns]

    

    # 实际上，pf_0 取值全为 1，pf_1,2,3,4 其实就是 high_imps 的四个特征

    drop_columns = [f'p{i+1}_pf{j}' for j in range(len(high_imps) + 1)]

    pf_df.drop(columns=drop_columns, inplace=True)

    pf_df_list.append(pf_df)

    

# 可以看到，生成的特征比较多，但估计只有部分有用。可以先全部加上，后面慢慢删（根据特征重要性等）。

pf_df_list[0].head()

# 将生成的特征拼接到原特征中(注：后面分类时不同 property 使用不同的多项式特征)

data = pd.concat([data] + pf_df_list, axis=1).reset_index(drop=True)

show_info(data)

import gc

import pickle

import datetime



import numpy as np

from catboost import CatBoostRegressor, Pool

from sklearn.model_selection import KFold

# 6 个 property 共用的初始 3177 特征

use_cols = [col for col in data.columns if col != 'id' and 'p' not in col]

print('Number of common used features:', len(use_cols))

%%time

# 交叉验证折数和待预测的 property 数

n_fold = 5  

n_property = 6  

train = data[data['p1'].notna()].reset_index(drop=True)

test = data[data['p1'].isna()].reset_index(drop=True)

print('Shape of train dataset:', train.shape)

print('Shape of test dataset:', test.shape)


# 测试集和训练集的预测结果

pred_test = np.zeros((len(test), n_property))

pred_train = np.zeros((len(train), n_property))


# 保存特征重要性的 DataFrame

imps = [pd.DataFrame() for _ in range(n_property)]

smape_score = np.zeros((n_fold, n_property))

all_index = range(len(train))

kfold = KFold(n_splits=n_fold, shuffle=True, random_state=1)

for fold_idx, (tr_idx, va_idx) in enumerate(kfold.split(all_index)):

    time_stamp = datetime.datetime.now()

    print('*' * 120)

    print(f"Fold [{fold_idx}]: " + time_stamp.strftime('%Y.%m.%d-%H:%M:%S'))

    temp_train = train.iloc[tr_idx]

    temp_valid = train.iloc[va_idx]

    

    for p_idx in range(n_property):

        print(f'Training for p{p_idx + 1} ......')

        

        # 当前 property 使用的特征

        temp_cols = use_cols + [col for col in data.columns if f'p{p_idx+1}_' in col]

        print('Length of temp_cols:', len(temp_cols))

        

        label = 'p' + str(p_idx + 1)

        x_train, x_valid = temp_train[temp_cols], temp_valid[temp_cols]

        y_train, y_valid = temp_train[label], temp_valid[label]

        cate_features = []

        train_pool = Pool(x_train, y_train, cat_features=cate_features)

        eval_pool = Pool(x_valid, y_valid, cat_features=cate_features)

    

        cbt_model = CatBoostRegressor(iterations=10000, # 注：baseline 提到的分数是用 iterations=60000 得到的，但运行时间有点久

                           learning_rate=0.1, # 注：事实上好几个 property 在 lr=0.1 时收敛巨慢。后面可以考虑调大

                           eval_metric='SMAPE',

                           use_best_model=True,

                           random_seed=42,

                           logging_level='Verbose',

                           task_type='GPU',

                           devices='0',

                           gpu_ram_part=0.5,

                           early_stopping_rounds=400,

                           )

        cbt_model.fit(train_pool, eval_set=eval_pool, verbose=1000)

        smape_score[fold_idx][p_idx] = cbt_model.best_score_['validation']['SMAPE']

        

        # 模型保存

        if not os.path.exists('models'):

            os.mkdir('models')

        with open(f'./models/cbt_f{len(use_cols)}_fold_{fold_idx}_p{p_idx + 1}.pkl', 'wb') as f:

            pickle.dump(cbt_model, f)

    

        if fold_idx == 0:

            imps[p_idx]['feat'] = temp_cols

        imps[p_idx][f'score{fold_idx}'] = cbt_model.feature_importances_

    

        # 注：pred_train 保存的是 oof 预测结果，可用于后面的 stacking 等操作。（baseline 里其实没有用到 pred_train）

        pred_train[va_idx, p_idx] = cbt_model.predict(x_valid)

        pred_test[:, p_idx] += cbt_model.predict(test[temp_cols]) / 5

        

        print()

        

    print('smape_score:', smape_score[fold_idx])

    print(f'fold {fold_idx} round {cbt_model.best_iteration_} : smape_mean: {np.mean(smape_score[fold_idx]):.6f}')



    del cbt_model, train_pool, eval_pool

    del x_train, y_train, x_valid, y_valid

    gc.collect()



print('*' * 120)

print('Mean smape in each fold:', smape_score.mean(1))

print()

print('Mean smape of each property:', smape_score.mean(0))

print()

print('Final mean smape:', np.mean(smape_score))

for i in range(len(imps)):

    # 对五折的 importance score 进行平均

    imps[i]['score_mean'] = imps[i].apply(lambda row: (row['score0'] + row['score1'] + row['score2'] + row['score3'] + \

                            row['score4']) / 5, axis=1)

    imps[i] = imps[i].sort_values(by='score_mean', ascending=False).reset_index(drop=True)

    

for imp in imps:

    print(imp.head())

print()

# 每个 property 重要性前四的特征

high_imps_list = []

for imp in imps:

    high_imps_list.append(imp['feat'].values[:4])



'''

使用原始 3177 特征跑出来的 high_imps_list 如下

[array(['25', '3175', '34', '3158'], dtype=object),

 array(['25', '16', '3175', '3157'], dtype=object),

 array(['34', '10', '3157', '15'], dtype=object),

 array(['3161', '3155', '3173', '3157'], dtype=object),

 array(['3155', '3173', '3175', '3151'], dtype=object),

 array(['3155', '3173', '3175', '3151'], dtype=object)]

'''

high_imps_list

pred_df = pd.DataFrame(pred_test, columns=[f'p{i}' for i in range(1, 7)])



result = pd.concat([test_data[['id']], pred_df], axis=1).reset_index(drop=True)



result.to_csv('result.csv', index=False)

result.head()