{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Problem statement\n\nQ: How much does it cost to cool a skyscraper in the summer?\nA: A lot! And not just in dollars, but in environmental impact.\n\nThankfully, significant investments are being made to improve building efficiencies to reduce costs and emissions. The question is, are the improvements working? That’s where you come in. Under pay-for-performance financing, the building owner makes payments based on the difference between their real energy consumption and what they would have used without any retrofits. The latter values have to come from a model. Current methods of estimation are fragmented and do not scale well. Some assume a specific meter type or don’t work with different building types.\n\nIn this competition, you’ll develop accurate models of metered building energy usage in the following areas: chilled water, electric, hot water, and steam meters. The data comes from over 1,000 buildings over a three-year timeframe. With better estimates of these energy-saving investments, large scale investors and financial institutions will be more inclined to invest in this area to enable progress in building efficiencies."},{"metadata":{},"cell_type":"markdown","source":"### Note about codestyle\n\nIn Jupyter notebooks, the emphasis is on quick experimentation. The quality of the code is not what we optimize for. So when \"productionizing\" this notebook, take everything with a grain of salt and rethink the structure of the code."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport lightgbm as lgb\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import GroupKFold\n\nsns.set_style('darkgrid')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Input data is available on the file system in `../input/ashrae-energy-prediction`. Let's just list it first."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"DATA_DIR = '../input/ashrae-energy-prediction'\nprint(os.listdir(DATA_DIR))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Description of files\n\n(pasted from https://www.kaggle.com/c/ashrae-energy-prediction/data)\n\n### train.csv\n* `building_id` - Foreign key for the building metadata.\n* `meter` - The meter id code. Read as `{0: electricity, 1: chilledwater, 2: steam, 3: hotwater}`. Not every building has all meter types.\n* `timestamp` - When the measurement was taken\n* `meter_reading` - The target variable. Energy consumption in kWh (or equivalent). Note that this is real data with measurement error, which we expect will impose a baseline level of modeling error. UPDATE: as discussed here, the site 0 electric meter readings are in kBTU.\n\n### building_meta.csv\n* `site_id` - Foreign key for the weather files.\n* `building_id` - Foreign key for training.csv\n* `primary_use` - Indicator of the primary category of activities for the building based on EnergyStar property type definitions\n* `square_feet` - Gross floor area of the building\n* `year_built` - Year building was opened\n* `floor_count` - Number of floors of the building\n\n### weather_[train/test].csv\nWeather data from a meteorological station as close as possible to the site.\n\n* `site_id`\n* `air_temperature` - Degrees Celsius\n* `cloud_coverage` - Portion of the sky covered in clouds, in oktas\n* `dew_temperature` - Degrees Celsius\n* `precip_depth_1_hr` - Millimeters\n* `sea_level_pressure` - Millibar/hectopascals\n* `wind_direction` - Compass direction (0-360)\n* `wind_speed` - Meters per second\n\n### test.csv\nThe submission files use row numbers for ID codes in order to save space on the file uploads. `test.csv` has no feature data; it exists so you can get your predictions into the correct order.\n\n* `row_id` - Row id for your submission file\n* `building_id` - Building id code\n* `meter` - The meter id code\n* `timestamp` - Timestamps for the test data period"},{"metadata":{},"cell_type":"markdown","source":"# Load up the data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Function to reduce the DF size\ndef reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)  \n            else:\n                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                    df[col] = df[col].astype(np.float16)\n                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)    \n    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df\n\ndef load_df(fname):\n    df = pd.read_csv(os.path.join(DATA_DIR, fname))\n    if 'timestamp' in df.columns:\n        # I guess fortunately all timestamp columns are called `timestamp`.\n        df['timestamp'] = pd.to_datetime(df['timestamp'])\n    return df\n\ntrain_df = reduce_mem_usage(load_df('train.csv'))\nbuilding_metadata_df = reduce_mem_usage(load_df('building_metadata.csv'))\nweather_train_df = reduce_mem_usage(load_df('weather_train.csv'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Exploratory data analysis, minor cleanup and feature generation"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def describe_df(df):\n    print('Shape of data: ', df.shape)\n    print('\\nBasic info:')\n    print(df.info())\n    print('\\nQuick peek at the data:')\n    print(df.head())\n    print('\\nBasic description of the data:')\n    print(df.describe())\n    print('\\nLooking at NAs')\n    print(df.isna().sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Analysis: `train_df`"},{"metadata":{"trusted":true},"cell_type":"code","source":"describe_df(train_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are no NAs in `train_df`, so that's nice. The timestamps have been treated well. The memory consumption is modest, so there's no need to mess around with that.\n\nIt's worth looking at the meter readings as timeseries data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# The duration of the training data\ntrain_df['timestamp'].min(), train_df['timestamp'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.relplot(x='timestamp', y='meter_reading', hue='meter', kind='line', \n            palette=sns.color_palette('hls', 4), aspect=16/9,\n            data=(train_df[train_df['timestamp'] > pd.to_datetime('2016-12-01')]\n                  .groupby(by=['meter', 'timestamp'])\n                  .agg({'meter_reading': 'median'}).reset_index()))\nplt.xticks(rotation=15)\nplt.title('Median meter readings over time for different meter types (1 month)')\nplt.gca().set(yscale='log')\nplt.ylabel('Meter reading (log-scale)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Note: The `meter_reading` column is in log-scale above. The outputs of meter type `2` are much higher than the rest, so it's easier to see patterns in log-scale.\n\nIn the plot above, we can see that the meter type is obviously an important feature. In addition, clearly the `meter_reading` has plenty of seasonality. The most evident examples of seasonality here are based on time of day and day of week and that seasonality is different for different meter types.\n\nFor example, notice that meter `0` is clearly affected by weekend vs. weekday dynamics more so than other meter types. There are probably also monthly seasonal effects based on heating needs being different in winter vs. summer, although that's not visible in the graph above. To see that aspect of seasonality look at the plot below."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.relplot(x='month_of_year', y='meter_reading', hue='meter', kind='line', \n            palette=sns.color_palette('hls', 4), aspect=16/9,\n            data=(train_df\n                  .assign(month_of_year=train_df['timestamp'].dt.month)\n                  .groupby(by=['meter', 'month_of_year'])\n                  .agg({'meter_reading': 'median'}).reset_index()))\nplt.xticks(rotation=15)\nplt.title('Median meter readings over months for different meter types (1 year)')\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So far, we've concluded that the meter type, the hour of day, the day of week, the day of year are potentially useful features. So let's add those to the `train_df` dataframe.\n\n## Feature generation: `train_df`"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df.assign(hour_of_day=train_df['timestamp'].dt.hour, \n                           day_of_week=train_df['timestamp'].dt.dayofweek,\n                           day_of_year=train_df['timestamp'].dt.dayofyear)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Analysis: `building_metadata.csv`\n\nThe `train_df` dataframe had a `building_id` column that we never investigated. Let's do that now, combined with the `building_metadata.csv` file."},{"metadata":{"trusted":true},"cell_type":"code","source":"describe_df(building_metadata_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One of the problems here is that the columns `year_built` and `floor_count` have plenty of null values. For now, I'm planning to leave them as is. Lightgbm can handle null values, so we'll rely on that for now. \n\nIt would be interesting to see how both `square_feet`, `primary_use` and `year_built` influence `energy_consumption`"},{"metadata":{},"cell_type":"markdown","source":"### Relationship between `square_feet` and energy consumption"},{"metadata":{"trusted":true},"cell_type":"code","source":"g = sns.jointplot(x='square_feet', y='meter_reading', height=8,\n                  data=(train_df.groupby(by='building_id')\n                        .agg({'meter_reading': 'median'})\n                        .join(building_metadata_df, on=['building_id'])))\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's interesting that there doesn't seem to be much of a linear relationship between `square_feet` and the `meter_reading`. What's more interesting is that the relationship is more evident in log-log space."},{"metadata":{"trusted":true},"cell_type":"code","source":"g = sns.jointplot(x='square_feet', y='meter_reading', height=8, kind='reg',\n                  data=(train_df.groupby(by='building_id')\n                        .agg({'meter_reading': 'median'})\n                        .join(building_metadata_df, on=['building_id'])\n                        .pipe(lambda df: df.assign(meter_reading=np.log1p(df.meter_reading),\n                                                   square_feet=np.log1p(df.square_feet)))))\ng.ax_joint.set_xlabel('Meter reading (log-scale)')\ng.ax_joint.set_ylabel('Square feet (log-scale)')\ng.fig.suptitle('Relationship between Square footage and energy consumption in log-log space')\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Relationship between `primary_use` and energy consumption"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(10, 20))\nsns.violinplot(x='meter_reading', y='primary_use', orient='h', scale='count',\n               data=(train_df.groupby(by='building_id')\n                     .agg({'meter_reading': 'median'})\n                     .join(building_metadata_df, on=['building_id'])))\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As could have been expected, the distributions of energy consumption for various types of buildings is different. Since the violin plots are scaled by the count of datapoints, it means that most of the data in the training set is for education buildings.\n\n### Relationship between `year_built` and energy consumption"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10, 5))\nsns.lineplot(x='year_built', y='meter_reading',\n             data=(train_df\n                   .merge(building_metadata_df, on=['building_id'])\n                   .groupby(by='year_built')\n                   .agg({'meter_reading': 'median'})\n                   .reset_index()))\nplt.title('Meter readings for buildings built in different years')\nplt.show()\nplt.close()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Clearly, the median meter reading for buildings built during different years is very different. So of course the year_built is an important feature. Perhaps it would be appropriate to keep it a categorical feature. "},{"metadata":{},"cell_type":"markdown","source":"## Feature generation: `building_metadata.csv`\n\nFor now, we definitely want to encode the 2 categorical variables: `year_built` and `primary_use`. We also want to log-transform the `square_feet` column because of our previous analysis."},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata_enc = {\n    'year_built': LabelEncoder(),\n    'primary_use': LabelEncoder(),\n}\nbuilding_metadata_df['year_built_enc'] = (building_metadata_enc['year_built']\n                                          .fit_transform(building_metadata_df['year_built']))\nbuilding_metadata_df['primary_use_enc'] = (building_metadata_enc['primary_use']\n                                           .fit_transform(building_metadata_df['primary_use']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata_df['square_feet_log'] = np.log1p(building_metadata_df['square_feet'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Analysis: `weather_train.csv`"},{"metadata":{"trusted":true},"cell_type":"code","source":"describe_df(weather_train_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are lot of nans here. We might have to do some imputation for weather."},{"metadata":{},"cell_type":"markdown","source":"# Prepare training data\n\nFirst, we'll just merge the different dataframes so it can be used for model training."},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_df = (train_df\n             .merge(building_metadata_df, on=['building_id'])\n             .merge(weather_train_df, on=['site_id', 'timestamp']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We established earlier that there is value in dealing with the meter_reading in log-space. So we'll do that transformation here as well."},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_df['meter_reading_log'] = np.log1p(merged_df['meter_reading'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Of course there are plenty of columns in the `merged_train_df` dataframe that are not required for the actual modeling. So let's select the columns we care for. Let's also make a separate dataframe for testing."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(sorted(merged_df.columns))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_df = merged_df.sort_values('timestamp')\nfeature_cols = ['building_id', 'day_of_week', 'day_of_year', 'floor_count', \n                'hour_of_day', 'meter', 'primary_use_enc', 'site_id', \n                'square_feet_log', 'year_built_enc']\ncategorical_features = ['building_id', 'day_of_week', 'day_of_year', 'hour_of_day',\n                        'meter', 'primary_use_enc', 'site_id', 'year_built_enc']\nX_df = merged_df[feature_cols]\ny_df = merged_df[['meter_reading_log']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# LightGBM model"},{"metadata":{},"cell_type":"markdown","source":"## Train/test split methodology\n\nThe training dataset ranges from `2016-01-01 00:00:00` to `2016-12-31 23:00:00` and the test dataset ranges from `2017-01-01 00:00:00` to `2018-12-31 23:00:00`.\n\nSo if we split the data for our own train/validation routine, we should also split the data with non-overlapping timestamps. One way to think about it is that the model should learn the day-to-day changes in energy consumption from the seasonality and weather features and not from knowing the changes in energy consumption for the same timestamps for other buildings in the same site."},{"metadata":{"trusted":true},"cell_type":"code","source":"# We can use `sklearn.model_selection.GroupKFold` and define \"groups\" to be based on\n# the month. Our training dataset has 12 months, and we can choose 3 splits.\nkfold = GroupKFold(n_splits=3)\ngroups = merged_df['timestamp'].dt.month","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\n\n# These are the usual ipython objects, including this one you are creating\nipython_vars = ['In', 'Out', 'exit', 'quit', 'get_ipython', 'ipython_vars']\n\n# Get a sorted list of the objects and their sizes\nsorted([(x, sys.getsizeof(globals().get(x))) for x in dir() if not x.startswith('_') and x not in sys.modules and x not in ipython_vars], key=lambda x: x[1], reverse=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get rid of other dataframes and save some memory before the real show begins\ndel train_df\ndel weather_train_df\ndel merged_df\n\nimport gc\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import sys\n\n# These are the usual ipython objects, including this one you are creating\nipython_vars = ['In', 'Out', 'exit', 'quit', 'get_ipython', 'ipython_vars']\n\n# Get a sorted list of the objects and their sizes\nsorted([(x, sys.getsizeof(globals().get(x))) for x in dir() if not x.startswith('_') and x not in sys.modules and x not in ipython_vars], key=lambda x: x[1], reverse=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"models = []\nparams = {\n    \"objective\": \"regression\",\n    \"boosting\": \"gbdt\",\n    \"num_leaves\": 1280,\n    \"learning_rate\": 0.05,\n    \"feature_fraction\": 0.85,\n    \"reg_lambda\": 2,\n    # The actual metric we'd be measured against is RMLSE \n    # (https://www.kaggle.com/c/ashrae-energy-prediction/overview/evaluation)\n    # but since we've already taken the log of the meter reading as the target\n    # we should just use the RMSE metric.\n    \"metric\": \"rmse\",\n}\nfor idx, (train_index, val_index) in enumerate(kfold.split(X_df, y_df, groups)):\n    print(f'Training fold {idx}')\n    X_train_df, y_train_df = X_df.loc[train_index], y_df.loc[train_index]\n    X_val_df, y_val_df = X_df.loc[val_index], y_df.loc[val_index]\n    train_dataset = lgb.Dataset(X_train_df, label=y_train_df, \n                                categorical_feature=categorical_features)\n    val_dataset = lgb.Dataset(X_val_df, label=y_val_df,\n                              categorical_feature=categorical_features)\n    model = lgb.train(params=params, train_set=train_dataset, num_boost_round=1000,\n                      valid_sets=[val_dataset],\n                      early_stopping_rounds=50, verbose_eval=25)\n    models.append(model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for model in models:\n    lgb.plot_importance(model)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.mean([model.best_score['valid_0']['rmse'] for model in models]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluating the model on the test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load the test dataframes first\ntest_df = reduce_mem_usage(load_df('test.csv'))\nweather_test_df = reduce_mem_usage(load_df('weather_test.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = test_df.assign(hour_of_day=test_df['timestamp'].dt.hour, \n                         day_of_week=test_df['timestamp'].dt.dayofweek,\n                         day_of_year=test_df['timestamp'].dt.dayofyear)\n\nmerged_test_df = (test_df\n                  .merge(building_metadata_df, how='left', on=['building_id'])\n                  .merge(weather_test_df, how='left', on=['site_id', 'timestamp']))\nX_test_df = merged_test_df[feature_cols]\nrow_ids = merged_test_df['row_id']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del test_df\ndel weather_test_df\ndel merged_test_df\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"describe_df(X_test_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results = sum(np.expm1(model.predict(X_test_df)) / len(models) \n              for model in models)\nresults_df = pd.DataFrame({'row_id': row_ids, \n                           'meter_reading': np.clip(results, 0, None)})\nresults_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results_df.to_csv('submission.csv', index=False, float_format='%.4f')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Notes\n\n* Why encode year_built? Isn't it already encoded?\n* Play around with imputing the empty features for the weather data\n* Play around with month of year vs. day of year\n* Consider converting site 0 readings to kWh just like the rest\n* Maybe floor count should be categorical?"}],"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":1}