<jupyter_start><jupyter_text>ASHRAE Energy Prediction> by Energetic Engineering, the tabular-1-online group in the 2021 fast.ai course- toc: true- branch: master- badges: true- comments: true- categories: [fastai, kaggle]- image: images/joey-kyber-Pihl8kTtX-s-unsplash.jpg- hide: false- search_exclude: false The [ASHRAE Great Energy Predictor III](https://www.kaggle.com/c/ashrae-energy-prediction) is a Kaggle competition to predict energy use in buildings.The training dataset is 3.2GB, containing 20 million energy meter readings for 1449 buildings over a year. Four energy meters are sampled: chilled water, electric, hot water, and steam. Each dependent measurement is given with independent variable metadata for the building and weather. The [Kaggle data page](https://www.kaggle.com/c/ashrae-energy-prediction/data) shows the available files, with histograms for each variable.Our goal was to create a model that predicts energy usage per building and per meter in that building, using fast.ai.We started by downloading the dataset from Kaggle, joining the tables, fixing timestamp inconsistencies and transforming some values into smaller datatypes to save memory. We aimed to work in Colab so memory was precious. We serialised this data preparation to a Parquet file on our Google Drives.<jupyter_code>#collapse
!pip install -Uqq fastbook dtreeviz
import fastbook
fastbook.setup_book()
#collapse
import os
import gc
import pandas as pd
import datetime as dt
import seaborn as sns
from tqdm.auto import tqdm
from pandas.api.types import is_string_dtype, is_numeric_dtype, is_categorical_dtype
from fastai.tabular.all import *
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from dtreeviz.trees import *
from IPython.display import Image, display_svg, SVG
from sklearn.metrics import mean_squared_error

#mpl.rcParams['figure.dpi'] = 100
sns.set()
pd.options.display.max_rows = 20
pd.options.display.max_columns = 8
#collapse
%cd '/content/gdrive/MyDrive/Colab Notebooks/ashrae'
train_valid = pd.read_parquet("feature_enhanced_train_combined.parquet.snappy")
#collapse
## Memory optimization

# Original code from https://www.kaggle.com/gemartin/load-data-reduce-memory-usage by @gemartin
# Modified to support timestamp type, categorical type
# Modified to add option to use float16

from pandas.api.types import is_datetime64_any_dtype as is_datetime
from pandas.api.types import is_categorical_dtype

def reduce_mem_usage(df, use_float16=False):
    """
    Iterate through all the columns of a dataframe and modify the data type to reduce memory usage.        
    """
    
    start_mem = df.memory_usage().sum() / 1024**2
    print("Memory usage of dataframe is {:.2f} MB".format(start_mem))
    
    for col in df.columns:
        if is_datetime(df[col]) or is_categorical_dtype(df[col]) or is_string_dtype(df[col]):
            continue
        col_type = df[col].dtype
        
        if col_type != object:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == "int":
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)  
            else:
                if use_float16 and c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
        else:
            df[col] = df[col].astype("category")

    end_mem = df.memory_usage().sum() / 1024**2
    print("Memory usage after optimization is: {:.2f} MB".format(end_mem))
    print("Decreased by {:.1f}%".format(100 * (start_mem - end_mem) / start_mem))
    
    return df

train_valid = reduce_mem_usage(train_valid, use_float16=True)<jupyter_output>Memory usage of dataframe is 3243.66 MB
Memory usage after optimization is: 1189.98 MB
Decreased by 63.3%<jupyter_text>Exploratory Data AnalysisWe started by browsing exploratory data analysis notebooks posted to Kaggle. We learnt that this is a gnarly dataset, with many missing values and outliers. For example, the 'meter_reading' dependent variable has mostly small values but also some very large values:<jupyter_code>train_valid[['meter_reading']].hist()<jupyter_output><empty_output><jupyter_text>Exploring further (with credit to https://www.kaggle.com/nroman/eda-for-ashrae) we discover that most of the energy in the dataset is consumed by a single building and single meter, building_id 1099 meter 2.The following plots show daily averaged meter readings for:* the full dataset* just building_id 1099 meter 2* the dataset with building_id 1099 removedIt is tempting to model this building independently from the rest of the dataset, to avoid numerical issues and overfitting to this building.<jupyter_code>#collapse
fig, axes = plt.subplots(3,1,figsize=(14, 20), dpi=100)

#train_valid[['timestamp', 'meter_reading']].set_index('timestamp').resample('H').mean()['meter_reading'].plot(ax=axes[1], alpha=0.8, label='By hour', color='tab:blue').set_ylabel('Mean meter reading', fontsize=13);
train_valid[['timestamp', 'meter_reading']].set_index('timestamp').resample('D').mean()['meter_reading'].plot(ax=axes[1], alpha=1, label='By day', color='tab:orange').set_ylabel('Mean meter reading', fontsize=13);

#train_valid[(train_valid['meter'] == 2) & (train_valid['building_id'] == 1099)][['timestamp', 'meter_reading']].set_index('timestamp').resample('H').mean()['meter_reading'].plot(ax=axes[0], alpha=0.8, label='By hour', color='tab:blue').set_ylabel('Mean meter reading', fontsize=13);
train_valid[(train_valid['meter'] == 2) & (train_valid['building_id'] == 1099)][['timestamp', 'meter_reading']].set_index('timestamp').resample('D').mean()['meter_reading'].plot(ax=axes[0], alpha=1, label='By day', color='tab:orange').set_ylabel('Mean meter reading', fontsize=13);

#train_valid[~((train_valid['meter'] == 2) & (train_valid['building_id'] == 1099))][['timestamp', 'meter_reading']].set_index('timestamp').resample('H').mean()['meter_reading'].plot(ax=axes[2], alpha=0.8, label='By hour', color='tab:blue').set_ylabel('Mean meter reading', fontsize=13);
train_valid[~((train_valid['meter'] == 2) & (train_valid['building_id'] == 1099))][['timestamp', 'meter_reading']].set_index('timestamp').resample('D').mean()['meter_reading'].plot(ax=axes[2], alpha=1, label='By day', color='tab:orange').set_ylabel('Mean meter reading', fontsize=13);

axes[0].set_title('Full dataset', fontsize=13);
axes[1].set_title('building_id==1099 and meter==2', fontsize=13);
axes[2].set_title('building_id 1099 excluded', fontsize=13);
plt.subplots_adjust(hspace=0.45)<jupyter_output><empty_output><jupyter_text>Splitting into Training and Validation setsGiven that we have a year's worth of time series data, we investigated several methods for splitting into training and validation sets. We ruled out random splits but considered periodic splits such as taking every 4th week for validation. In the end, we opted for a simple 2-fold time-based cross-validation, taking the first 6 months of data to train one model and the second 6 months for a second model, averaging predictions from each. Data TransformsWe started with the default trio of fast.ai transforms: Categorify, FillMissing, Normalize. We experimented with omitting each and saw our model performance degrade. We also experimented with dropping rows that were missing the key independent variables, since there aren't many such rows. This approach gave good performance so we would consider this when tuning our model. Feature EngineeringWe engineered some additional time features from the 'timestamp' variable, namely:* year_fraction* week_fraction* day_fractionWe found through feature clustering that 'year_fraction' was highly correlated with 'air_temperature' which makes sense because air temperature is a good proxy for month of the year. We omitted 'year_fraction' to simplify the model.<jupyter_code>#collapse
# Get a subset of the data, a few buildings and all meters

#train = train_valid[(train_valid['building_id']>=1126)]
train_valid['meter_reading'] = np.log1p(train_valid['meter_reading']) # convert target to log
train = train_valid
del train_valid
gc.collect()<jupyter_output><empty_output><jupyter_text>The code below shows how we created the DataLoaders, using only a subset of the most important features:<jupyter_code># Use an splitter and load using TabularPandas then to a dataloaders

# 2-fold cross-validation. Build 2 models, each with half the year as training data
splits1=MaskSplitter(train['timestamp'] >= '2016-07-01')(range_of(train))
cat_names = ['building_id', 'meter', 'site_id', 'primary_use']
cont_names = ['air_temperature', 'dew_temperature', 'square_feet', 'week_fraction', 'day_fraction']
procs = [Categorify, FillMissing, Normalize]
to1 = TabularPandas(train, procs=procs,
                   cat_names = cat_names,
                   cont_names = cont_names, 
                   y_names='meter_reading',
                   splits=splits1)
dls1 = to1.dataloaders(bs=1024)
del splits1, to1
gc.collect()
# Then do the same for the second fold

splits2=MaskSplitter(train['timestamp'] < '2016-07-01')(range_of(train))
cat_names = ['building_id', 'meter', 'site_id', 'primary_use']
cont_names = ['air_temperature', 'dew_temperature', 'square_feet', 'week_fraction', 'day_fraction']
procs = [Categorify, FillMissing, Normalize]
to2 = TabularPandas(train, procs=procs,
                   cat_names = cat_names,
                   cont_names = cont_names, 
                   y_names='meter_reading',
                   splits=splits2)
dls2 = to2.dataloaders(bs=1024)
del splits2, to2
gc.collect()<jupyter_output><empty_output><jupyter_text>Feature ImportanceWe trained a scikit-learn DecisionTreeRegressor to compare against our neural network and to assess feature importance.We discovered that the most significant features influencing energy usage for a particular building are:* 'square_feet' representing the size of the building* 'air_temperature' which is unsurprising* 'primary_use' surprisingly, whether a building is used for education, office, retail etc has a large influence on energy use* 'dew_temperature' is not as correlated with air temperature as you would expect* 'day_fraction' meaning that energy is consumed at regular times every day* 'week_fraction' indicating that energy use changes through the week<jupyter_code>#collapse
# Train a fast.ai decision tree
m = DecisionTreeRegressor(min_samples_leaf=2000)
m.fit(to1.train.xs, to1.train.y);
#collapse
def rf_feat_importance(m, df):
    return pd.DataFrame({'feature':df.columns, 'importance':m.feature_importances_}
                       ).sort_values('importance', ascending=False)

fi = rf_feat_importance(m, to1.train.xs)
fi[:10]
#collapse

print("Training set RMSLE: " + str(mean_squared_error(to1.train.y, m.predict(to1.train.xs), squared=False)))
print("Validation set RMSLE: " + str(mean_squared_error(to1.valid.y, m.predict(to1.valid.xs), squared=False)))

building=2
meter=2
#plt.scatter(valid_xs[(valid_xs['building_id']==building) & (valid_xs['meter']==meter)]['air_temperature'], valid_y.loc[valid_xs[(valid_xs['building_id']==building) & (valid_xs['meter']==meter)].index])
#plt.scatter(valid_xs[(valid_xs['building_id']==building) & (valid_xs['meter']==meter)]['air_temperature'], pd.DataFrame(valid_y_pred, index=valid_y.index).loc[valid_xs[(valid_xs['building_id']==building) & (valid_xs['meter']==meter)].index], color='orange')<jupyter_output>Training set RMSLE: 1.0176237141131312
Validation set RMSLE: 1.4695399113467222<jupyter_text>Learning RatesWe tried using lr_find() to choose a learning rate. However, the curve is upside down.<jupyter_code># Train a tabular model
learn1 = tabular_learner(dls1, metrics=rmse)
#learn1.summary()

# Exploring learning rates
learn1.lr_find(stop_div=False, num_it=200)

learn2 = tabular_learner(dls2, metrics=rmse)<jupyter_output><empty_output><jupyter_text>Model TrainingSo, we trained a bunch of models, experimenting with all the hyperparameters discussed, including:* Features* Batch sizes* Learning rates* Weight decay* Number of buildings included<jupyter_code>learn1.fit_one_cycle(5, lr_max=.00001)
#collapse
learn2.fit_one_cycle(5, lr_max=.00001)
#collapse 
# Pickling didn't work due to an unspecified google drive error
with open('learn1.pkl','wb') as outfile:
  save_pickle(outfile, learn1)
#collapse
infile = open('learn1.pkl','rb')
learn1 = load_pickle(infile)
#collapse 
# Pickling didn't work due to an unspecified google drive error
with open('learn2.pkl','wb') as outfile:
  save_pickle(outfile, learn2)
#collapse
infile = open('learn2.pkl','rb')
learn2 = load_pickle(infile)<jupyter_output><empty_output><jupyter_text>MetricsThe competition stated that that chosen metric was root mean squared log error: RMSLE=sqrt(1/𝑛∑(log(𝑝+1)−log(𝑎+1))**2)fast.ai doesn't have this metric but the metric can be achieved by taking the log of the dependent variable, 'meter_reading' and then using root mean squared error (RMSE).We found that using the log of 'meter_reading' had the added bonus of much better performance and more stable training. This is because 'meter_reading' had a large dynamic range, with some small values and some very large values. Taking the log is a form of normalisation, which is more numerically stable. We achieved a validation set RMSLEs of 1.38 and 1.52 for our two models, compared to validation set RMSLE of 1.47 for our DecisionTreeRegressor. The winning Kaggle entry scored RMSLE of 1.23. Assuming that our cross-validated sets are representative of the test set (a big assumption!) this puts us around 2700th on the Kaggle leaderboard of 3500 entries. That's reasonable given our first attempt with no data cleaning.<jupyter_code>from sklearn.metrics import mean_squared_error

#collapse
if 'to1' in locals():
  del to1
if 'dls1' in locals():
  del dls1

# y_targs1 is the same as to1.train.y
# valid_y_targs1 is the same as to1.valid.y

y_pred1,y_targs1 = learn1.get_preds(ds_idx=0) # training set predictions
valid_y_pred1,valid_y_targs1 = learn1.get_preds() # validation set predictions
print("Model 1 training set RMSLE: " + str(mean_squared_error(y_pred1.detach().numpy(), y_targs1, squared=False)))
print("Model 1 Validation set RMSLE: " + str(mean_squared_error(valid_y_pred1.detach().numpy(), valid_y_targs1, squared=False)))
if 'to2' in locals():
  del to2
if 'dls2' in locals():
  del dls2

y_pred2,y_targs2 = learn2.get_preds(ds_idx=0) # training set predictions
valid_y_pred2,valid_y_targs2 = learn2.get_preds() # validation set predictions
print("Model 2 training set RMSLE: " + str(mean_squared_error(y_pred2.detach().numpy(), y_targs2, squared=False)))
print("Model 2 Validation set RMSLE: " + str(mean_squared_error(valid_y_pred2.detach().numpy(), valid_y_targs2, squared=False)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11190522181443564, 'fraction_numerical': 0.02440280728864851, 'mean_word_length': 5.1610472034906785, 'pattern_counts': {'":': 1, '<': 42, '<?xml version=': 0, '>': 44, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '478600', 'license': 'Apache-2.0', 'n_tokens_mistral': 5219, 'n_tokens_neox': 4773, 'n_words': 1438, 'path': '_notebooks/2021-07-23-tabular1online-presentation.ipynb', 'repo_name': 'fastai-energetic-engineering/ashrae', 'size': '15531'}