[校企合作]oneAPI的使用之预测淡水质量

一：问题简介

1.问题描述

  淡水是我们最重要和最稀缺的自然资源之一，仅占地球总水量的 3%。它几乎触及我们日常生活的方方面面，从饮用、游泳和沐浴到生产食物、电力和我们每天使用的产品。获得安全卫生的供水不仅对人类生活至关重要，而且对正在遭受干旱、污染和气温升高影响的周边生态系统的生存也至关重要。

2.预期解决方案

  通过参考英特尔的类似实现方案，预测淡水是否可以安全饮用和被依赖淡水的生态系统所使用，从而可以帮助全球水安全和环境可持续性发展。这里分类准确度和推理时间将作为评分的主要依据。

3.数据集来源

https://filerepo.idzcn.com/hack2023/datasetab75fb3.zip

二：数据处理

1.数据预处理

导所需要的包，若没有则用pip install下载。

import modin.pandas as pd
import os
import xgboost
from xgboost import XGBClassifier
import time
import warnings

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.io as pio
import plotly.graph_objects as go
from sklearn.utils import resample

from sklearn.model_selection import train_test_split, StratifiedKFold, GridSearchCV, RandomizedSearchCV
from sklearn.preprocessing import RobustScaler
from sklearn.metrics import roc_auc_score, roc_curve, auc, accuracy_score, f1_score
from sklearn.preprocessing import StandardScaler
import sklearn
from sklearn.metrics import precision_recall_curve, average_precision_score

为了方便查出输出结果，对运行环境进行设置,使其支持中文显示，并配置 pandas 最大列数为 100,避免其对列数进行隐藏,这样利于后续查看各样本列的情况,代码如下:

os.environ['NLSLANG']='SIMPLIFIED CHINESE_CHINA.UTF8'
pd.set_option( 'display.max_columns', 100)

并导入我们此次intel加速组件所需要的包，来加速训练和预测速度。

import daal4py as d4p
os.environ["MODIN_ENGINE"]="dask"
import modin.pandas as pd
from modin.config import Engine
Engine.put("dask")

from sklearnex import patch_sklearn

patch_sklearn()

2 数据探索

我们首先数据集导入，并且对数据规模和数据集情况进行查看。

import pandas
data=pandas.read_csv('./dataset.csv')
print('数据规模:{}\n'.format(data.shape))
display(data.head())
data.info()

结果如下：

通过数据的结果我们可以看到，我们该数据集的规模为有5956842行，24列，对应每列的列名和属性也可以看到，其中float64类型占 19 列,int64类型字段占 2 列,object类型占 3 列，总共占用内存空间大约 1.1GB。

将Color， Source，Month等英文字段进行因子化处理为数字型变量，代码如下：

display(data.head())
factor = pd.factorize(data['Color'])
print(factor)
data.Color = factor[0]
factor = pd.factorize(data['Source'])
print(factor)
data.Source = factor[0]
factor = pd.factorize(data['Month'])
print(factor)
data.Month = factor[0]

data

 输出的结果如下：

3.特征选择 

我们首先对数据集中的列也就是属性进行相关性分析，并输出对应的柱状图，代码如下：

# 相关性分析
bar = data.corr()['Target'].abs().sort_values(ascending=False)[1:]

plt.bar(bar.index, bar, width=0.5)
# 设置figsize的大小
pos_list = np.arange(len(data.columns))
params = {
    'figure.figsize': '20, 10'
}
plt.rcParams.update(params)
plt.xticks(bar.index, bar.index, rotation=-60, fontsize=10)
plt.show()

结果如下：

通过柱状图我们可以看到最后面的几列是不相关的列，所以我们就可以直接将其删除，代码如下：

 # 删除不相关的列
 
 data = data.drop(
    columns=['Index', 'Day', 'Time of Day', 'Month', 'Water Temperature', 'Source', 'Conductivity', 'Air Temperature'])

之后我们就对数据集中的缺失值和重复值进行查看，代码如下：

display(data.isna().sum())
missing = data.isna().sum().sum()
duplicates = data.duplicated().sum()
print("\n数据集中 {:,.0f} 缺失值.".format(missing))
print("数据集中 {:,.0f} 重复值.".format(duplicates))

输出的结果如下：

可以看到我们的数据中有1603858条缺失值和143032条重复值，占比很大，所以我们要对缺失值和重复值进行处理，代码如下：

data = data.fillna(0)
data = data.drop_duplicates()

from scipy.stats import pearsonr

variables = data.columns
data = data

var = data.var()
numeric = data.columns
df = data.fillna(0)
for i in range(0, len(var) - 1):
    if var[i] <= 0.1:  # 方差大于10％
        print(variables[i])
        data = data.drop(numeric[i], 1)
variables = data.columns

for i in range(0, len(variables)):
    x = data[variables[i]]
    y = data[variables[-1]]
    if pearsonr(x, y)[1] > 0.05:
        print(variables[i])
        data = data.drop(variables[i], 1)

variables = data.columns
print(variables)
print(len(variables))

处理之后我们在此对缺失值和重复值进行查看，即查看处理结果，代码如下：

display(data.isna().sum())
missing = data.isna().sum().sum()
duplicates = data.duplicated().sum()
print("\n数据集中 {:,.0f} 缺失值.".format(missing))
print("数据集中 {:,.0f} 重复值.".format(duplicates))

输出的结果如下：

可以看到处理过后，我们的数据集中已经没有缺失值和重复值了。

之后我们就对数据进行数据的不平衡处理，我们首先查看一下数据的平衡情况，这里是输出了target的具体分布，并用饼图进行了直观的查看，代码如下：

print(data.Target.value_counts())
target = data.Target.value_counts()
target.rename(index={1: 'state 1', 0: 'state o'}, inplace=True)
plt.pie(target, [0, 0.05], target.index, autopct='%1.1f%%')
plt.show()

输出的结果如下：

 

 

我们这里使用下采样对数据进行平衡处理。

  首先，将数据框df中的特征和目标变量分别赋值给X和y。然后，使用imblearn库中的RandomUnderSampler()函数进行下采样，使得正负样本比例接近1:1。这里使用了随机数种子（random_state=21）来确保每次采样的结果相同。接下来，使用train_test_split()函数将X和y分割成训练集和测试集，其中测试集占总数据的20%。同时，使用StandardScaler()函数对训练集和测试集进行标准化处理，使得每个特征的均值为0，标准差为1，以提高模型的性能。最后，输出训练集和测试集的形状。

代码如下：

from imblearn.under_sampling import RandomUnderSampler
import datetime

X = data.iloc[:, 0:len(df.columns.tolist()) - 1].values
y = data.iloc[:, len(df.columns.tolist()) - 1].values

# # 下采样
under_sampler = RandomUnderSampler(random_state=21)
X, y = under_sampler.fit_resample(X, y)

X = data.drop('Target', axis=1)
y = data['Target']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=21)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
print("Train Shape: {}".format(X_train_scaled.shape))
print("Test Shape: {}".format(X_test_scaled.shape))

X_train, X_test = X_train_scaled, X_test_scaled

输出的结果如下：

三：模型使用

1.模型定义

模型采用XGBoost，并使用随机网格搜索（RandomizedSearchCV）进行优化，模型参数定义如下

1.树的最大深度(max_depth),其取值范围是10，15，20;
2.评估器数量参数(n _estimators),可选取值为 10、50、80、100;
3.gamma参数 （min_split_loss），可选值为0，1，2， 默认是0，分裂节点时，损失函数减小值只有大于等于gamma节点才分裂，gamma值越大，算法越保守，越不容易过拟合，但性能就不一定能保证，需要平衡。
4.colsample_bytree：可选值为0.3, 0.5, 1，默认= 1，列采样率，也就是特征采样率。范围为（0，1]
5.subsample：可选值为0.9, 1，默认值= 1，构建每棵树对样本的采样率，如果设置成0.5，XGBoost会随机选择一半的样本作为训练集。范围：（0,1]
6.min_child_weight：可选值为4, 6, 8，默认值= 1，如果新分裂的节点的样本权重和小于min_child_weight则停止分裂 。这个可以用来减少过拟合，但是也不能太高，会导致欠拟合。范围：[0，∞]
7.alpha（reg_alpha）：可选值为3, 4, 5， 默认= 0，权重的L1正则化项。增加此值将使模型更加保守。代码如下：

from sklearn.metrics import make_scorer, precision_score, recall_score, accuracy_score,f1_score, roc_auc_score
param_grid = {
    'max_depth': [10, 15, 20],
    "gamma": [0, 1, 2], # -> 0
    "subsample": [0.9, 1], # -> 1
    "colsample_bytree": [0.3, 0.5, 1], # -> 1
    'min_child_weight': [4, 6, 8], # -> 6
    "n_estimators": [10,50, 80, 100], # -> 80
    "alpha": [3, 4, 5] # -> 4
}

scorers = {
    'precision_score': make_scorer(precision_score),
    'recall_score': make_scorer(recall_score),
    'accuracy_score': make_scorer(accuracy_score),
    'f1_score': make_scorer(f1_score),
    'roc_auc_score': make_scorer(roc_auc_score),
}

之后对定义XGBoost分类器，代码如下：

xgb = XGBClassifier(
    learning_rate=0.1,
    n_estimators=15,
    max_depth=12,
    min_child_weight=6,
    gamma=0,
    subsample=1,
    colsample_bytree=1,
    objective='binary:logistic', # 二元分类的逻辑回归，输出概率
    nthread=4,
    alpha=4,
    scale_pos_weight=1,
    seed=27)

2模型优化并训练 

我们首先使用RandomizedSearchCV从scikit-learn库来进行超参数优化的，特别是针对XGBoost分类器。

refit_score = "f1_score"
 
start_time = datetime.datetime.now()
print(start_time)
rd_search = RandomizedSearchCV(xgb, param_grid, n_iter=10, cv=3, refit=refit_score, scoring=scorers, verbose=10, return_train_score=True)
rd_search.fit(X_train, y_train)
print(rd_search.best_params_)
print(rd_search.best_score_)
print(rd_search.best_estimator_)
print(datetime.datetime.now() - start_time)

参数定义如下： 

• xgb: 要优化的模型，这里假设是一个XGBoost分类器实例。
• param_grid: 超参数的搜索范围
• n_iter=10：随机采样的参数组合数量。
• cv=3: 使用3折交叉验证。
• refit=refit_score: 使用F1得分来重新拟合模型。
• scoring=scorers: 定义评分函数或评分函数的列表。
• verbose=10: 输出详细的搜索信息。
• return_train_score=True: 在返回结果中包括训练集上的得分。

最后输出的结果如下：

 

3.模型使用

我们现在已经有经过超参数优化后的模型了，之后我们就将模型使用在测试集上，并输出相应的混淆矩阵。代码如下：

from dask_ml.xgboost import XGBClassifier
from sklearn.metrics import confusion_matrix
y_pred = rd_search.best_estimator_.predict(X_test)

# confusion matrix on the test data.
print('\nConfusion matrix of Random Forest optimized for {} on the test data:'.format(refit_score))
print(pd.DataFrame(confusion_matrix(y_test, y_pred),
                   columns=['pred_neg', 'pred_pos'], index=['neg', 'pos']))

输出结果如下：

 

四：结果输出

 

from datetime import datetime
# 记录开始时间
inference_start_time = datetime.now()
 
# 模型推理代码
y_pred = rd_search.best_estimator_.predict(X_test)
 
# 计算模型推理时间
inference_time = datetime.now() - inference_start_time
print("模型推理时间:", inference_time)

print(accuracy_score(y_test, y_pred))
print(recall_score(y_test, y_pred))

print(roc_auc_score(y_test, rd_search.best_estimator_.predict_proba(X_test)[:, 1]))
print(f1_score(y_test, y_pred))

 运行的时间为：02.114928秒

准确率、查全率、AUC 值、F1值分别如下： 

五：英特尔组件的使用 

在本项目中主要使用了以下的组件

1.daal4py组件
可将训练好的xgboost转换为 daal4py 模型，以便进一步改进预测时间性能, 利用底层的英特尔® 高级矢量扩展指令集（英特尔® AVX-512）硬件，最大限度地提高英特尔® 至强® 处理器上的梯度提升性能。

2.modin
Modin 使用Ray、Dask或Unidist提供一种轻松的方式来加速您pandas 、脚本和库。与其他分布式 DataFrame 库不同，Modin 提供与现有 pandas 代码的无缝集成和兼容性。即使使用 DataFrame 构造函数也是相同的。

3.sklearnex
借助面向 Scikit-learn* 的英特尔® 扩展，加速 Scikit-learn ，并且仍然完全符合所有 Scikit-Learn API 和算法。英特尔® Extension for Scikit-learn* 是一款免费软件 AI 加速器，可为各种应用程序带来超过 10-100 倍的加速。

六：项目总结

 本次的项目相较于上次的项目来说数据量更大、并且才用到的英特尔组件也与之前大不相同，让我更进一步的了解和学习到了英特尔组件和机器学习相关的知识。也希望后续我可以更加努力的学习并了解各种模型，并在后续的项目中利用到更多的oneAPI组件来加速。