{"cells":[{"metadata":{},"cell_type":"markdown","source":"# **ASHRAE III - Metering Data Analysis - EXPLORE & TRAIN**\nAs an engineer with a background in Building Automation, I have some first-hand experience dealing with building metering data, particularly the analysis of it for reducing energy/water consumption. Some insights from wider industry around metering data:\n* If a meter is offline, it will generally read zero consumption (as opposed to empty values)\n* Water meters and steam meters generate consumption data by 'pulsing' every time it accumulates a certain volume of water; the L/pulse calibration needs to be set correctly otherwise the meter will read potentially orders of magnitude out (e.g. is it 1 L/pulse vs 10 L/pulse vs 1m3/pulse)\n* It is common for a building's meter profile to change with variances in tenancy occupation and building construction works\n* Public holidays will have signiciant impacts on building meter profiles","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Load Data\nFirst we start with standard imports and merging metering readings with the building meta data and weather data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import dependencies\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport xgboost as xgb\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_log_error\nfrom sklearn.metrics import mean_squared_error\nimport gc\nimport os\nfrom pandas.tseries.holiday import USFederalHolidayCalendar as calendar\n\n# Define file + folder names\nfolder = '/kaggle/input/ashrae-energy-prediction/'\ntrain_filename = 'train.csv'\nweather_train_filename = 'weather_train.csv'\nbuilding_meta_data_filename = 'building_metadata.csv'\ntest_filename = 'test.csv'\nsample_submission_filename = 'sample_submission.csv'\nweather_test_filename = 'weather_test.csv'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# Load train data into dataframes\ntrain = pd.read_csv(folder + train_filename)\nweather_train = pd.read_csv(folder + weather_train_filename)\nbuilding_meta_data = pd.read_csv(folder + building_meta_data_filename)\n\n# Merge data\ntrain = train.merge(building_meta_data, left_on = 'building_id', right_on = 'building_id')\ntrain = train.merge(weather_train, left_on = ['site_id','timestamp'], right_on = ['site_id','timestamp'])\n\n# Delete redundant DataFrames from memory\ndel weather_train, building_meta_data\n\ntrain.head()\ntrain.info()\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Convert Categorical Information\n\nAs per above data information, the only categorical variable requiring encoding is a building's 'primary_use'. The below shows the counts of buildings according to their primary use.\n\nWe then convert the 'primary_use' feature to be a numeric variable","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# show counts of buildings according to their primary use\nprint(train.groupby('primary_use')['building_id'].unique().apply(len).sort_values(ascending=False))\n\n# convert primary_use (type 'object') to numerical\ndef convert_to_numeric(df,cols):\n    df_cats = pd.DataFrame()\n    for feature in cols:\n        df[feature] = df[feature].astype('category').cat.codes\n        df_cats[feature] = dict(enumerate(df[feature].astype('category').cat.categories))\n        print(feature)\n    return df, df_cats\n\ntrain, category_dict = convert_to_numeric(train, ['primary_use'])\n\ncategory_dict","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Date-Time Features\n\nPerhaps most importantly we need to create date & time features. Some speculations on how these features can impact the accuracy include:\n* Month: Changes in building setpoints during winter months to account for people wearing warmer clothes\n* Hour: Occupancy hours will differ for each building\n* Day Of Week & Holiday: Some buildings will be occupied on the weekend, whilst others will not","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# Add new features for Date & Time\ntrain['timestamp'] = pd.to_datetime(train['timestamp'], format='%Y-%m-%d %H:%M:%S')\ntrain['date'] = train['timestamp'].dt.date\ntrain['year'] = train['timestamp'].dt.year\ntrain['month'] = train['timestamp'].dt.month\ntrain['hour'] = train['timestamp'].dt.hour\ntrain['week'] = train['timestamp'].dt.week\ntrain['day_of_week'] = train['timestamp'].dt.dayofweek\ntrain['week_day'] = [0 if day==5 or day==6 else 1 for day in train['day_of_week']]\n\n# Add feature for US holidays\ndates_range = pd.date_range(start='2015-12-31', end='2019-01-01')\nus_holidays = calendar().holidays(start=dates_range.min(), end=dates_range.max())\ntrain['holiday'] = (train['timestamp'].dt.date.astype('datetime64').isin(us_holidays)).astype(np.int8)\n    \ndel train['timestamp']\n  \nprint(train.columns)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Reduce Memory Consumption\n\nThe memory usage of the data is reduced using kernel from https://www.kaggle.com/arjanso/reducing-dataframe-memory-size-by-65.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Based on this great kernel https://www.kaggle.com/arjanso/reducing-dataframe-memory-size-by-65\ndef reduce_mem_usage(df):\n    start_mem_usg = df.memory_usage().sum() / 1024**2 \n    print(\"Memory usage of properties dataframe is :\",start_mem_usg,\" MB\")\n    NAlist = [] # Keeps track of columns that have missing values filled in. \n    for col in df.columns:\n        if df[col].dtype != object:  # Exclude strings            \n            # Print current column type\n            print(\"******************************\")\n            print(\"Column: \",col)\n            print(\"dtype before: \",df[col].dtype)            \n            # make variables for Int, max and min\n            IsInt = False\n            mx = df[col].max()\n            mn = df[col].min()\n            print(\"min for this col: \",mn)\n            print(\"max for this col: \",mx)\n            # Integer does not support NA, therefore, NA needs to be filled\n            if not np.isfinite(df[col]).all(): \n                NAlist.append(col)\n                df[col].fillna(mn-1,inplace=True)  \n                   \n            # test if column can be converted to an integer\n            asint = df[col].fillna(0).astype(np.int64)\n            result = (df[col] - asint)\n            result = result.sum()\n            if result > -0.01 and result < 0.01:\n                IsInt = True            \n            # Make Integer/unsigned Integer datatypes\n            if IsInt:\n                if mn >= 0:\n                    if mx < 255:\n                        df[col] = df[col].astype(np.uint8)\n                    elif mx < 65535:\n                        df[col] = df[col].astype(np.uint16)\n                    elif mx < 4294967295:\n                        df[col] = df[col].astype(np.uint32)\n                    else:\n                        df[col] = df[col].astype(np.uint64)\n                else:\n                    if mn > np.iinfo(np.int8).min and mx < np.iinfo(np.int8).max:\n                        df[col] = df[col].astype(np.int8)\n                    elif mn > np.iinfo(np.int16).min and mx < np.iinfo(np.int16).max:\n                        df[col] = df[col].astype(np.int16)\n                    elif mn > np.iinfo(np.int32).min and mx < np.iinfo(np.int32).max:\n                        df[col] = df[col].astype(np.int32)\n                    elif mn > np.iinfo(np.int64).min and mx < np.iinfo(np.int64).max:\n                        df[col] = df[col].astype(np.int64)    \n            # Make float datatypes 32 bit\n            else:\n                df[col] = df[col].astype(np.float32)\n            \n            # Print new column type\n            print(\"dtype after: \",df[col].dtype)\n            print(\"******************************\")\n    # Print final result\n    print(\"___MEMORY USAGE AFTER COMPLETION:___\")\n    mem_usg = df.memory_usage().sum() / 1024**2 \n    print(\"Memory usage is: \",mem_usg,\" MB\")\n    print(\"This is \",100*mem_usg/start_mem_usg,\"% of the initial size\")\n    return df, NAlist\n\ntrain, nalist = reduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data Exploration & Visualization\n\nTo ensure we have an even distribution of readings we visualise the frequency histogram of the various date & time features.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualize distribution of Date-Time features\ndatetime_features = ['year','month','hour','day_of_week']\nfor idx, feature in enumerate(datetime_features):\n    fig, axs = plt.subplots(1,1,constrained_layout=True)\n    axs.hist(train[feature],bins=len(np.unique(train[feature])))\n    axs.set_title(feature)\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The distribution of meter readings is plotted per meter type, which approximately follows a power-law relationship. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Visualize distribution of meter readings\nmeter_types = {0: \"electricity\",\n               1: \"chilledwater\",\n               2: \"steam\",\n               3: \"hotwater\"}\nfor meter_idx in meter_types:\n    fig, axs = plt.subplots(1,1,constrained_layout=True)\n    axs.hist(train[train['meter']==meter_idx]['meter_reading'],bins=1000)\n    axs.set_title(meter_types[meter_idx])\n    axs.set_yscale('log')\n    axs.set_xscale('log')\n    axs.set_xlabel('meter_readings')\n    axs.set_ylabel('count')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To get a general idea of the yearly trends per meter type, we group by day and then plot the median values for each day. Using the median helps to be more resilient to some of the higher consumptions exhibited by outlier buildings. The following insights can be gleaned:\n* chilledwater, hotwater and steam are all heavily dependent on seasonality\n* electricity has a fairly consistent weekly profile","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Group to find daily consumption (kWh), then group by building_id to find median daily consumption\ntrain_median_daily_consumption = train.groupby(['meter','date','building_id']).sum()['meter_reading'].groupby(level=[0,1]).median()\nmeters_daily = {0: train_median_daily_consumption.loc[0,:].droplevel(0),\n                1: train_median_daily_consumption.loc[1,:].droplevel(0),\n                2: train_median_daily_consumption.loc[2,:].droplevel(0),\n                3: train_median_daily_consumption.loc[3,:].droplevel(0)}\ndel train_median_daily_consumption\n\n# Plot median daily values\nfor meter_idx in meter_types:\n    fig, axs = plt.subplots(1,1,constrained_layout=True)\n    axs.plot(meters_daily[meter_idx])\n    axs.set_title('Median Daily Meter Readings (' + str(meter_types[meter_idx]) + ')')\n    axs.set_xlabel('date')\n    axs.set_ylabel('meter reading (kWh)')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data Cleansing","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Cleaning the data readings firstly requires addressing any zero values according to the following conditions based on meter type:\n 1. Electricity - if zero for that day, then remove that day\n 2. Chilled Water - if zero for that week, and week is in the winter months (i.e. week<10 or week>47)\n 3. Steam - if zero for that week, and week is in the summer months (i.e. 22<week<35)\n 4. Hot Water - if zero for that week, and week is in the summer months (i.e. 22<week<35)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define function for handling advanced filtering by multiple columns\ndef advanced_filter(df, conditions, invert=False):\n    cols = conditions.columns\n    df_as_str = df[cols].astype(str).sum(axis = 1) # cast as string\n    conditions_as_str = conditions.astype(str).sum(axis = 1) # cast as string\n    return df[df_as_str.isin(conditions_as_str)!=invert]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get multi-index of date + building_id where daily consumption is zero\nelectricity_daily = train[train['meter']==0].groupby(['meter','date','building_id']).sum()['meter_reading']\noffline_electricity_meters = electricity_daily[electricity_daily==0]\nprint('Percent of days where electricity meter deemed to be offline: ' + str(round(offline_electricity_meters.size/electricity_daily.size*100,1)) + '%')\noffline_electricity_meters = offline_electricity_meters.reset_index().drop('meter_reading',axis=1)\noffline_electricity_meters.head()\n\n# Apply filter\ntrain = advanced_filter(train, offline_electricity_meters, invert=True)\n\n# Get multi-index of week + building_id where weekly consumption is zero\nmeters_weekly = train.groupby(['meter','week','building_id']).sum()['meter_reading']\noffline_meters_weekly = meters_weekly[meters_weekly==0]\nprint('Percent of weeks where meter reading found to be offline: ' + str(round(offline_meters_weekly.size/meters_weekly.size*100)) + '%') # Note this may also include electricity meters offline days\noffline_meters_weekly = offline_meters_weekly.reset_index().drop('meter_reading',axis=1)\noffline_meters_weekly.head()\n\n# Apply filter\ntrain = advanced_filter(train, offline_meters_weekly, invert=True)\n    \ntrain\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Determine outlier buildings that will require separate models built","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Group by building_id and calculate mean meter reading\nmeter_readings_by_building_id = train.groupby(['building_id','meter']).mean()['meter_reading'].sort_values(ascending=False)\n# Print first few rows per meter type\nfor meter_idx in meter_types:\n    print('\\n' + meter_types[meter_idx] + 'readings: ')\n    print(meter_readings_by_building_id.loc[:,meter_idx].head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot mean meter read per building\nfor meter_idx in meter_types:\n    fig, axs = plt.subplots(1,1)\n    axs.bar(meter_readings_by_building_id.loc[:,meter_idx].index, meter_readings_by_building_id.loc[:,meter_idx], width=5)\n    axs.set_title('mean ' + meter_types[meter_idx] + ' meter readings per building')\n    axs.set_xlabel('building_id')\n    axs.set_ylabel('mean meter reading (kWh)')\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As a result of the above breakdown showing evidence that some buildings have readings orders of magnitude above the rest, the following buildings will have their own separate models built:\n* building 778: chilledwater\n* building 1099: steam\n* building 1021: hotwater","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define outlier buildings & associated meters\noutlier_meters = pd.DataFrame({'building_id': [778, 1099, 1021],\n                               'meter': [1, 2, 3]})\n\n# Create training data for these outlier meters\noutlier_train = advanced_filter(train, outlier_meters)\n\n# Remove from the other train data\ntrain = advanced_filter(train, outlier_meters, invert=True)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Convert site 0 electrical meter from kBTU to kWh for training as per https://www.kaggle.com/c/ashrae-energy-prediction/discussion/119261 (Multiply by 0.2931).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"site_0_elec_bool = (train['meter']==0) & (train['site_id']==0)\ntrain.loc[site_0_elec_bool]['meter_reading'] = train.loc[site_0_elec_bool]['meter_reading']*0.2931","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define X / Y columns\nX_features = ['building_id', 'meter', 'site_id', 'primary_use', 'square_feet', 'year_built', 'floor_count',  'air_temperature', 'dew_temperature', 'precip_depth_1_hr', 'month', 'hour', 'day_of_week', 'holiday']\ny_target = ['meter_reading']\n# Create data for main model\nX_train = train[X_features]\ny_train = train[y_target]\n# Create data for outlier model\nX_outlier_train = outlier_train[X_features]\ny_outlier_train = outlier_train[y_target]\n# Remove redundant data\ndel train, outlier_train","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Step 1: Find optimum number of trees for some base initial conditions**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Define XGBoost Wrapper for Cross-Validation\ndef xgb_wrapper(learning_rate, max_depth, gamma, n_estimators, subsample):\n    # Split data for cross-validation\n    kf = KFold(n_splits=3, random_state=0, shuffle=True)\n    # Loop through splits\n    rmsle = []\n    for train_split_index, test_split_index in kf.split(X_train):\n        # Create Data Matrixes\n        dtrain_split = xgb.DMatrix(data=X_train.iloc[train_split_index], label=y_train.iloc[train_split_index])\n        dtest_split = xgb.DMatrix(X_train.iloc[test_split_index])\n        # Define XGBRegressor parameters\n        xgb_params = {\n            'eval_metric': 'rmse',\n            'tree_method': 'gpu_hist',\n            \"learning_rate\": learning_rate,\n            \"max_detph\": max_depth,\n            \"gamma\": gamma,\n            \"subsample\": subsample}\n        # Train XGBoost Model & make predictions\n        print('\\nCreating new model: ')\n        print(xgb_params)\n        xgb_reg = xgb.train(xgb_params, dtrain_split, n_estimators)\n        predictions = xgb_reg.predict(dtest_split)\n        predictions[predictions < 0] = 0 # set negative predictions to zero\n        # Calculate Evaluation Metrix (Root Mean Squared Logarithmic Error - RMSLE)\n        model_rmsle = np.sqrt(mean_squared_log_error(y_train.iloc[test_split_index], predictions))\n        print('RMSLE: '+ str(model_rmsle))\n        rmsle.append(model_rmsle)\n        # Clear memory\n        del dtrain_split, dtest_split, xgb_reg\n        gc.collect()\n    # Return evaluation metrix\n        print(rmsle)\n    cv_result = np.array(rmsle).mean()\n    return cv_result\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cv_results = {}\n# Find optimum n_estimators\nn_estimators_range = list(range(100,1001,100))\nfor n_estimators in n_estimators_range:\n    gc.collect()\n    print('\\n****************** Performing CV for n_estimators=' + str(n_estimators)+' ******************')\n    cv_result = xgb_wrapper(0.1, 5, 0.1, n_estimators, 0.8)\n    cv_results.update({n_estimators: cv_result})\n    print(cv_results)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nplt.plot(list(cv_results.keys()),list(cv_results.values()))\nplt.xlabel('n_estimators')\nplt.ylabel('rmse')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Step 2: Find optimum max depth & min_child_weight for some best n_estimators=1,000**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Find optimum value for max_depth\nmax_depth_range = list(range(3,7))\nmax_depth_cv_result = {}\nfor max_depth in max_depth_range:\n    gc.collect()\n    print('\\n****************** Performing CV for max_depth=' + str(max_depth)+' ******************')\n    cv_result = xgb_wrapper(0.1, max_depth, 0.1, 1000, 0.8)\n    max_depth_cv_result.update({max_depth: cv_result})\n    print(max_depth_cv_result)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nplt.plot(list(max_depth_cv_result.keys()),list(max_depth_cv_result.values()))\nplt.xlabel('max_depth')\nplt.ylabel('rmse')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Find optimum value for min_child_weight\ngamma_range = list(range(0,21,5))\ngamma_cv_result = {}\nfor gamma in gamma_range:\n    gamma=gamma/100\n    gc.collect()\n    print('\\n****************** Performing CV for gamma=' + str(gamma)+' ******************')\n    cv_result = xgb_wrapper(0.1, 3, gamma, 1000, 0.8)\n    gamma_cv_result.update({gamma: cv_result})\n    print(gamma_cv_result)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nplt.plot(list(gamma_cv_result.keys()),list(gamma_cv_result.values()))\nplt.xlabel('gamma')\nplt.ylabel('rmse')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Step 3: Drop learning_rate and increase n_estimators and train on whole data set**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()\n# Train a model on whole data set\ndtrain = xgb.DMatrix(data=X_train, label=y_train)\n# Define XGBRegressor parameters\nxgb_params = {\n    'eval_metric': 'rmse',\n    'tree_method': 'gpu_hist',\n    \"learning_rate\": 0.05,\n    \"max_detph\": 3,\n    \"gamma\": 0.15,\n    \"subsample\": 0.8}\n# Train XGBoost Model & make predictions\nxgb_reg = xgb.train(xgb_params, dtrain, 3000)\n# Save model\nxgb_reg.save_model('0001.model')\n\n# Train outlier model\ndtrain = xgb.DMatrix(data=X_outlier_train, label=y_outlier_train)\nxgb_outlier_reg = xgb.train(xgb_params, dtrain, 3000)\nxgb_outlier_reg.save_model('0001_outlier.model')\n","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}