【CCF BDCI】2020 滴滴路况状态时空预测 冠军方案详解

参考信息

1. 比赛平台链接：https://www.datafountain.cn/competitions/466
2. Datafountain平台并未保留原始数据且滴滴的数据外链已失效，多方检索找到的原始数据集，源自博主https://blog.csdn.net/tangxianyu，感谢！
   比赛数据集链接：https://pan.baidu.com/s/1Zwd1JK3sf_szCGykEbZPCQ
   提取码：xcul
3. 冠军方案源码链接：https://github.com/shyoulala/CCF_BDCI_2020_DIDI_rank1_solution
   冠军方案阐述公众号：OTTO Data Lab
   公众号方案文章链接：https://mp.weixin.qq.com/s/h2rn2t7T79mFkczuYOVHBw

方案详解

基本信息

1. 任务：
   比赛提供滴滴平台2019年7月1日至2019年7月30日西安市的实时和历史路况信息, 以及道路属性和路网拓扑信息，预测未来一段时间内道路小段的路况状态(即畅通, 缓行和拥堵几类状态）

2. 评估指标
   采用加权 F1 Score 作为算法评价指标。其中畅通权重0.2, 缓行权重0.2, 拥堵权重0.6

3. 数据集介绍
   Datafountain官网的数据集介绍和真实数据对应不太好，忽略；直接从实际数据进行理解介绍。

   Link: 对完整道路按照拓扑切分后得到的小段, 由唯一id标识
   时间片: 对时间的离散化描述,以2分钟为一个单位,2分钟内认为道路的路况状态是统一的。

   1）道路属性特征表
   其中 length、width、speed_limit是连续特征，其余为类别特征

   link    length    direction    road_class    speed_class   LaneNum        speed_limit    level    width
   0        19        1            5            7                1            4.166667        5        30
   1        19        1            5            7                1            4.166667        5        30
   

   2）路网拓扑特征表
   道路Link和它的下游道路Link

   link           next_link
   611897	    630844,611898,611691
   

   3）历史与实时路况
   以第一行数据为例

   link    label   当前时间片    未来时间片
   353495    1        236           245;
   
   # 近期数据
   时间片  速度  eta速度 标签状态 车辆数
   232:   29.80, 32.40,    1,     4       233:31.60,32.20,1,2 234:20.00,21.90,2,2 235:22.20,25.90,2,5 236:21.30,26.30,2,4;
   
   # 28天前数据
   245:30.00,32.70,0,9 246:30.00,36.10,0,10 247:27.40,35.20,1,12 248:26.90,35.70,1,10 249:28.90,37.00,1,9;
   # 21天前
   245:36.10,37.30,1,7 246:29.30,38.50,1,7 247:27.70,39.70,1,6 248:28.60,40.20,1,3 249:29.60,38.70,1,4;
   # 14天前
   245:30.40,40.10,1,6 246:32.30,40.10,1,6 247:30.60,41.10,1,5 248:29.60,39.20,1,4 249:28.00,37.90,1,4;
   # 7天前
   245:28.30,38.40,1,7 246:28.20,39.40,1,6 247:28.80,35.10,1,3 248:30.00,35.60,1,4 249:29.40,37.20,1,5
   

数据分析

1. 道路属性特征拥有最全的Link为686105个，训练集为14523 ， 测试集为12224 ， 路网拓扑含下游节点为684813 ；
   训练集和测试集Link为交集状态 ，且两者有少量节点没有路网拓扑信息。
2. 冠军方案经数据分析发现：

   • 训练集中只有0.07%不存在当前时间片最近的信息，而测试集存在11%的数据没有当前时间片的信息，所以测试集分布不一致，这就导致预测会更具有难度，因为对于路况预测来讲，待预测时间片近期时间片的信息非常重要。

   • 对于训练集的数据而言，两分钟一个时刻，训练数据未来时间片分布如下，而测试集的数据却集中在凌晨，分布不一样。
     

   • 待预测时间-当前时间越小，所提供的信息的价值就越大，从训练集测试集分布来看，训练集趋于平稳，而测试集集中大于20分钟的情况，这说明测试集预测更加困难，最终线下线上差距会比较大，若使用传统特征工程+树模型，需要对原始训练数据进行采样构建。(图为时间差占各自的占比)
     

LGB方案

采样方案

基于上述分析，在使用树模型时需要保持训练集和测试集的分布一致。

• 根据测试集待预测时间分布，选取训练集5-30日（5日之前的数据用来做target encoding，防止泄露）待预测时间片在10-40内的数据和30日待预测时间片大于40的数据5万；
  上文数分提到，测试集的未来时间片集中在凌晨部分，所以主要采样数据集中在(10,40)

if d != 20190730:
      tmp_file = tmp_file[(tmp_file['future_slice_id'] > 10) & (tmp_file['future_slice_id'] < 40)]
  else:
      tmp_file_1 = tmp_file[(tmp_file['future_slice_id'] > 10) & (tmp_file['future_slice_id'] < 40)]
      tmp_file_2 = tmp_file[tmp_file['future_slice_id'] >= 40].sample(n=50000, random_state=1)
      tmp_file = tmp_file_1.append(tmp_file_2)
            
# 猜测为了符合测试集的分布，该方案先从历史数据中获得凌晨部分的数据；其次从最接近测试集时间的0730当天抽取其余长尾的数据

• 根据测试集分布预测间隔分布采样，得到训练集D，后续特征工程主要在D上进行

# 根据测试集curr_state 即当前时间片的状态分布情况 从训练集种进行相应的采样
# {1: 328063, 0: 78498, 2: 5863, 3: 1149, 4: 166}
curr_state_dict = (origin_test['curr_state'].value_counts(normalize=True) * origin_train.shape[0]).astype(
    int).to_dict()
sample_train = pd.DataFrame()
for t, group in origin_train.groupby('curr_state'):
    if t == 0:
        sample_tmp = group
    else:
        sample_tmp = group.sample(n=curr_state_dict[t], random_state=1)
    sample_train = sample_train.append(sample_tmp)

# 根据time_diff分布采样
sample_df = pd.DataFrame()
for t, group in sample_train.groupby('time_diff'):
    sample_tmp = group.sample(n=time_dict[t], random_state=1)
    sample_df = sample_df.append(sample_tmp)

for j, i in enumerate(range(11, 16)):
    sample_df.loc[sample_df['time_diff'] == i, 'time_diff'] = sample_df.loc[
                                                                  sample_df['time_diff'] == i, 'time_diff'] * (
                                                                      10 + j)
# 将测试集突出的time_diff部分--11-15 进行乘积变大操作
# 这样在升排序时，由于time_diff 值越大，越靠后，再通过去重降采样时，会倾向保留time_diff值最大的那部分
sample_df = sample_df.sort_values('time_diff').drop_duplicates(subset=['linkid', 'future_slice_id'], keep='last')
# time_diff 还原
for j, i in enumerate(range(11, 16)):
    sample_df.loc[sample_df['time_diff'] == i * (10 + j), 'time_diff'] = sample_df.loc[
                                                                             sample_df['time_diff'] == i * (
                                                                                     10 + j), 'time_diff'] / (
                                                                                 10 + j)

特征工程方案

整体特征工程方案如图所示

• 当前和历史状态的统计特征，包含速度、eta速度、车辆数、状态

 numeric_list = ['recent_speed', 'recent_eta', 'recent_vichles_num',
                 'history_speed_28', 'history_speed_21', 'history_speed_14', 'history_speed_7',
                 'history_eta_28', 'history_eta_21', 'history_eta_14', 'history_eta_7',
                 'history_vichles_num_28', 'history_vichles_num_21', 'history_vichles_num_14','history_vichles_num_7'
                 ]
 for col in tqdm(numeric_list):
     data[col + '_mean'] = data[col].apply(lambda x: np.mean(x))
     data[col + '_max'] = data[col].apply(lambda x: np.max(x))
     data[col + '_min'] = data[col].apply(lambda x: np.min(x))
     data[col + '_std'] = data[col].apply(lambda x: np.std(x))
     data[col + '_median'] = data[col].apply(lambda x: np.median(x))

# status 取值有 0 1 2 3 4
status_list = ['recent_status', 'history_status_28', 'history_status_21', 'history_status_14', 'history_status_7']
for col in tqdm(status_list):
    data[col + '_mean'] = data[col].apply(lambda x: np.mean(x))
    data[col + '_max'] = data[col].apply(lambda x: np.max(x))
    data[col + '_std'] = data[col].apply(lambda x: np.std(x))
    # before 特征
    for i in range(0, 5):
        data[col + '_bef{}'.format(i + 1)] = data[col].apply(lambda x: x[i])
    # 状态计数特征
    for i in range(0, 5):
        data[col + '_{}_cnt'.format(i)] = data[col].apply(lambda x: x.count(i))

• 目标编码特征
  方案计算了link_id 、 future_slice_id 、 link_id + hour 、 link_id + hour + weekday 四种形式的编码特征
  该部分核心逻辑如下：

def gen_ctr_features(d, key):
	# 目标日期前的数据
    his_ = his_df[his_df['day'] < d].copy()
    his_ = his_.drop_duplicates(subset=['link_id', 'future_slice_id', 'day', 'label'], keep='last')
    
    # 标签的独热编码和源数据拼接
    dummy = pd.get_dummies(his_['label'], prefix='label')
    his_ = pd.concat([his_, dummy], axis=1)
    
    # ctr即目标编码的计算方式为，针对key进行分组，统计不同标签的平均值
    ctr = his_.groupby(key)[['label_0', 'label_1', 'label_2']].mean().reset_index()
    ctr.columns = key + ['_'.join(key) + '_label_0_ctr', '_'.join(key) + '_label_1_ctr', '_'.join(key) + '_label_2_rt']
    ctr['day'] = d
    del his_
    gc.collect()
    return ctr

• 近期数据中，不同标签对应时间片的特征；即在近期数据中，发生标签1的最大、最小时间片以及发生标签1的时间片之间的时间差值得平均值。
  具体核心逻辑分别为：

# time_slice是link近期状态中时间片的list , 例如 [16, 17, 18, 19, 20]
# status_list是近期状态的list , 例如 [1.0, 1.0, 1.0, 1.0, 1.0]
# status 是目标的标签 label的取值是 1、2、3、4
def gen_label_timedelta(time_slice, status_list, status):
    timedelta = [time_slice[i] for i, f in enumerate(status_list) if f in status]
    if len(timedelta) > 0:
        timedelta = np.max(timedelta)
    else:
        timedelta = np.nan
    return timedelta


def gen_label_timedelta_min(time_slice, status_list, status):
    timedelta = [time_slice[i] for i, f in enumerate(status_list) if f in status]
    if len(timedelta) > 0:
        timedelta = np.min(timedelta)
    else:
        timedelta = np.nan
    return timedelta

# 这里使用np.diff计算相邻元素值之间得差值，元素至少为1
def gen_label_timedelta_diff(time_slice, status_list, status):
    timedelta = [time_slice[i] for i, f in enumerate(status_list) if f in status]
    if len(timedelta) > 1:
        timedelta = np.mean(np.diff(timedelta))
    else:
        timedelta = np.nan
    return timedelta

• 该部分特征工程，使用路网拓扑数据，主要实现每个节点上\下游节点状态的统计计算
  核心逻辑如下：

# 上下游节点的路况信息
def get_topo_info(df, topo_df, slices=30, mode='down'):
    if mode == 'down':
        flg = 'down_target_state'
    else:
        flg = 'up_target_state'

    # 筛选topo_df中属于df linkid的数据，且筛选target_link_list中属于df linkid的数据 保证至少有一个节点
    use_ids = set(df['linkid'].unique())
    topo_df = topo_df[topo_df['linkid'].isin(use_ids)]
    topo_df['target_link_list'] = topo_df['target_link_list'].apply(
        lambda x: [int(c) for c in x.split(',') if int(c) in use_ids])
    topo_df['len'] = topo_df['target_link_list'].apply(len)
    topo_df = topo_df[topo_df['len'] > 0]
    del topo_df['len']

    # 将列表展开
    # 以统计下游节点模式为例
    # 先匹配linkid的future_slice_id ， 然后匹配其下游节点的current_slice_id和curr_state
    curr_df = []
    for i in topo_df.values:
        for j in i[1]:
            curr_df.append([i[0], j])
    curr_df = pd.DataFrame(curr_df, columns=['linkid', 'target_id'])
    curr_df = curr_df.merge(df[['linkid', 'future_slice_id']], on='linkid', how='left')

    tmp_df = df[['linkid', 'current_slice_id', 'curr_state']]
    tmp_df.columns = ['target_id', 'current_slice_id', 'curr_state']

    curr_df = curr_df.merge(tmp_df, on='target_id', how='left')

    # linkid需要预测的future_slice_id - 下游节点的current_slice_id
    # 这个差值控制在[0,slices]
    curr_df['{}_diff_slice'.format(flg)] = curr_df['future_slice_id'] - curr_df['current_slice_id']
    curr_df = curr_df[(curr_df['{}_diff_slice'.format(flg)] >= 0) & (curr_df['{}_diff_slice'.format(flg)] <= slices)]

    curr_df = curr_df.drop_duplicates()
    tmp_list = ['{}_diff_slice'.format(flg)]
    # 权重 根据下游节点当前时间片距离future_slice_id的远近赋予权重
    curr_df['{}_diff_slice'.format(flg)] = (slices - curr_df['{}_diff_slice'.format(flg)]) / slices
    # 4 个标签 0、1、2、3、4
    for s in range(5):
        curr_df['{}_{}'.format(flg, s)] = curr_df['curr_state'].apply(lambda x: 1 if x == s else 0)
        curr_df['{}_{}'.format(flg, s)] = curr_df['{}_{}'.format(flg, s)] * curr_df['{}_diff_slice'.format(flg)]
        tmp_list.append('{}_{}'.format(flg, s))
    curr_df = curr_df.groupby(['linkid', 'future_slice_id'])[tmp_list].agg('sum').reset_index()
    return curr_df

模型训练细节

• 通过使用 kfold 针对 train[‘linkid’].unique() 进行分折

train_user_id = train['linkid'].unique()
#output_preds = []
#feature_importance_df = pd.DataFrame()
#offline_score = []
#train['oof_pred_0'] = 0
#train['oof_pred_1'] = 0
#train['oof_pred_2'] = 0
for i, (train_idx, valid_idx) in enumerate(folds.split(train_user_id), start=1):
    train_X, train_y = train.loc[train['linkid'].isin(train_user_id[train_idx]), feats], train.loc[
        train['linkid'].isin(train_user_id[train_idx]), target]
    test_X, test_y = train.loc[train['linkid'].isin(train_user_id[valid_idx]), feats], train.loc[
        train['linkid'].isin(train_user_id[valid_idx]), target]

• 在每折训练过程中，使用scipy.optimize.minimize寻找权重，使得如下评价指标最优。

  • SciPy的scipy.optimize.minimize和Optuna都是用于优化问题的工具，‌虽然两者都用于优化问题，‌但scipy.optimize.minimize更侧重于函数优化，‌包括无约束和有约束的最优化问题，‌而Optuna则专注于超参数优化，‌特别是在机器学习模型的调优方面表现出色。‌

  • 应用场景：‌

    1. scipy.optimize.minimize主要用于函数优化，‌包括无约束和有约束的最优化问题。‌它适用于各种优化算法，‌如Nelder-Mead、‌BFGS、‌L-BFGS-B、‌SLSQP等，‌可以用于解决最小化（‌或最大化）‌目标函数的问题1。‌
    2. Optuna则更专注于超参数优化，‌特别是在机器学习模型的超参数调优方面表现出色。‌它通过模拟退火算法来寻找最优的超参数组合，‌适用于大规模的超参数搜索空间。‌

  • 优化方法：‌

    1. scipy.optimize.minimize提供了多种优化算法，‌可以根据问题的特性选择合适的算法进行优化。‌例如，‌对于非线性最小二乘问题，‌可以使用scipy.optimize.least_squares()方法。‌此外，‌它还支持有约束的优化问题，‌通过constraints参数来定义约束条件。‌
    2. Optuna则主要使用基于梯度的优化算法，‌通过模拟退火策略来探索超参数空间，‌寻找最优的超参数组合。‌它特别适合于具有大量超参数和复杂依赖关系的优化问题。‌

  • 使用目的：‌

    1. scipy.optimize.minimize的主要目的是最小化（‌或最大化）‌给定的目标函数，‌适用于各种数学和科学计算中的优化问题。‌
    2. Optuna旨在帮助机器学习从业者优化模型的超参数，‌以提高模型性能。‌它可以帮助确定哪些超参数对模型的整体性能有最显著的影响，‌并提供可视化的工具来查看搜索历史和超参数的重要性。‌

  • scipy.optimize.minimize一些主要参数及其含义：

    1. `fun`：需要最小化的目标函数。
    2. `x0`：初始值，是一个数值或数组，或者是一个具有多个初始值的列表。
    3. `method`：优化方法，可以是`Nelder-Mead`、`Powell`、`CG`、`BFGS`、`Newton-CG`等。
    4. `jac`：目标函数的梯度，如果提供了这个参数，那么`method`参数必须选择使用梯度的方法，如`BFGS`、`Newton-CG`等。
    5. `hessp`：目标函数的Hessian矩阵，如果提供了这个参数，那么`method`参数必须选择使用Hessian矩阵的方法，如`Newton-CG`。
    6. `bounds`：变量的边界，形式为`(min, max)`，可以是一个元组，也可以是一个包含多个元组的序列。
    7. `constraints`：约束条件，可以是`dict`类型或者`Constraint`对象。
    8. `tol`：停止准则的容差。
    9. `callback`：一个可调用的函数，它可以在每次迭代时被调用。
    10. `options`：优化选项，可以通过`scipy.optimize.OptimizeResult`类的`message`属性来查看。
    
    例如，一个简单的使用示例：
    
    ```python
    import numpy as np
    from scipy.optimize import minimize
    
    def func(x):
        return x[0] ** 2 + x[1] ** 2
    
    x0 = np.array([1, 1])
    res = minimize(func, x0, method='Nelder-Mead')
    print(res.x)
    

本方案优化目标：
score = report[‘0.0’][‘f1-score’] * 0.2 + report[‘1.0’][‘f1-score’] * 0.2 + report[‘2.0’][‘f1-score’] * 0.6
1）. 输入验证集的三类标签的预测结果、真实标签
2）. 优化的目标是：分别给三类标签的概率值以一个初始权重1，1，1；将三个权重分别乘以各自概率以优化上述评价指标

op = Optimizedkappa()
op.fit( train.loc[train['linkid'].isin(train_user_id[valid_idx]) 
			, ['oof_pred_0', 'oof_pred_1', 'oof_pred_2']].values,
        test_y )
        
class Optimizedkappa(object):
	def __init__(self):
	    self.coef_ = []
	
	def _kappa_loss(self, coef, X, y):
	    """
	    y_hat = argmax(coef*X, axis=-1)
	    :param coef: (1D array) weights
	    :param X: (2D array)logits
	    :param y: (1D array) label
	    :return: -f1
	    """
	    X_p = np.copy(X)
	    X_p = coef * X_p
	    report = classification_report(y, np.argmax(X_p, axis=1), digits=5, output_dict=True)
	    score = report['0.0']['f1-score'] * 0.2 + report['1.0']['f1-score'] * 0.2 + report['2.0']['f1-score'] * 0.6
	
	    return -score
	
	def fit(self, X, y):
	    loss_partial = partial(self._kappa_loss, X=X, y=y)
	    initial_coef = [1. for _ in range(len(set(y)))]
	    self.coef_ = sp.optimize.minimize(loss_partial, initial_coef, method='Powell')

模型预测

• 使用每折的模型对测试集进行预测并追加至lgb_preds列表中，使用np.mean对应位置求均值，输出最终三列标签的平均结果

    for i in range(0, 3):
        lgb_sub['lgb_pred_{}'.format(i)] = np.mean(lgb_preds, axis=0)[:, i]
    lgb_sub.to_csv('../user_data/lgb_sub_prob.csv', index=False)

NN方案

• 使用深度学习方法进行建模主要考虑三个方向：

  1. 如何利用实时/历史同期路况特征
     - 时序数据使用Lstm建模
  2. 如何利用路网拓扑关系
     - 使用图神经网络GNN，该方案使用GCN，图卷积神经网络
  3. 如何提取道路属性和日期属性多种类别特征组合
     - 使用推荐算法NFM ， 推荐算法擅长捕捉类别特征的高阶组合特征

• 完整的模型结构设计：
  

方案中的代码细节

• 类别特征

# 获取周形式的周期特征
test['day']=32
train['week_day']=train['day'].apply(lambda x:x%7)
test['week_day']=4

# train数据中current_slice_id会有负数情况，意义为前一天的时间片 ， 例如 0 、2 为当天的第一二个时间片 ， -1 则为昨天的最后时间片
# 但是统计发现future_slice_id最小值为0 ， 不存在else后面的可能
train['hour']=train['future_slice_id'].apply(lambda x:x//30 if x>=0 else (720+x)//30)
test['hour']=test['future_slice_id'].apply(lambda x:x//30 if x>=0 else (720+x)//30)

• 构建 全局的特征-Id的映射
  col_thre_dict 中每个key的值，本意设定的是次数的阈值，但是设置的太低了基本起不到过滤作用，可无视；
  后面所有涉及的操作都将用映射的Id进行

1. 这里主要输出是两个字典 mp_col_ids_indexs 以列名的形式存储每个特征下的特征->Id的映射关系，其中686110是给direction空值预留的编码，每个特征都进行了预留。
   
2. ids_indexs 按照col+‘_’+str(ids)、col+'_unknow’的形式存储所有的映射关系

col_thre_dict={'link':0.5, 'direction':0.5, 'path_class':0.5, 'speed_class':0.5
               , 'LaneNum':0.5, 'level':0.5, 'continuous_width':0.5
               ,'continuous_length':0.5,'continuous_speed_limit':0.5,'time_gap':0
               ,'current_slice_id':0.5,'future_slice_id':0.5,'week_day':0,'hour':0}
len(col_thre_dict)


mp_col_ids_indexs={} # 存储不同col的编码情况

ids_indexs={}    # 存储所有编码情况
ids_indexs['padding']=0
for col,thre in tqdm(col_thre_dict.items()):
    mp={}
    unknow=None
    #连续特征一个field
    if 'continuous_' in col:
        mp[col]=len(ids_indexs)
        ids_indexs[col]=len(ids_indexs)
        unknow=len(ids_indexs)
        ids_indexs[col+'_unknow']=len(ids_indexs)
        mp_col_ids_indexs[col]=[mp,unknow]
        continue
    if col=='link':
        curr_len=len(ids_indexs)
        ###use attr
        for i,ids in enumerate(attr_df['link'].values):
            ids_indexs[col+'_'+str(ids)]=i+curr_len
            mp[ids]=i+curr_len
        if thre!=0:
            unknow=len(ids_indexs)
            ids_indexs[col+'_unknow']=len(ids_indexs)
        mp_col_ids_indexs[col]=[mp,unknow]
        continue
        
    t=train[col].value_counts().reset_index()
    print(col+' or:',len(t))
    all_ids=t[t[col]>thre]['index'].values
    print(col+' new:',len(all_ids))
    
    print(col+' test or:',test[col].nunique())
    print(col+' test new:',test[test[col].isin(all_ids)][col].nunique())
    print('*'*50)
    
    curr_len=len(ids_indexs)
    for i,ids in enumerate(all_ids):
        ids_indexs[col+'_'+str(ids)]=i+curr_len
        mp[ids]=i+curr_len
    if thre!=0:
        unknow=len(ids_indexs)
        ids_indexs[col+'_unknow']=len(ids_indexs)
    mp_col_ids_indexs[col]=[mp,unknow]

• 子图拓扑获取
  这段代码非常精彩 ，因为整个路网拓扑涉及到的节点达60w，实际需要运用的较少；所以采用子图的思路进行局部图构建。
  有些依赖的小函数并未贴出。

import scipy.sparse as sp
link_toop_sub_graph={} #link_id:matrix
NUM_ADD_GRAPPH=4 #添加图的范围
MAX_NUMBER=200 #最大的子图数量
nums=[]
for link_id in tqdm(train_test_link):
    
    link_info=[] #记录边
    link_set=set([link_id]) # 记录相关点
    num=1 # 记录点的数量
    
    next_link=link_next_dict.get(link_id,[]) #下游
    before_link=link_before_dict.get(link_id,[]) #上游
    
    ###当前上下游link
    link_set,n=add_link_set(next_link,link_set)
    num+=n
    link_set,n=add_link_set(before_link,link_set)
    num+=n
    link_info.extend([[link_id,next_link_id] for next_link_id in next_link])
    link_info.extend([[link_id,before_link_id]for before_link_id in before_link])
    
    link_add_info=[]#其他边
    #扩展NUM_ADD_GRAPPH阶
    for graph_add in range(NUM_ADD_GRAPPH):
        if num==MAX_NUMBER:break
        next_link_next=[]
        before_link_before=[]
        #next link
        for sub_link_id in next_link:
            # 边关系  下游点数量   下游点列表
            sub_info,sub_num,sub_next_link= get_next_link(sub_link_id)
            # 遍历 sub_info 边关系 
            # 1.两个节点都已经在link_set 中 则追加 边关系至link_add_info
            # 2.两个节点都不在 则link_set追加两个点 同时也增加边关系
            # 3.某一个不存在 则link_set追加这个点 同时也增加边关系
            link_add_info,num,link_set=add_link_info(link_add_info,sub_info,num,MAX_NUMBER,link_set)
            # 将下游点追加至列表 ， 后续会将列表重新赋给next_link
            next_link_next.extend(sub_next_link)
            
            sub_info,sub_num,sub_next_link= get_before_link(sub_link_id)
            link_add_info,num,link_set=add_link_info(link_add_info,sub_info,num,MAX_NUMBER,link_set)
            next_link_next.extend(sub_next_link)
        #before link（bug：next_link_next：应该是before_link_before 但基本无影响）
        for sub_link_id in before_link:
            sub_info,sub_num,sub_before_link=get_before_link(sub_link_id)
            link_add_info,num,link_set=add_link_info(link_add_info,sub_info,num,MAX_NUMBER,link_set)
            next_link_next.extend(sub_before_link)
            
            sub_info,sub_num,sub_before_link= get_next_link(sub_link_id)
            link_add_info,num,link_set=add_link_info(link_add_info,sub_info,num,MAX_NUMBER,link_set)
            next_link_next.extend(sub_before_link)
        
        next_link=next_link_next
        before_link=before_link_before
    
    nums.append(num)
    link_info.extend(link_add_info)#所有边
    
    
    #转换成id
    edges=np.array(link_info) 
    """
    这里很关键 在每个子图中，当前link_id编号为0
    """
    link_mp={link_id:0} # 当前link_id编号为0  
    k=0
    for sub_link_id in set(edges.flatten()):
        if sub_link_id!=link_id:
            link_mp[sub_link_id]=k+1
            k+=1
    number=len(link_mp) #当前图 节点个数
    # 这里edges重新映射了边关系 ，且当前节点编号为0
    edges = np.array(list(map(link_mp.get, edges.flatten())),
                     dtype=np.int32).reshape(edges.shape)
    
    #转换成稀疏邻接矩阵
    """
    ( np.ones(edges.shape[0]), (edges[:, 0], edges[:, 1]) ) 
    np.ones(edges.shape[0]) :一维数组，包含矩阵中的非零元素,具体数值
    (edges[:, 0], edges[:, 1]) :(row ,col ) 一维数组，表示每个非零元素的行\列索引
    """
    # 这个将edges的第一列当作行，第二列当作列，构建边的连接关系 ,由于上面讲当前节点映射为0
    # 所以在这个子图当中 第一行代表着当前节点
    adj = sp.coo_matrix( ( np.ones(edges.shape[0]), (edges[:, 0], edges[:, 1]) ) 
                        ,shape=(number, number)
                        , dtype=np.float32
                       ).toarray()
    #获得现在的稀疏矩阵id
    link_ids=list(link_mp.keys())
    
    #paddding后转换成对称矩阵
    # 将adj数组在行和列方向上分别扩展至大小为MAX_NUMBER x MAX_NUMBER，扩展的部分用0填充。
    """
    (0, MAX_NUMBER-number)：表示在数组的第一个维度（行）上，前面填充0个元素，后面填充MAX_NUMBER-number个元素
    (0, MAX_NUMBER-number)：表示在数组的第二个维度（列）上，前面填充0个元素，后面填充MAX_NUMBER-number个元素
    mode='constant'：表示填充的模式为常数模式，即使用指定的常数值进行填充
    constant_values=(0)：表示填充的常数值为0。
    
    上面这种连接图 其实并不能是严格的无重复关系
    且本身是无向图的设定
    所以需要convert_symmetric转换成对称矩阵，函数里面增加了单位阵，并和转置矩阵相加，确保节点间的对称连接
    这样就会产生一个问题，假如某些节点间已经正常产生了连接，上述操作会增加数值由1->2等
    这个问题通过(adj>0).astype(np.int)
    """
    adj=np.pad(adj, ((0, MAX_NUMBER-number), (0, MAX_NUMBER-number)), mode='constant',constant_values=(0))
    adj=convert_symmetric(adj) #对称矩阵
    adj=(adj>0).astype(np.int)
    
    #D-1/2*A*D-1/2
    adj=normalize_adj(adj)
    #pad link_id
    link_ids=link_ids+[0]*(MAX_NUMBER-number)
    #save matrix
    link_toop_sub_graph[link_id]=(adj,link_ids)

上面代码 整体逻辑是 ：
1). 遍历训练集和测试集中所有的Link id , 通过下游找下游\上游；上游找下游\上游的方式 检索4阶范围内的不超过200的所有节点及边关系。

• 为什么不超过200？
  通过下图可以发现，经统计4阶范围内每个节点的关系节点99%都只有144个 ，200已经能较全面的覆盖，再大收益不高。
  

2). 在每个遍历节点的过程中 ， 重构节点图关系 ， 在该张图关系上使用GCN计算当前节点下图红线部分
GCN的计算逻辑：

• LSTM模型数据准备
  recent_feature 序列 即 近期的5个时间片的4特征数据
  history_feature序列 即 前28、21、14、7天的各5个时间片4特征数据的的堆叠

'recent_feature 序列'
recent_cols=[]
for i in range(1,6):
    recent_col=[ col for col in train.columns if 'recent_feature_{}'.format(i) in col]
    recent_cols.extend(recent_col)
train['recent_split_info']=[ i.reshape(4,5).reshape(5,4) for i in train[recent_cols].values] # 似乎不需要两次reshpe , 直接.reshape(5,4)
test['recent_split_info']=[ i.reshape(4,5).reshape(5,4) for i in test[recent_cols].values]

'history_feature 序列'
his_cols=[]
for i in range(1,6):
    his_col=[ col for col in train.columns if 'history_feature_cycle{}'.format(i) in col]
    his_cols.extend(his_col)
len(his_cols)
train['his_split_info']=[ i.reshape(4,20).reshape(20,4) for i in train[his_cols].values]
test['his_split_info']=[ i.reshape(4,20).reshape(20,4) for i in test[his_cols].values]

• 整个Model搭建

class DiDi_Model(nn.Module):
    def __init__(self,embedding_num,embedding_dim,field_dims=None):
        '''
        field_dims：the number of fileds
        embedding_num : sum of the index of all fields
        embedding_dim : the dim of embedding
        '''
        super(DiDi_Model, self).__init__()
        self.model_name = 'zy_Model'
        self.field_dims=field_dims  # 并未实际使用
        self.embedding_num=embedding_num  # len(ids_indexs) ,前文所说 ids_indexs编码了所有特征
        self.embedding_dim=embedding_dim  # 想要的维度 例如 32
        
        
        #link原始特征embedding
        # link_embeddding_matrix_cate 这个是link的道路属性原有数值
        self.link_or_em=nn.Embedding(link_embeddding_matrix_cate.shape[0],link_embeddding_matrix_cate.shape[1])
        self.link_or_em.weight.data.copy_(torch.from_numpy(link_embeddding_matrix_cate))
        self.link_or_em.requires_grad = False
        
        self.width_em=nn.Sequential(nn.Linear(1,self.embedding_dim))
        self.length_em=nn.Sequential(nn.Linear(1,self.embedding_dim))
        self.speed_em=nn.Sequential(nn.Linear(1,self.embedding_dim))
        
        #FM的一阶
        self.first_em=nn.Embedding(self.embedding_num,1)
        #FM的二阶
        self.embdedding_seq=nn.Embedding(self.embedding_num,embedding_dim)
        self.bi = Bi_interaction() # NFM的核心计算逻辑 , 实现高阶的特征交叉衍生
        
        #lstm
        # 这里定义两个lstm 分别用于rencet和history
        input_dim=4
        output_put_dim=8
        self.lstm_seq=nn.LSTM(input_dim,output_put_dim,2,batch_first=True,bidirectional=False)
        self.lstm_seq_1=nn.LSTM(input_dim,output_put_dim,2,batch_first=True,bidirectional=False)
        
        #GAE
        self.gcn=GCN(K=2,feat_dim=9*32)
        

        self.embed_output_dim = embedding_dim+output_put_dim*2+output_put_dim*2*4+9*32

        self.mlp=nn.Sequential(
                               nn.Dropout(0.5),
                               nn.Linear(self.embed_output_dim,self.embed_output_dim//2),
                               nn.BatchNorm1d(self.embed_output_dim//2),
                               nn.ReLU(True),
                               nn.Dropout(0.3),
                               nn.Linear(self.embed_output_dim//2,3))
     # mean pooling   
    def mask_mean(self,x,mask=None):
        if mask!=None:
            mask_x=x*(mask.unsqueeze(-1))
            x_sum=torch.sum(mask_x,dim=1)
            re_x=torch.div(x_sum,torch.sum(mask,dim=1).unsqueeze(-1))
        else:
            x_sum=torch.sum(x,dim=1)
            re_x=torch.div(x_sum,x.size()[1])
        return re_x
    
    # max pooling
    def mask_max(self,x,mask=None):
        if mask!=None:
            mask=mask.unsqueeze(-1)
            mask_x=x-(1-mask)*1e10
            x_max=torch.max(mask_x,dim=1)
        else:
            x_max=torch.max(x,dim=1)
        return x_max[0]

    def forward(self,category_index,category_value
    ,recent_split_info,his_split_info,link_seq,link_graph,label,is_test=False):
        
        batch_size=category_index.size()[0]
        
        ##FM二阶
        # 这里将连续特征也进行了向量嵌入
        # category_value在上文的处理过程中全设置为 , 此处就无作用了
        seq_em=self.embdedding_seq(category_index)*category_value.unsqueeze(2)
        # 进行高阶衍生
        x2=self.bi(seq_em)
        
        #lstm
        # 近期特征
        f,_=self.lstm_seq_1(recent_split_info)
        fmax=self.mask_max(f)
        fmean=self.mask_mean(f)
        recent_features=torch.cat([fmean,fmax],dim=1)
        
        # 历史特征
        # 历史特征共用一个lstm模型
        his_features=[]
        for i in range(4):
            his_split_info_sample=his_split_info[:,i*5:(i+1)*5,:]
            f,_=self.lstm_seq(his_split_info_sample)
            fmax=self.mask_max(f)
            fmean=self.mask_mean(f)
            his_features.append(fmax)
            his_features.append(fmean)
        his_features=torch.cat(his_features,dim=-1)
        
        
        ###GCN
        #######node feat
        number_graph_node=link_seq.size(1)
        # embdedding_seq 是所有id的嵌入向量表
        link_feat_link_id=self.embdedding_seq(link_seq)
        
        #  每个节点原始那些特征 ，长度宽度啥的
        link_feat=self.link_or_em(link_seq)# B*number_graph_node*8
        
        link_feat_cate=link_feat[:,:,:5].long() #类别特征  B*number_graph_node*5
        link_feat_cate=self.embdedding_seq(link_feat_cate).view(batch_size,number_graph_node,-1)
        
        link_feat_width=self.width_em(link_feat[:,:,5].float().unsqueeze(2)) # 连续特征使用width_em进行计算映射
        link_feat_length=self.length_em(link_feat[:,:,6].float().unsqueeze(2)) 
        link_feat_speed=self.speed_em(link_feat[:,:,7].float().unsqueeze(2)) 
        # dim = -1 沿着最后一个维度拼接
        # 相当于 列 形式的特征
        link_feat=torch.cat([link_feat_cate,link_feat_width,link_feat_length,link_feat_speed,link_feat_link_id],dim=-1) # B*node*(dim*9)
        
        
        """
        link_graph size 为 bs,node,node
        link_feat  size 为 bs,node,(dim*9)
        gcn_out    size 为 bs,node,(dim*9)
        之所以node取0,是因为在前面对每张子图进行编码时，当前节点的编码为0,即第一行的向量代表着当前节点
        """
        gcn_out=self.gcn(link_graph,link_feat)
        gcn_out=gcn_out[:,0,:].squeeze()
        

        
        #DNN全连接
        x2=torch.cat([x2,recent_features,his_features,gcn_out],dim=1)
        x3=self.mlp(x2)
        
        out=x3
        if not is_test:
        	# 类别加权
            loss_fun=nn.CrossEntropyLoss(torch.tensor([0.1,0.3,0.6]).to(DEVICE))
            loss=loss_fun(out,label)
            return loss,F.softmax(out,dim=1)
        else:
            loss_fun=nn.CrossEntropyLoss()
            loss=loss_fun(out,label)
            return loss,F.softmax(out,dim=1)
net = DiDi_Model(field_dims=len(col_thre_dict),embedding_num=len(ids_indexs),embedding_dim=32)
print('# Model parameters:', sum(param.numel() for param in net.parameters()))

• 训练细节
  通过validation_fn (里面实现了上文自定评价指标的计算方式 )，来控制最佳模型的保存

-- validation_fn定义代码
def metric_fn(preds,real_labels):
    ##for log
    preds=np.argmax(preds,axis=1)
    f1=f1_score(real_labels, preds,average=None)
    print(f1)
    return 0.2*f1[0]+0.2*f1[1]+0.6*f1[2]

def validation_fn(model,val_loader,is_test=False):
    model.eval()
    bar = tqdm(val_loader)
    preds=[]
    labels=[]
    weights=[]
    loss_all=[]
    for i,feat in enumerate(bar):
        category_index,category_value,recent_split_info,his_split_info,link_seq,link_graph,label=(_.to(DEVICE) for _ in feat)
        loss,p= model(category_index,category_value,recent_split_info,his_split_info,link_seq,link_graph,label,is_test=True)
        preds.append(p.detach().cpu().numpy())
        labels.append(label.detach().cpu().numpy())
        loss_all.append(loss.item())
    preds=np.concatenate(preds)
    labels=np.concatenate(labels)
    if not is_test:
        score=metric_fn(preds.squeeze(),labels)
        return np.mean(loss_all),score
    else:
        return preds.squeeze()

-- 训练代码
 val_loss,mse=validation_fn(model,val_loader)
# losses.append( 'train_loss:%.5f, score: %.5f, best score: %.5f\n val loss: %.5f\n' %
#     (np.array(loss_ages).mean(),mse,best_vmetric,val_loss))
# print(losses[-1])
 if mse>=best_vmetric:
     torch.save(model.state_dict(),model_save_path)
     best_vmetric=mse
     no_improve_epochs=0
     print('improve save model!!!')
 else:
     no_improve_epochs+=1
 if no_improve_epochs==early_stop_epochs:
     print('no improve score !!! stop train !!!')
     break

模型融合

融合策略是每个类别的概率计算都是 lgb0.3+nn0.7

test_sub['lgb_pred_0']=test_sub['lgb_pred_0']*0.3+nn_prob[:,0]*0.7
test_sub['lgb_pred_1']=test_sub['lgb_pred_1']*0.3+nn_prob[:,1]*0.7
test_sub['lgb_pred_2']=test_sub['lgb_pred_2']*0.3+nn_prob[:,2]*0.7
test_sub['label']=np.argmax(test_sub[['lgb_pred_0','lgb_pred_1','lgb_pred_2']].values,axis=1)+1