大数据文摘出品
高能质子对撞中会产生大量粒子团喷注（jet），喷注可以根据其不同内在特性分为胶体喷注、轻夸克喷注、魅夸克喷注、美夸克喷注。
但当前缺乏可靠方法分类所测量的真实喷注，开发一种稳健的算法来识别喷注味道，有助于学界更直接地比较实验观测和基本粒子理论。
目前，由北京智源人工智能研究院和 biendata 共同发布的“高能对撞粒子分类挑战赛”（2019 年 11 月 - 2020 年 3 月）正在火热进行中，总奖金为 10 万元。比赛要求选手要求选手根据喷注的性质（如喷注所含的粒子数、喷注能量、喷注质量、喷注方向），以及喷注中所有粒子的特征和对应的碰撞事件，把喷注分成四类中的一类。比赛和数据可于下方链接查看，或点击“阅读原文”，欢迎所有感兴趣的报名读者参赛。
目前，赛事已接近半程，为帮助选手进一步提升模型预测结果，biendata 邀请排行榜前列的CChan、FastCloud、hengheng选手分享三个分数在0.75左右的baseline，希望可以为大家在数据探索、预处理、模型选择、参数调优等方面提供参考和思路。
比赛地址：
https://www.biendata.com/competition/jet/
Baseline 1：
https://www.biendata.com/models/detail/4002/L_notebook/
Baseline 2：
https://www.biendata.com/models/category/4067/L_notebook/
Baseline 3：
https://biendata.com/models/category/4273/L_notebook/
Baseline 概述
（一）选手CChan基于数据标签的分布和特征选择直接对碰撞事件进行分类，以减少训练的数据量，提高效率。在特征工程方面，CChan的方案主要包括3类特征：对事件包含的喷注的物理属性求统计值、对事件包含的粒子的物理属性求统计值、对喷注包含的粒子的物理属性求统计值，再对事件求统计值。 在模型方面，CChan 使用了选用了 CatBoost 模型，五折交叉验证，该 baseline 线下验证分数为 0.765，线上分数为 0.766。
具体代码实现欢迎到页面查看或浏览下文：
https://www.biendata.com/models/detail/4002/L_notebook/
（二）选手FastCloud首先提取出数据集中 min、max、mean、count 等统计特征。具体而言：对于 jet 数据,就 event_id 做 groupby，对数值特征抽取 min、max 等特征，以 event_id 为键做 merge；对于 particle数据，就 event_id 做 groupby，对数值特征抽取 min、max 等特征，对类别特征抽取 count、nunique 特征，以event_id 为键做 merge；对于 particle 数据，就 jet_id 做 groupby，对数值特征抽取 min、max 等特征，对类别特征抽取 count、nunique 特征，以jet_id为键作 merge。 在模型选择上，FastCloud 使用 xgboost 模型，5 折交叉验证，线上分数为 0.747，
具体代码实现欢迎访问页面：
https://www.biendata.com/models/category/4067/L_notebook/
（三）选手hengheng首先对数据进行预处理和格式转化，并通过减小精度来减小内存使用，在读取文件后将训练集和测试集合并。在特征工程上，进行简单的特征线性组合,将 all_jet 文件中的除 xyz 空间特征外的其他特征分配(相除)到 x、y、z 方向上，在 all_jet 文件中以 event_id 为主键分组做相应的统计特征；对 all_particle 文件进行简单的特征线性组合，以 jet_id 为主键进行分组得到每个 particle 的统计特征。然后划分数据集。并将 label 进行 one-hot 处理，最后建立模型和参数设置。 选手 hengheng 使用的模型是 catboost，他认为大多数情况下 catboost 效果是三个 boost（xgb,lgb,cat）中最好的，且需要设置的参数较少并支持 GPU，该 baseline 线上得分 0.73、线下得分 0.75。
具体代码实现欢迎访问页面：
https://biendata.com/models/category/4273/L_notebook/
Baseline 详情
CChan 版
本 baseline 线下验证分数为0.765，线上分数为0.766。
1.赛题引入与问题分析
宇宙中大多数物质由原子构成，原子又由原子核和电子组成。其中，电子是基本粒子，但原子核又可分为质子和中子，并可进一步分为夸克和胶子。通过高能粒子碰撞，可以产生包括夸克和中子在内的大量粒子（particle）。在这个过程（碰撞事件，event）中，产生的呈喷射状的粒子团被称为喷注（jet）。根据喷注的不同内在特性（如质量和色量子数），可以将喷注划分为四类，包括：1）胶体喷注，2）轻夸克喷注，3）魅夸克喷注，4）美夸克喷注。
本次比赛要求选手根据碰撞事件中产生的喷注及粒子的物理属性，对喷注类型进行分类。提供的数据包括了事件、喷注、粒子三组数据，评价指标为ROC AUC。
对于多分类问题，AUC的求取是先用one-vs-all方法将多分类问题转换为 4 个二分类问题（对于其中任意一类，我们可以将其他类别归为另一大类），再对每个二分类问题求AUC，最后取其算术平均数。sklearn库中的roc_auc_score函数支持多分类计算，但需要注意的是，本次比赛评价指标是用预测标签求AUC，因此线下验证时，应该先将概率转为标签再进行验证。
2.数据探索性分析和可视化
2.1 数据预览
比赛数据包括了三组：
- 事件： 碰撞事件id，包含的喷注数
- 喷注： 喷注id，包含的粒子数，喷注的物理属性，所属的事件id，类别标签
- 粒子： 粒子类别，粒子的物理属性，所属的喷注id
- 物理属性包括： 能量、质量、方向（x，y，z）
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
# 事件
train_event = pd.read_csv('../jet_complex_data/complex_train_R04_event.csv')
test_event = pd.read_csv('../jet_complex_data/complex_test_R04_event.csv')
# 喷注
train_jet = pd.read_csv('../jet_complex_data/complex_train_R04_jet.csv')
test_jet = pd.read_csv('../jet_complex_data/complex_test_R04_jet.csv')
# 粒子
train_particle = pd.read_csv('../jet_complex_data/complex_train_R04_particle.csv')
test_particle = pd.read_csv('../jet_complex_data/complex_test_R04_particle.csv')
train_event.head(3)
train_jet.head(3)
train_particle.head(3)2.2 数据组成
训练数据包括：371377个碰撞事件，1134555个喷注，24297352个粒子。
测试数据包括：176720个碰撞事件，537949个喷注，11493800个粒子。
事件、喷注数据文件都以其id作为主键，事件、喷注文件中每一条数据代表不同的事件、喷注。因此我们可以利用其id直接进行merge操作。
此外，数据十分完整，没有缺失值，并且每个事件都有喷注数据，每个喷注都有粒子数据。
代码语言：javascript复制print('--- Data Size')
print('event: train %d, test %d' % (len(train_event), len(test_event)))
print('jet: train %d, test %d' % (len(train_jet), len(test_jet)))
print('particle: train %d, test %d' % (len(train_particle), len(test_particle)))
print('--- Amount')
print('event: train %d, test %d' % (train_event.event_id.nunique(), test_event.event_id.nunique()))
print('jet: train %d, test %d' % (train_jet.jet_id.nunique(), test_jet.jet_id.nunique()))
print('event in jet: train %d, test %d' % (train_jet.event_id.nunique(), test_jet.event_id.nunique()))
print('jet in particle: train %d, test %d' % (train_particle.jet_id.nunique(), test_particle.jet_id.nunique()))
print('--- NaN')
print('event: train %d, test %d' % (train_event.isnull().sum().sum(), test_event.isnull().sum().sum()))
print('jet: train %d, test %d' % (train_jet.isnull().sum().sum(), test_jet.isnull().sum().sum()))
print('particle: train %d, test %d' % (train_particle.isnull().sum().sum(), test_particle.isnull().sum().sum()))
--- Data Size
event: train 371377, test 176720
jet: train 1134555, test 537949
particle: train 24297352, test 11493800
--- Amount
event: train 371377, test 176720
jet: train 1134555, test 537949
event in jet: train 371377, test 176720
jet in particle: train 1134555, test 537949
--- NaN
event: train 0, test 0
jet: train 0, test 0
particle: train 0, test 02.3 问题标签
问题标签已经经过匿名处理，分别为 21，1，4，5。喷注标签分布较为均匀，其中喷注类别21的数量最多。
赛题虽然要求对喷注进行分类，但是我们可以发现，在这份数据中，同一碰撞事件中喷注的类别都是相同的，因此我们也可以将这个分类问题转化为对碰撞事件的分类。
可以看到，事件标签的分布十分均匀。
代码语言：javascript复制ax = train_jet.label.value_counts(normalize=True).plot(kind='bar', title='Distribution of jet label')
event_label = train_jet.groupby('event_id')['label'].agg('nunique')
print('Max number of jet types in a event: ', event_label.max())
Max number of jet types in a event: 1
ax = train_jet.groupby('event_id')['label'].nth(0).value_counts(normalize=True).plot(kind='bar',
title='Distribution of event label')
2.4 物理属性
2.4.1 喷注属性
通过绘制喷注不同属性的箱线图，可以看出类别21在x方向和能量等属性上，与其他类别的区分度较大。
代码语言：javascript复制plt.figure(figsize=(12, 4))
ax = plt.subplot(1,2,1)
ax.set_title('boxplot: jet direction')
sns.boxplot(data=train_jet[['jet_px', 'jet_py', 'jet_pz']])
plt.subplot(1,2,2)
ax = sns.boxplot(y='jet_px',x='label',data=train_jet)
ax.set_title('boxplot: jet_x of different jet type')
plt.show()
plt.figure(figsize=(12, 4))
ax = plt.subplot(1,2,1)
ax.set_title('boxplot: jet energy & mass')
sns.boxplot(data=train_jet[['jet_energy', 'jet_mass']])
plt.subplot(1,2,2)
ax = sns.boxplot(y='jet_energy',x='label',data=train_jet)
ax.set_title('boxplot: jet_energy of different jet type')
plt.show()
2.4.2 粒子属性
粒子的类型包含14种，类型分布和质量如下。
粒子类型22的数量最多，其质量为0，因此可以判断为胶子。
代码语言：javascript复制ax = train_particle.particle_category.value_counts(normalize=True).plot(kind='bar', title='Distribution of particle type')
print('粒子质量统计值')
train_particle.groupby('particle_category')['particle_mass'].agg(['min', 'max', 'mean'])
粒子质量统计值
3.思路分析
根据前文分析可知，同一碰撞事件中喷注的类别都是相同的，因此我们可以 1）构造喷注的特征，对喷注进行分类，或者 2）构造事件的特征，对事件进行分类。本baseline将采取方案二，直接对事件进行分类。这样做的好处是减少训练的数据量，提高效率，因为事件数量仅是喷注数量的三分之一。
特征工程方面，主要包括3类特征：
- 对事件包含的喷注的物理属性求统计值。
- 对事件包含的粒子的物理属性求统计值。
- 对喷注包含的粒子的物理属性求统计值，再对事件求统计值。
模型可以选择支持GPU的机器学习模型，提高训练效率。本baseline选用了CatBoost模型，五折交叉验证。
4.代码与代码解释
代码语言：javascript复制import warnings
warnings.filterwarnings('ignore')
import os
import numpy as np
import pandas as pd
from scipy import stats
import pickle
import time
import catboost as cbt
from sklearn.preprocessing import LabelEncoder, LabelBinarizer
from sklearn.metrics import roc_auc_score, accuracy_score, f1_score
from sklearn.model_selection import StratifiedKFold4.1 读取数据
代码语言：javascript复制train_jet = pd.read_csv('../jet_complex_data/complex_train_R04_jet.csv')
test_jet = pd.read_csv('../jet_complex_data/complex_test_R04_jet.csv')
train_event = pd.read_csv('../jet_complex_data/complex_train_R04_event.csv')
test_event = pd.read_csv('../jet_complex_data/complex_test_R04_event.csv')
train_particle = pd.read_csv('../jet_complex_data/complex_train_R04_particle.csv')
test_particle = pd.read_csv('../jet_complex_data/complex_test_R04_particle.csv')4.2 构造标签和特征
代码语言：javascript复制# 构造事件标签
def gen_event_label(event, jet):
assert jet.groupby('event_id')['label'].nunique().max() == 1
event_label = jet.groupby('event_id')['label'].first().reset_index()
event_label = event[['event_id']].merge(event_label, 'left', 'event_id')
return event_label
# 特征：1）对事件包含的喷注的物理属性求统计值。
def gen_jet_feat(event, jet):
feat = jet[['event_id', 'jet_px', 'jet_py', 'jet_pz', 'jet_energy', 'jet_mass', 'number_of_particles_in_this_jet']]
feat['jet_p'] = (feat['jet_px'] ** 2 feat['jet_py'] ** 2 feat['jet_pz'] ** 2) ** 0.5
feat['jet_cos(x)'] = feat['jet_px'] / feat['jet_p']
feat['jet_cos(y)'] = feat['jet_py'] / feat['jet_p']
feat['jet_cos(z)'] = feat['jet_pz'] / feat['jet_p']
feat['jet_angle(x)'] = np.arccos(feat['jet_cos(x)'])
feat['jet_angle(y)'] = np.arccos(feat['jet_cos(y)'])
feat['jet_angle(z)'] = np.arccos(feat['jet_cos(z)'])
feat['jet_energy/jet_mass'] = feat['jet_energy'] / feat['jet_mass']
cols = ['jet_px', 'jet_py', 'jet_pz', 'jet_energy', 'jet_mass', 'number_of_particles_in_this_jet', 'jet_energy/jet_mass',
'jet_p', 'jet_cos(x)', 'jet_cos(y)', 'jet_cos(z)', 'jet_angle(x)', 'jet_angle(y)', 'jet_angle(z)']
st = ['min', 'max', 'mean', 'std', 'sum']
st_cols = [(c '_' s) for c in cols for s in st]
feat = feat.groupby('event_id')[cols].agg(st).reset_index()
feat.columns = ['event_id'] st_cols
feat = event[['event_id', 'number_of_jet_in_this_event']].merge(feat, 'left', 'event_id')
feat = feat.drop(columns=['event_id'])
return feat
# 特征：2）对事件包含的粒子的物理属性求统计值。
def gen_particle_feat(event, jet, particle):
# 事件的粒子属性特征
particle = particle.copy().merge(jet[['jet_id', 'event_id']], 'left', 'jet_id')
particle['particle_p'] = (particle['particle_px'] ** 2 particle['particle_py'] ** 2 particle['particle_pz'] ** 2) ** 0.5
particle['particle_cos(x)'] = particle['particle_px'] / particle['particle_p']
particle['particle_cos(y)'] = particle['particle_py'] / particle['particle_p']
particle['particle_cos(z)'] = particle['particle_pz'] / particle['particle_p']
particle['particle_angle(x)'] = np.arccos(particle['particle_cos(x)'])
particle['particle_angle(y)'] = np.arccos(particle['particle_cos(y)'])
particle['particle_angle(z)'] = np.arccos(particle['particle_cos(z)'])
particle['particle_energy/particle_mass'] = particle['particle_energy'] / particle['particle_mass']
cols = ['particle_px', 'particle_py', 'particle_pz', 'particle_energy', 'particle_mass',
'particle_p', 'particle_cos(x)', 'particle_cos(y)', 'particle_cos(z)',
'particle_angle(x)', 'particle_angle(y)', 'particle_angle(z)', 'particle_energy/particle_mass']
st = ['min', 'max', 'mean', 'std', 'sum']
st_cols = [(c '_e_' s) for c in cols for s in st]
particle_st = particle.groupby('event_id')[cols].agg(st).reset_index()
particle_st.columns = ['event_id'] st_cols
feat = event[['event_id']].merge(particle_st, 'left', 'event_id')
# 事件的粒子类别特征
particle['1'] = 1
cat_cnt = particle.pivot_table(index='event_id', columns='particle_category', values='1', aggfunc='sum').reset_index()
cat_cols = ['cat_e_%d' % i for i in range(14)]
cat_cnt.columns = ['event_id'] ['cat_e_%d' % i for i in range(14)]
cat_cnt['cat_e_sum'] = cat_cnt[cat_cols].sum(axis=1)
cat_cnt_rate = cat_cnt[cat_cols] / cat_cnt['cat_e_sum'].values.reshape((-1, 1))
cat_cnt_rate.columns = ['%s_e_rate' % c for c in cat_cols]
cat_cnt = pd.concat([cat_cnt, cat_cnt_rate], axis=1)
feat = feat.merge(cat_cnt, 'left', 'event_id')
feat = feat.drop(columns=['event_id'])
return feat
# 特征：3）对喷注包含的粒子的物理属性求统计值，再对事件求统计值。
def gen_jet_particle_feat(event, jet, particle):
# 喷注的粒子属性特征
particle = particle.copy() # merge(jet[['jet_id', 'event_id']], 'left', 'jet_id')
particle['particle_p'] = (particle['particle_px'] ** 2 particle['particle_py'] ** 2 particle['particle_pz'] ** 2) ** 0.5
particle['particle_cos(x)'] = particle['particle_px'] / particle['particle_p']
particle['particle_cos(y)'] = particle['particle_py'] / particle['particle_p']
particle['particle_cos(z)'] = particle['particle_pz'] / particle['particle_p']
particle['particle_angle(x)'] = np.arccos(particle['particle_cos(x)'])
particle['particle_angle(y)'] = np.arccos(particle['particle_cos(y)'])
particle['particle_angle(z)'] = np.arccos(particle['particle_cos(z)'])
particle['particle_energy/particle_mass'] = particle['particle_energy'] / particle['particle_mass']
cols = ['particle_px', 'particle_py', 'particle_pz', 'particle_energy', 'particle_mass',
'particle_p', 'particle_cos(x)', 'particle_cos(y)', 'particle_cos(z)',
'particle_angle(x)', 'particle_angle(y)', 'particle_angle(z)', 'particle_energy/particle_mass']
st = ['min', 'max', 'mean', 'std', 'sum']
st_cols = [(c '_j_' s) for c in cols for s in st]
particle_st = particle.groupby('jet_id')[cols].agg(st).reset_index()
particle_st.columns = ['jet_id'] st_cols
# 喷注的粒子类别特征
particle['1'] = 1
cat_cnt = particle.pivot_table(index='jet_id', columns='particle_category', values='1', aggfunc='sum').reset_index()
cat_cols = ['cat_j_%d' % i for i in range(14)]
cat_cnt.columns = ['jet_id'] ['cat_j_%d' % i for i in range(14)]
cat_cnt['cat_j_sum'] = cat_cnt[cat_cols].sum(axis=1)
cat_cnt_rate = cat_cnt[cat_cols] / cat_cnt['cat_j_sum'].values.reshape((-1, 1))
cat_cnt_rate.columns = ['%s_j_rate' % c for c in cat_cols]
cat_cnt = pd.concat([cat_cnt, cat_cnt_rate], axis=1)
feat = jet[['event_id', 'jet_id']].merge(particle_st, 'left', 'jet_id').merge(cat_cnt, 'left', 'jet_id')
# 对事件求上述喷注特征的统计值
cols = [c for c in feat.columns if not c in ['event_id', 'jet_id']]
st = ['min', 'max', 'mean', 'std', 'sum']
st_cols = [(c '_e_' s) for c in cols for s in st]
feat = feat.groupby('event_id')[cols].agg(st).reset_index()
feat.columns = ['event_id'] st_cols
feat = event[['event_id']].merge(feat, 'left', 'event_id')
feat = feat.drop(columns=['event_id'])
return feat
# 事件标签
event_label = gen_event_label(train_event, train_jet)
# 特征 1
train_jet_feat = gen_jet_feat(train_event, train_jet)
test_jet_feat = gen_jet_feat(test_event, test_jet)
# 特征 2
train_particle_feat = gen_particle_feat(train_event, train_jet, train_particle)
test_particle_feat = gen_particle_feat(test_event, test_jet, test_particle)
# 特征 3
train_jet_particle_feat = gen_jet_particle_feat(train_event, train_jet, train_particle)
test_jet_particle_feat = gen_jet_particle_feat(test_event, test_jet, test_particle)
# 合并特征
train_data = pd.concat([train_jet_feat, train_particle_feat, train_jet_particle_feat], axis=1)
test_data = pd.concat([test_jet_feat, test_particle_feat, test_jet_particle_feat], axis=1)
used_feats = list(train_data.columns)
x_train = train_data
x_test = test_data
print(len(used_feats))
print(used_feats)
print(x_train.shape)
print(x_test.shape)
label = 'label'
lbl_enc = LabelEncoder()
event_label[label] = lbl_enc.fit_transform(event_label[label])
print(lbl_enc.classes_)
print(lbl_enc.transform(lbl_enc.classes_))
y_train = event_label[label]
del train_jet_feat, train_particle_feat, train_jet_particle_feat
del test_jet_feat, test_particle_feat, test_jet_particle_feat4.3 模型训练预测
代码语言：javascript复制lbr = LabelBinarizer().fit(y_train)
def auc_metric(y_true, y_pred):
y_true = lbr.transform(y_true)
y_pred = y_pred.reshape((4, -1)).T
y_pred = lbr.transform(np.argmax(y_pred, axis=1))
score = roc_auc_score(y_true=y_true, y_score=y_pred, average='macro')
return 'auc', score, Truetic = time.time()
test_pred = np.zeros((len(x_test), 4))
oof_pred = np.zeros((len(x_train), 4))
feat_imp = pd.DataFrame(used_feats, columns=['feat'])
feat_imp['imp'] = 0
scores = []
kfold = StratifiedKFold(n_splits=5, random_state=12306, shuffle=True)
for i, (trn_idx, val_idx) in enumerate(kfold.split(X=x_train, y=y_train)):
print('-' * 88)
print('Fold %d:' % i)
x_trn, y_trn = x_train.iloc[trn_idx], y_train.iloc[trn_idx]
x_val, y_val = x_train.iloc[val_idx], y_train.iloc[val_idx]
trn_pool = cbt.Pool(x_trn, y_trn)
val_pool = cbt.Pool(x_val, y_val)
model = cbt.CatBoostClassifier(iterations=100000, learning_rate=0.1, eval_metric='MultiClass',# depth=10,
use_best_model=True, random_seed=2020, logging_level='Verbose',
task_type='GPU', devices='0', early_stopping_rounds=200, loss_function='MultiClass',
)
# model.set_params(**params)
model.fit(trn_pool, eval_set=val_pool, verbose=100)
pickle.dump(file=open('./models/cbt_model_%d.pkl' % i, 'wb'), obj=model)
feat_imp['imp'] = (model.feature_importances_ / 5)
scores.append(model.best_score_['validation']['MultiClass'])
test_pred = (model.predict_proba(x_test) / 5)
oof_pred[val_idx] = model.predict_proba(x_val)
del x_trn, y_trn, x_val, y_val
del trn_pool, val_pool
toc = time.time()
print('times: %f' % (toc - tic))print('loss：%s' % scores)
print('mean loss: %f' % np.mean(scores))
print('acc: %f' % accuracy_score(y_train, np.argmax(oof_pred, axis=1)))
print('auc: %f' % roc_auc_score(y_true=lbr.transform(y_train),
y_score=lbr.transform(np.argmax(oof_pred, axis=1)),
average='macro'))# 以喷注为主体重新计算AUC
event_pred = train_event[['event_id']]
event_pred['pred'] = np.argmax(oof_pred, axis=1)
jet_pred = train_jet[['event_id']].merge(event_pred, 'left', 'event_id')['pred'].values
jet_label = lbl_enc.transform(train_jet['label'])
print('acc: %f' % accuracy_score(jet_label, jet_pred))
print('auc: %f' % roc_auc_score(y_true=lbr.transform(jet_label),
y_score=lbr.transform(jet_pred),
average='macro'))
# 0.766特征重要性
代码语言：javascript复制import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=(5,10))
feat_imp= feat_imp.sort_values(by='imp', ascending=True)[-50:]
plt.barh(feat_imp['feat'], feat_imp['imp'])4.4 生成结果
代码语言：javascript复制event_pred = test_event[['event_id']]
event_pred['pred'] = np.argmax(test_pred, axis=1)
jet_pred = test_jet[['event_id']].merge(event_pred, 'left', 'event_id')['pred'].valuessubmission = pd.DataFrame()
submission['id'] = test_jet['jet_id']
submission['label'] = lbl_enc.inverse_transform(jet_pred)
submission.head()
submission.to_csv('./result.csv', index=False)FastCloud版
由于文章篇幅限制，该baseline请查看：
https://www.biendata.com/models/category/4067/L_notebook/
hengheng 版
由于文章篇幅限制，该baseline请查看：
https://biendata.com/models/category/4273/L_notebook/
