{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# 新洋教育Kaggle零基础教学计划 - 数据挖掘项目\n# 预测建筑物的能源消耗\n\nASHRAE（American Society of Heating, Refrigerating and Air-Conditioning Engineers），中文名称“美国采暖、制冷与空调工程师学会”，于1894年在美国纽约成立，是由暖通空调（HVAC）工程师所组成的学会，全球拥有超过54,000名成员。协会及其成员专注于建筑系统、能源效率、室内空气质量、制冷和行业内的可持续性。通过调研、标准编写、出版和继续教育，ASHRAE发展至现在的规模。\n\n\n\n问：夏天给大楼降温需要多少钱？\n\n答：非常多！政府正在进行投资，以降低能源成本，减少排放。但是问题是，这些改进是否真的有效？\n\n在这次竞赛中，我们通过**预测冷水表、电表、热水表和蒸汽表的读数**来对这些节能投资进行更好的估计。数据来自近三年来1000栋建筑中的各表读数。大型投资者和金融机构将更倾向于在这一领域投资，以提高建筑能源使用效率。\n\n>**提示：**Code 和 Markdown 区域可通过 **Shift + Enter** 快捷键运行。此外，Markdown可以通过双击进入编辑模式。\n\n我们将这个notebook分为不同的步骤，你可以使用下面的链接来浏览此notebook。\n\n* [Step 1](#step1): 导入数据\n* [Step 2](#step2): 探索性数据分析\n* [Step 3](#step3): 数据预处理\n* [Step 4](#step4): CatBoost\n* [Step 5](#step5): 结果预测\n\n在该项目中包含了如下的问题：\n\n* [问题 1](#question1): 回顾课上内容并查阅资料，归纳总结数据预处理需要的步骤。\n* [问题 2](#question2): 思考此处为何要进行对数转换。\n* [问题 3](#question3): 查阅资料，总结LightGBM与CatBoost的差异。"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"#导入必要的库\nimport pandas as pd\nimport numpy as np\nimport os\nimport gc\nimport copy\nimport warnings\n\nimport catboost as cb\nfrom catboost import CatBoostRegressor\nfrom sklearn.metrics import mean_squared_log_error\nfrom sklearn.model_selection import StratifiedKFold, KFold\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nwarnings.filterwarnings('ignore')\npd.set_option(\"max_columns\", 500)\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# 1. 导入数据\n\nASHRAE给出的数据包含大量的特征数据，包括仪表读数，天气和建筑的数据等等。该问题为典型的监督学习问题。比赛举办方提供了6个csv文件，包括5个数据集与1个提交样本。其中数据集的字段含义为：\n\n`[train/test].csv`\n- building_id：建筑原数据的外键\n- meter : 仪表的id码, {0: 电表 , 1: 冷水表, 2: 蒸汽表, 3: 热水表}，不是每栋建筑都有全部类型的仪表\n- timestamp：读表的时间\n- meter_reading：目标变量, 用千瓦时（或等效值）表示的能耗。这是带有测量误差的真实数据，其中site0的电表读数出现问题，单位是千英热\n\n`building_meta.csv`\n- site_id: 天气文件的外键\n- building_id: training.csv对应的外键\n- primary_use: 基于EnergyStar property type definitions的建筑物活动的主要类别的指标（education, office…)\n- square_feet: 建筑物的总建筑面积\n- year_built: 建筑完成的年份\n- floor_count: 建筑物层数\n\n`weather_[train/test].csv`：气象站提供的气象数据,尽可能接近现场。\n- site_id: 天气文件的外键\n- air_temperature: 气温，单位为摄氏度\n- cloud_coverage: 天空中被云层覆盖的部分，单位为oktas\n- dew_temperature: 露点温度，单位为摄氏度\n- precip_depth_1_hr: 降水深度，单位为毫米\n- sea_level_pressure: 海平面压力，单位为毫巴/公顷\n- wind_direction: 风向，使用的是指南针方向（0-360）\n- wind_speed: 风速，单位为米每秒"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/ashrae-energy-prediction/train.csv\", parse_dates=[\"timestamp\"])\ntest = pd.read_csv(\"../input/ashrae-energy-prediction/test.csv\", parse_dates=[\"timestamp\"])\nbuilding = pd.read_csv('../input/ashrae-energy-prediction/building_metadata.csv')\nweather_train = pd.read_csv('../input/ashrae-energy-prediction/weather_train.csv', parse_dates=[\"timestamp\"])\nweather_test = pd.read_csv(\"../input/ashrae-energy-prediction/weather_test.csv\", parse_dates=[\"timestamp\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# 2. 探索性数据分析\n\n- 根据主办方提供的数据，\"meter\"代表仪表的类型，对应关系为 {0: 电表 , 1: 冷水表, 2: 蒸汽表, 3: 热水表}\n- 观察各个建筑的仪表读数，可以发现，其中一些建筑的读数在某些区间出现了持续为0的异常情况,也有读数异常的高的值\n- 可以通过修改`site`,`meter_type`,`primary_use`三个参数选择想要绘制的数据\n\n关联train表格和building表格用于作图"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_plot = train.merge(building, on='building_id', how='left')\n\nsite = 0 #建筑物的地点\nmeter_type = 1 #仪表的类型\nprimary_use = 'Education' #建筑物的用途\n\nr = int(np.ceil(len(train_plot[(train_plot['site_id'] == site) & (train_plot['primary_use'] == primary_use) & (train_plot['meter'] == meter_type)]['building_id'].value_counts(dropna=False).index.to_list())/2))\nfig, axes = plt.subplots(r,2,figsize=(14, 36), dpi=100)\nfor i, building_id in enumerate(train_plot[(train_plot['site_id'] == site) & (train_plot['primary_use'] == primary_use) & (train_plot['meter'] == meter_type)]['building_id'].value_counts(dropna=False).index.to_list()):\n train_plot[(train_plot['site_id'] == site) & (train_plot['primary_use'] == primary_use) & (train_plot['meter'] == meter_type) & \n (train_plot['building_id'] == building_id )][['timestamp', 'meter_reading']].set_index('timestamp').resample('H').mean()['meter_reading'].plot(ax=axes[i%r][i//r], \n alpha=0.8, label='By hour', color='tab:blue').set_ylabel('Mean meter reading', fontsize=13);\n train_plot[(train_plot['site_id'] == site) & (train_plot['primary_use'] == primary_use) & (train_plot['meter'] == meter_type) & \n (train_plot['building_id'] == building_id )][['timestamp', 'meter_reading']].set_index('timestamp').resample('D').mean()['meter_reading'].plot(ax=axes[i%r][i//r], \n alpha=1, label='By day', color='tab:orange').set_xlabel('');\n axes[i%r][i//r].legend();\n axes[i%r][i//r].set_title('building_id: ' + str(building_id ), fontsize=13);\n plt.subplots_adjust(hspace=0.45)\n \ndel train_plot,fig,axes,r\ngc.collect();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# 3. 数据预处理\n\n\n### __问题 1:__\n\n回顾课上内容并查阅资料，归纳总结数据预处理需要的步骤。\n\n__回答:__ "},{"metadata":{},"cell_type":"markdown","source":"## 3.1 数据类型转换"},{"metadata":{"trusted":true},"cell_type":"code","source":"def compress_dataframe(df):\n '''将所有数据的类型都转换为数值型'''\n result = df.copy()\n for col in result.columns:\n col_data = result[col]\n dn = col_data.dtype.name\n if dn == \"object\":\n result[col] = pd.to_numeric(col_data.astype(\"category\").cat.codes, downcast=\"integer\")\n elif dn == \"bool\":\n result[col] = col_data.astype(\"int8\")\n elif dn.startswith(\"int\") or (col_data.round() == col_data).all():\n result[col] = pd.to_numeric(col_data, downcast=\"integer\")\n else:\n result[col] = pd.to_numeric(col_data, downcast='float')\n return result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.2 缺失值填充与特征扩展"},{"metadata":{"trusted":true},"cell_type":"code","source":"def set_time(df):\n df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n #这里将timestamp转换成了16年1月1日0点开始计算的小时数‘//’代表除法运算后取整\n return df\n\n# 根据分析得出各个site来自哪个时区，来修正时间\n# https://www.kaggle.com/patrick0302/locate-cities-according-weather-temperature\nsite_GMT_offsets = [-5, 0, -7, -5, -8, 0, -5, -5, -5, -6, -7, -5, 0, -6, -5, -5]\n\n#转换天气数据表格中的时间,并填充缺失值\ndef weather_set_time(df,time_zone):\n df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n \n GMT_offset_map = {site: offset for site, offset in enumerate(site_GMT_offsets)}\n df.timestamp = df.timestamp + df.site_id.map(GMT_offset_map)\n #根据时区的不同，统一时间\n site_dfs = []\n for site_id in df.site_id.unique():\n # 确保包括所有可能的小时数\n site_df = df[df.site_id == site_id].set_index(\"timestamp\").reindex(time_zone)\n site_df.site_id = site_id\n for col in [c for c in site_df.columns if c != \"site_id\"]:\n site_df[f\"had_{col}\"] = ~site_df[col].isna()\n site_df[col] = site_df[col].interpolate(limit_direction='both', method='linear')\n # 这里使用中位数来填充缺失值\n site_df[col] = site_df[col].fillna(df[col].median())\n site_dfs.append(site_df)\n df = pd.concat(site_dfs).reset_index() # make timestamp back into a regular column\n for col in df.columns:\n if df[col].isna().any(): df[f\"had_{col}\"] = ~df[col].isna()\n #如果某列其中有缺失值，就增加一列新的特征：had_xxx 表示这一行在xxx这一列是否有记录\n return df\n \n#增加星期，月份，时间的特征\ndef _add_time_features(X):\n return X.assign(tm_day_of_week=((X.timestamp // 24) % 7), tm_hour_of_day=(X.timestamp % 24))\n\nbuilding = compress_dataframe(building.fillna(-1)).set_index(\"building_id\")\n\ntrain = compress_dataframe(set_time(train))\ntest = compress_dataframe(set_time(test)).set_index(\"row_id\")\nweather_train = compress_dataframe(weather_set_time(weather_train,range(8784))).set_index([\"site_id\", \"timestamp\"])\nweather_test = compress_dataframe(weather_set_time(weather_test,range(8784,26304))).set_index([\"site_id\", \"timestamp\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.3 关联数据"},{"metadata":{"trusted":true},"cell_type":"code","source":"def combined_data(df,weather):\n df = compress_dataframe(df.join(building, on=\"building_id\").join(weather,\n on=[\"site_id\", \"timestamp\"]).fillna(-1))\n return df.drop(columns=[\"meter_reading\"]),df.meter_reading","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.4 异常值处理 "},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_is_bad_zero(Xy_subset, min_interval=48, summer_start=3000, summer_end=7500):\n #夏天，3000/24=125，7500/24=312.5,第125天到第312.5天为夏天。\n\n meter = Xy_subset.meter_id.iloc[0]\n is_zero = Xy_subset.meter_reading == 0 #返回读数为0的电表的indices\n if meter == 0:\n #电表的度数不应该为0，所以电表（meter为0）读数为0的行从training dataframe中drop掉\n return is_zero\n\n transitions = (is_zero != is_zero.shift(1))#出现0和非0变化的位置\n all_sequence_ids = transitions.cumsum()#到各位置出现的变化的和，是一个pd.Seires\n ids = all_sequence_ids[is_zero].rename(\"ids\")#将其中读数为0的提取出来\n if meter in [2, 3]:\n # 蒸汽和热水有可能在夏天被关闭\n keep = set(ids[(Xy_subset.timestamp < summer_start) |\n (Xy_subset.timestamp > summer_end)].unique())#不在夏天的indices\n is_bad = ids.isin(keep) & (ids.map(ids.value_counts()) >= min_interval) \n #将不在夏天却被关闭的蒸汽和热水表提取出来，至少被关闭了48小时以上的\n elif meter == 1:\n time_ids = ids.to_frame().join(Xy_subset.timestamp).set_index(\"timestamp\").ids#将ids和timestamp对应起来\n is_bad = ids.map(ids.value_counts()) >= min_interval#关闭时间大于48小时的\n\n # 冷水在冬天可能被关闭\n jan_id = time_ids.get(0, False)#一月份的开始的id\n dec_id = time_ids.get(8283, False)#十二月份开始的id\n if (jan_id and dec_id and jan_id == time_ids.get(500, False) and\n dec_id == time_ids.get(8783, False)):\n #如果一月500小时和十二月500小时的读表都为0的话\n is_bad = is_bad & (~(ids.isin(set([jan_id, dec_id]))))\n #将这一部分的的行从is_bad中删除\n else:\n raise Exception(f\"Unexpected meter type: {meter}\")\n\n result = is_zero.copy()\n result.update(is_bad)\n return result\n\ndef find_bad_zeros(X, y):\n \"\"\"返回仅包含应该删除的行的Index\"\"\"\n Xy = X.assign(meter_reading=y, meter_id=X.meter)\n is_bad_zero = Xy.groupby([\"building_id\", \"meter\"]).apply(make_is_bad_zero)\n return is_bad_zero[is_bad_zero].index.droplevel([0, 1])\n\ndef find_bad_sitezero(X):\n \"\"\"返回Site 0 读数异常的行的index.\"\"\"\n return X[(X.timestamp < 3378) & (X.site_id == 0) & (X.meter == 0)].index\n\ndef find_bad_building1099(X, y):\n \"\"\"返回建筑1099的读数异常高的行的index .\"\"\"\n return X[(X.building_id == 1099) & (X.meter == 2) & (y > 3e4)].index\n\ndef find_bad_rows(X, y):\n return find_bad_zeros(X, y).union(find_bad_sitezero(X)).union(find_bad_building1099(X, y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X, y = combined_data(train,weather_train)\n\nbad_rows = find_bad_rows(X, y)\n#输出异常值的index\npd.Series(bad_rows.sort_values()).to_csv(\"rows_to_drop.csv\", header=False, index=False)\n\nX = X.drop(index=bad_rows)\ny = y.reindex_like(X)\n\nX = _add_time_features(X)\nX = compress_dataframe(X)\nX = X.drop(columns=\"timestamp\") # drop掉原本的timestamp\n\ndel bad_rows,train,weather_train\ngc.collect();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.5 评价函数\n\n由于需要预测连续值，因此需要采用回归模型。由于该项目是Kaggle赛题，测试集是使用均方根对数误差 RMSLE（Root Mean Squared Logarithmic Error, RMSLE)评测的，因此这里只能使用RMSLE。RMSLE的计算公式为：\n\n$${\\rm RMSLE} = \\sqrt{\\frac{1}{n} \\sum_{i=1}^n (\\log(p_i + 1) - \\log(a_i+1))^2 }$$\n\n其中\n- $n$（public/private）数据集中的样本总数,\n- $p_i$ 是目标的预测值\n- $a_i$ 第i个目标的真实值.\n- $\\log(x)$ 是自然对数\n\n我们只需要对目标值进行$y = \\log(y+1)$的变换，就可以使用常见的RMSE作为评价函数，我们使用numpy中的log1p就可以实现。\n\n注意：进行预测的时候需要使用$y = e^y-1$将目标值转换回去，可以使用 y = np.exp1m(y)。"},{"metadata":{"trusted":true},"cell_type":"code","source":"#对目标值进行变换\ny = np.log1p(y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n### __问题 2:__\n\n思考此处为何要进行对数转换。\n\n__回答:__ "},{"metadata":{},"cell_type":"markdown","source":"\n# 4. CatBoost\n## 4.1 模型参数\n\nCatBoost主要调节的参数包括：\n- `learning_rate`：迭代步长,学习率；\n- `max_depth`：用于调节树的复杂度，较小导致欠拟合，较大导致过拟合；\n- `l2_leaf_reg`：L2正则化系数，用来控制过拟合；\n- `n_estimators`：迭代步数。"},{"metadata":{"trusted":true},"cell_type":"code","source":"cbr = CatBoostRegressor(n_estimators=100, \n use_best_model=True,\n objective='RMSE',\n loss_function='RMSE',\n boosting_type='Plain',\n eval_metric='RMSE',\n l2_leaf_reg=1,\n learning_rate=0.03,\n max_depth=8,\n task_type = 'CPU'\n )\n\n#设置分类变量\ncategorical_features=['building_id', 'site_id', 'primary_use', 'had_air_temperature', 'had_cloud_coverage', \n 'had_dew_temperature', 'had_precip_depth_1_hr','had_sea_level_pressure', 'had_wind_direction',\n 'had_wind_speed', 'tm_day_of_week', 'tm_hour_of_day']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4.2 模型训练"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_splits = 3\n\n# CatBoost的运行时间远大于LightGBM\nfor val in X['meter'].unique():\n X1 = X[X['meter'] == val].drop(columns=['meter'])\n kf = StratifiedKFold(n_splits=n_splits,random_state=42)\n t = 0\n for train_index,test_index in kf.split(X1,X1['tm_hour_of_day']):\n train_features = X1.iloc[train_index]\n train_target = y[X1.iloc[train_index].index]\n \n test_features = X1.iloc[test_index]\n test_target = y[X1.iloc[test_index].index]\n \n print(\"Building model meter :\",val,'fold:',t) \n md = cbr.fit(train_features, train_target,eval_set=(test_features, test_target),\n cat_features=categorical_features,early_stopping_rounds=500, verbose=10)\n md.save_model('cat_val{}_fold{}.bin'.format(val,t))\n t += 1\ndel train_features, test_features, X, y\ngc.collect();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n### __问题 3:__\n\n查阅资料，总结LightGBM与CatBoost的差异。\n\n__回答:__ \n"},{"metadata":{},"cell_type":"markdown","source":"\n# 5. 结果预测"},{"metadata":{"trusted":true},"cell_type":"code","source":"X = compress_dataframe(test.join(building, on=\"building_id\").join(weather_test, on=[\"site_id\", \"timestamp\"]).fillna(-1))\nX = compress_dataframe(_add_time_features(X))\nX = X.drop(columns=\"timestamp\") # drop掉原本的timestamp\n\ndel test, weather_test\ngc.collect();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#输出预测结果\nresult = np.zeros(len(X))\nfor val in X['meter'].unique():\n ix = np.nonzero((X['meter'] == val).to_numpy())\n for i in tqdm(range(n_splits)):\n #加载刚才保存的模型\n model = cbr.load_model('cat_val{}_fold{}.bin'.format(val,i))\n result[ix] += model.predict(X.iloc[ix].drop(columns=['meter']))/n_splits\n del model\n gc.collect();\n \npredictions = pd.DataFrame({\n \"row_id\": X.index,\n \"meter_reading\": np.clip(np.expm1(result), 0, None)\n})\n\n# float_format设置保留四位小数，减少文件大小，为文件上传节省时间\npredictions.to_csv(\"submission.csv\", index=False, float_format=\"%.4f\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}