{"cells":[{"metadata":{},"cell_type":"markdown","source":"\n<h1 align=center><font size = 4>ASHRAE - Great Energy Predictor III </font></h1>\n<h1 align=center><font size = 5>How much energy will a building consume?</font></h1>"},{"metadata":{},"cell_type":"markdown","source":"# Table of Contents\n* [Introduction/Business Problem](#introduction)\n* [Setup](#setup)\n* [Get the Data](#get_data)\n* [Take a Quick Look at the Data Structure](#data_structure)\n* [Explore the Data to Gain Insights](#explore)\n* [Prepare Data for ML](#preparation)\n* [Select and Train a Model](#selection)\n* [Make Predictions](#predictions)"},{"metadata":{},"cell_type":"markdown","source":"<a id = \"introduction\"></a>\n# Introduction/Business Problem"},{"metadata":{},"cell_type":"markdown","source":"Q: 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":"<a id=\"setup\"></a>\n# Setup"},{"metadata":{},"cell_type":"markdown","source":"Import a few common modules and ensure MatplotLib plots figures inline"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import seaborn as sns\nsns.set(style=\"darkgrid\")\n\n# to make this notebook's output stable across runs\nnp.random.seed(42)\n\n# To plot figures\n%matplotlib inline\nimport matplotlib\nimport matplotlib.pyplot as plt\nplt.rcParams['axes.labelsize'] = 14\nplt.rcParams['xtick.labelsize'] = 12\nplt.rcParams['ytick.labelsize'] = 12\n\nprint('Libraries imported.')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Memory optimization\n# Original code from https://www.kaggle.com/gemartin/load-data-reduce-memory-usage by @gemartin\n\ndef reduce_mem_usage(df):\n    \"\"\" iterate through all the columns of a dataframe and modify the data type\n        to reduce memory usage.        \n    \"\"\"\n    start_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))\n    \n    for col in df.columns:\n        col_type = df[col].dtype\n        \n        if col_type != object:\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                if 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        else:\n            df[col] = df[col].astype('object')\n\n    end_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))\n    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"get_data\"></a>\n# Get the Data"},{"metadata":{},"cell_type":"markdown","source":"Assessing the value of energy efficiency improvements can be challenging as there's no way to truly know how much energy a building would have used without the improvements. The best we can do is to build counterfactual models. Once a building is overhauled the new (lower) energy consumption is compared against modeled values for the original building to calculate the savings from the retrofit. More accurate models could support better market incentives and enable lower cost financing.\n\nThis competition challenges you to build these counterfactual models across four energy types based on historic usage rates and observed weather. The dataset includes three years of hourly meter readings from over one thousand buildings at several different sites around the world."},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_ashrae_energy_data(filename, ashrae_path = '../input/ashrae-energy-prediction/'):\n    csv_path = os.path.join(ashrae_path, filename)\n    return reduce_mem_usage(pd.read_csv(csv_path))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"data_structure\"></a>\n# Take a Quick Look at the Data Structure"},{"metadata":{},"cell_type":"markdown","source":"### Weather Dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train = load_ashrae_energy_data('weather_train.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# missing data\ntotal = weather_train.isnull().sum().sort_values(ascending=False)\npercent = (weather_train.isnull().sum()/weather_train.isnull().count()).sort_values(ascending=False)\nmissing_data = pd.concat([total, percent], axis=1, keys=['Total', 'Percent'])\nmissing_data.head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train.drop(['cloud_coverage', 'precip_depth_1_hr'], axis=1, inplace=True)\nweather_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Create lags for weather features**"},{"metadata":{"trusted":true},"cell_type":"code","source":"attributes = [\"air_temperature\", \"dew_temperature\", \"sea_level_pressure\", \"wind_direction\", \"wind_speed\"]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin\n\nclass LagWeatherFeatureCalculator(BaseEstimator, TransformerMixin):\n    def __init__(self, frequency='W', shift=1, attributes=['air_temperature']):\n        self.frequency = frequency\n        self.shift = shift\n        self.attributes = attributes\n    \n    def fit(self, X, y=None):\n        print('LagFeatureCalculator fit')\n        return self\n    \n    def transform(self, X, y=None):\n        print('LagFeatureCalculator transform')\n        print(\"Frequency is: {}\".format(self.frequency))\n        \n        X['timestamp'] = pd.to_datetime(X['timestamp'])\n        \n        frame = X.set_index(keys=['timestamp', 'site_id'])\n        frame_shifted = frame[self.attributes].unstack().resample(self.frequency).mean().shift(self.shift,freq=self.frequency).resample('H').ffill()\n        \n        columns_shifted = [col+'_'+ self.frequency +'%s' % 1 for col in frame.columns]\n        frame_shifted.columns.set_levels(columns_shifted, level=0, inplace=True)\n        \n        return frame.merge(frame_shifted.stack(), left_index=True, right_index=True, how='left').reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.pipeline import Pipeline\n\nlag_pipeline = Pipeline([\n    (\"hist_D1\", LagWeatherFeatureCalculator(frequency='D',shift=1, attributes=attributes)),\n    (\"hist_W1\", LagWeatherFeatureCalculator(frequency='W',shift=1, attributes=attributes)),\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lag_pipeline.fit_transform(weather_train).tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train_full = lag_pipeline.fit_transform(weather_train)\nweather_train_full.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Building Metadata"},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata = load_ashrae_energy_data('building_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# missing data\ntotal = building_metadata.isnull().sum().sort_values(ascending=False)\npercent = (building_metadata.isnull().sum()/building_metadata.isnull().count()).sort_values(ascending=False)\nmissing_data = pd.concat([total, percent], axis=1, keys=['Total', 'Percent'])\nmissing_data.head(20)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata.drop(['floor_count', 'year_built'], axis=1, inplace=True)\nbuilding_metadata.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Train Dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = load_ashrae_energy_data('train.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.info(verbose=True, null_counts=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['timestamp'] = pd.to_datetime(train['timestamp'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Merge Datasets"},{"metadata":{"trusted":true},"cell_type":"code","source":"ashrae = train.merge(building_metadata, on=['building_id'], how='left').merge(weather_train_full, on=['timestamp', 'site_id'], how='left')\nashrae.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ashrae.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del building_metadata, train, weather_train, weather_train_full","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Weird Values"},{"metadata":{},"cell_type":"markdown","source":"It is reported in this discussion by @barnwellguy that all electricity meter is 0 until May 20 for site_id == 0. I will remove these data from training data."},{"metadata":{"trusted":true},"cell_type":"code","source":"ashrae['meter_reading_log1p'] = np.log1p(ashrae['meter_reading'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_meter0_site0 = ashrae.query('meter == 0 & site_id ==0')\ntrain_data_meter0_site0.groupby('timestamp').sum()['meter_reading_log1p'].plot(figsize=(10, 5))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ashrae.drop(['meter_reading_log1p'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_meter0_site0.building_id.sort_values().unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(ashrae.query('site_id==0 & meter==0 & timestamp<=\"2016-05-20\"').shape)\nashrae = ashrae.query('not (site_id==0 & meter==0 & timestamp<=\"2016-05-20\")').reset_index(drop=True)\nashrae.info(verbose=True, null_counts=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Zero readings for electical meters: There's no reason for a building to ever have zero electrical usage, so I simply drop them all away."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(ashrae.query(\"(meter==0 & meter_reading==0)\").shape)\nashrae = ashrae.query(\"not (meter==0 & meter_reading==0)\").reset_index(drop=True)\nashrae.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Abnormally high readings from building 1099: These values are just absurdly high and don't fit. I remove them from the dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"ashrae.query(\"building_id==1099 and meter==2\").set_index(keys=['timestamp'])['meter_reading'].hist()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(ashrae.query(\"building_id==1099 & meter==2 & meter_reading>=3e4\").shape)\nashrae = ashrae.query(\"not (building_id==1099 & meter==2 & meter_reading>=3e4)\").reset_index(drop=True)\nashrae.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_data_meter0_site0, missing_data, total, percent","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id = \"explore\"></a>\n# Explore the Data to Get Insights"},{"metadata":{},"cell_type":"markdown","source":"Create a clean copy of the dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = ashrae.copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Target Variable"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data['meter_reading_log1p'] = np.log1p(train_data['meter_reading'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data[['meter_reading', 'meter_reading_log1p']].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_data[['meter_reading']])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#skewness and kurtosis\nprint(\"Skewness: %f\" % train_data['meter_reading'].skew())\nprint(\"Kurtosis: %f\" % train_data['meter_reading'].kurt())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_data[['meter_reading_log1p']])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#skewness and kurtosis\nprint(\"Skewness: %f\" % train_data[['meter_reading_log1p']].skew())\nprint(\"Kurtosis: %f\" % train_data[['meter_reading_log1p']].kurt())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Datetime Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"date_attributes = ['timestamp']\ntrain_data[date_attributes].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data['hour'] = train_data['timestamp'].dt.hour\ntrain_data['weekday'] = train_data['timestamp'].dt.weekday\ntrain_data['month'] = train_data['timestamp'].dt.month\ntrain_data[['timestamp','hour','weekday','month']].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"date_features = ['hour', 'weekday', 'month']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Categorical Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_features = ['meter', 'primary_use', 'site_id']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\n\nax1=plt.subplot(131)\nsns.countplot(x='meter', data=train_data[['meter']].replace(\n    {0:'electricity', 1:'chilledwater', 2:'steam', 3:'hotwater'}), ax=ax1)\nplt.xlabel('meter')\nplt.xticks(rotation=90)\n\nax2=plt.subplot(132)\nsns.countplot(x='primary_use', data=train_data, ax=ax2)\nplt.xlabel('primary_use')\nplt.xticks(rotation=90)\n\nax3=plt.subplot(133)\nsns.countplot(x='site_id', data=train_data, ax=ax3)\nplt.xlabel('site_id')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Numerical Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.impute import SimpleImputer\n\nfrom statsmodels.stats.outliers_influence import variance_inflation_factor\n\nclass ReduceVIF(BaseEstimator, TransformerMixin):\n    def __init__(self, thresh=5.0, impute=True, impute_strategy='median'):\n        # From looking at documentation, values between 5 and 15 are \"okay\".\n        # Above 10 is too high and so should be removed.\n        self.thresh = thresh\n        \n        # The statsmodel function will fail with NaN values, as such we have to impute them.\n        # By default we impute using the median value.\n        if impute:\n            self.imputer = SimpleImputer(strategy=impute_strategy)\n\n    def fit(self, X, y=None):\n        print('ReduceVIF fit')\n        print(self.imputer)\n        if hasattr(self, 'imputer'):\n            self.imputer.fit(X)\n        return self\n\n    def transform(self, X, y=None):\n        print('ReduceVIF transform')\n        columns = X.columns.tolist()\n        if hasattr(self, 'imputer'):\n            X = pd.DataFrame(self.imputer.transform(X), columns=columns)\n        return ReduceVIF.calculate_vif(X, self.thresh)\n\n    @staticmethod\n    def calculate_vif(X, thresh=5.0):\n        dropped=True\n        while dropped:\n            variables = X.columns\n            dropped = False\n            vif = [variance_inflation_factor(X[variables].values, X.columns.get_loc(var)) for var in X.columns]\n            \n            max_vif = max(vif)\n            if max_vif > thresh:\n                maxloc = vif.index(max_vif)\n                print(f'Dropping {X.columns[maxloc]} with vif={max_vif}')\n                X = X.drop([X.columns.tolist()[maxloc]], axis=1)\n                dropped=True\n        print(X.columns)\n        return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(24)\nm = 10000000\nidx = np.random.permutation(len(train_data))[:m]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#dropper = ReduceVIF(thresh=10)\n#num_features = list(dropper.fit_transform(train_data.drop(['meter_reading','meter_reading_log1p'], axis=1).select_dtypes(include=['int32','float32']).iloc[idx,:]).columns.values)\n#del dropper\n\nnum_features = ['square_feet', 'air_temperature', \n                'dew_temperature', 'wind_direction', 'wind_speed', 'dew_temperature_D1', 'wind_direction_D1', 'wind_speed_D1', 'dew_temperature_W1', 'wind_speed_W1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(24)\nm = 30000\nidx = np.random.permutation(len(train_data))[:m]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(data=train_data.loc[idx, num_features+['meter_reading','meter_reading_log1p']].dropna())\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for feature in (num_features+['meter_reading','meter_reading_log1p']):\n    print('Feature:{}\\tSkewness:{:.3f}\\tKurtosis:{:.3f}'.format(feature, train_data[feature].skew(), train_data[feature].kurt()))\n    \ntransf_features = [feature for feature in num_features if abs(train_data[feature].skew()) > 1 and abs(train_data[feature].kurt()) > 1]\ntransf_features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"non_transf_features = list(set(num_features).difference(set(transf_features)))\nnon_transf_features","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_transform = pd.concat([train_data[non_transf_features + ['meter_reading','meter_reading_log1p']],\n                                  train_data[transf_features].apply(lambda x: np.sign(x) * np.log(1 + np.abs(x)))], axis=1)\ntrain_data_transform","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for feature in num_features+['meter_reading','meter_reading_log1p']:\n    print('Feature:{}\\tSkewness:{:.3f}\\tKurtosis: {:.3f}'.format(feature, train_data_transform[feature].skew(), train_data_transform[feature].kurt()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data_transform.drop(['meter_reading'], axis=1, inplace=True)\ncorrmat = train_data_transform.corr()\nf, ax = plt.subplots(figsize=(7, 7))\nsns.heatmap(corrmat, vmax=.8, square=True);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"k = 10 #number of variables for heatmap\ncols = corrmat.nlargest(k, 'meter_reading_log1p')['meter_reading_log1p'].index\ncm = np.corrcoef(train_data_transform[cols].dropna().values.T)\nf, ax = plt.subplots(figsize=(7, 7))\nsns.set(font_scale=1.25)\nhm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_data, train_data_transform, cm, corrmat","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"preparation\"></a>\n# Prepare Data for ML"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin\n\n# A class to select numerical, categorical or datetime columns \nclass DataFrameSelector(BaseEstimator, TransformerMixin):\n    def __init__(self, attribute_names):\n        self.attribute_names = attribute_names\n        \n    def fit(self, X, y=None):\n        print(self.attribute_names)\n        return self\n    \n    def transform(self, X):\n        return X[self.attribute_names]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Handling Categorical Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_features = ['meter', 'primary_use', 'site_id', 'building_id']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class CategoricalImputer(BaseEstimator, TransformerMixin):\n    def fit(self, X, y=None):\n        self.most_frequent_ = pd.Series([X[attr].value_counts().index[0] for attr in X], index=X.columns)\n        return self\n    \n    def transform(self, X, y=None):\n        return X.fillna(self.most_frequent_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import OrdinalEncoder\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn.pipeline import Pipeline\n\ncat_pipeline = Pipeline([\n    (\"selector\", DataFrameSelector(cat_features)),\n    ('imputer', CategoricalImputer()),\n    #('encoder', OneHotEncoder(sparse=True))\n    (\"encoder\", OrdinalEncoder())\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Handling Datetime Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"date_attributes = ['timestamp']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"date_features = ['hour', 'weekday', 'month']\n\nclass TimeInfoExtractor(BaseEstimator, TransformerMixin):\n    def fit(self, X, y=None):\n        return self\n    \n    def transform(self, X, y=None):\n        return np.c_[X['timestamp'].dt.hour.astype(int),\n                     X['timestamp'].dt.weekday.astype(int),\n                     X['timestamp'].dt.month.astype(int)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"date_pipeline = Pipeline([\n    (\"selector\", DataFrameSelector(date_attributes)),\n    ('extractor', TimeInfoExtractor()),\n    #('encoder', OneHotEncoder(sparse=True))\n    ('encoder', OrdinalEncoder())\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Handling Numerical Attributes"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(num_features)\nprint(transf_features)\nprint(non_transf_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"transf_features_idx = [num_features.index(elem) for elem in transf_features]\nnon_transf_features_idx = [num_features.index(elem) for elem in non_transf_features]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(transf_features_idx)\nprint(non_transf_features_idx)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class LogModulusTransformer(BaseEstimator, TransformerMixin):\n    def fit(self, X, y=None):\n        return self\n    \n    def transform(self, X, y=None):\n        return np.c_[X[:,non_transf_features_idx], \n                     np.apply_along_axis(lambda x: np.sign(x) * np.log(1 + np.abs(x)), 1, X[:,transf_features_idx])]\n        # return X.apply(lambda x: np.sign(x) * np.log(1 + np.abs(x)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.pipeline import Pipeline\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.preprocessing import StandardScaler\n\nnum_pipeline = Pipeline([\n    ('selector', DataFrameSelector(num_features)),\n    ('imputer', SimpleImputer(strategy='mean')),\n    ('transformer', LogModulusTransformer()),\n    ('scaler', StandardScaler()),\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Transformation Pipelines"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = ashrae.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.pipeline import FeatureUnion\n\npreprocess_pipeline = FeatureUnion(transformer_list=[\n    (\"num_pipeline\", num_pipeline),\n    (\"cat_pipeline\", cat_pipeline),\n    (\"date_pipeline\", date_pipeline)\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = preprocess_pipeline.fit_transform(train_data)\ny = train_data['meter_reading'].apply(lambda x: np.log1p(x)).values\n\ndel  train_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(X.shape, y.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"total_features = num_features + cat_features + date_features\nprint(total_features)\n\ncategorical_features = cat_features + date_features\nprint(categorical_features)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"selection\"></a>\n# Select and Train a Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"np.random.seed(24)\nm = 5000000\nidx = np.random.permutation(len(X))[:m]\n\nX_subset = X[idx, :]\ny_subset = y[idx]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nX_train, X_valid, y_train, y_valid = train_test_split(X_subset, y_subset, test_size=0.2, random_state=42)\n\ndel X_subset, y_subset\n\nprint(X_train.shape, y_train.shape)\nprint(X_valid.shape, y_valid.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Use validation set for early stopping"},{"metadata":{"trusted":true},"cell_type":"code","source":"fit_params={\"early_stopping_rounds\": 5, \n            \"eval_metric\" : 'rmse', \n            \"eval_set\" : [(X_valid, y_valid)],\n            'eval_names': ['validation'],\n            'verbose': 100,\n            'categorical_feature': categorical_features, \n            'feature_name':total_features}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Set Hyperparameter Search"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats import randint, uniform, expon\n\nparams ={'num_leaves': randint(10, 50), \n         'min_child_samples': randint(500, 1000), \n         'colsample_bytree': uniform(loc=0.4, scale=0.6),\n         'learning_rate' : [0.001, 0.01, 0.1, 0.9, 1.5],\n         'subsample': uniform(loc=0.2, scale=0.8), \n         'colsample_bytree': uniform(loc=0.4, scale=0.6),\n         'reg_alpha': expon(scale=1.0)}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Randomized Search"},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nfrom sklearn.model_selection import RandomizedSearchCV\n\nlgb_reg = lgb.LGBMRegressor(max_depth=-1, random_state=42, silent=True, metric='None', n_jobs=4, n_estimators=1000)\n\nrandom_grd = RandomizedSearchCV(estimator=lgb_reg, param_distributions=params,\n                                n_iter=10, scoring='neg_mean_squared_error', cv=3, refit=True, random_state=42, verbose=2)\n\nrandom_grd.fit(X_train, y_train, **fit_params)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Best score reached: {} with params: {} '.format(random_grd.best_score_, random_grd.best_params_))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"optimal_params = random_grd.best_params_\n#optimal_params = {'colsample_bytree': 0.7905330837693118, 'learning_rate': 0.9, 'min_child_samples': 757, 'num_leaves': 33, 'reg_alpha': 1.786429543354675, 'subsample': 0.36987128854262097} ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cvres = random_grd.cv_results_\nfor mean_score, params in zip(cvres[\"mean_test_score\"], cvres[\"params\"]):\n    print(np.sqrt(-mean_score), params)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Build the Final Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Use the full dataset\nfrom sklearn.model_selection import train_test_split\nX_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, random_state=42)\n\ndel X, y\n\nprint(X_train.shape, y_train.shape)\nprint(X_valid.shape, y_valid.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lgb_opt = lgb.LGBMRegressor(max_depth=-1, random_state=42, silent=True, metric='None', n_jobs=4, n_estimators=1000)\n\n# set optimal parameters\nlgb_opt.set_params(**optimal_params)\n\nprint(lgb_opt)\n\nfit_params={\"early_stopping_rounds\": 10, \n            \"eval_metric\" : 'rmse', \n            \"eval_set\" : [(X_valid, y_valid)],\n            'eval_names': ['validation'],\n            'verbose': 100,\n            'categorical_feature': categorical_features, \n            'feature_name':total_features}\n\nt0, t1 = 900, 1000\ndef learning_schedule(t):\n    return t0 / (t + t1)\n\n#lgb_opt.fit(X_train, y_train, **fit_params, callbacks=[lgb.reset_parameter(learning_rate=learning_schedule)])\nlgb_opt.fit(X_train, y_train, **fit_params)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del X_train, y_train, X_valid, y_valid, ashrae","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Feature Importance"},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_importance = pd.DataFrame()\nfeature_importance[\"features\"] = total_features\nfeature_importance[\"importance\"] = lgb_opt.feature_importances_\n\n\nplt.figure(figsize=(10, 5))\nsns.barplot(x=\"importance\", y=\"features\", data=feature_importance.sort_values(by=\"importance\", ascending=False))\nplt.title(\"LightGBM Feature Importance\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"predictions\"></a>\n# Make Predictions"},{"metadata":{},"cell_type":"markdown","source":"### Weather Test Dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_test = load_ashrae_energy_data('weather_test.csv')\nweather_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_test.drop(['cloud_coverage', 'precip_depth_1_hr'], axis=1, inplace=True)\nweather_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_test_full = lag_pipeline.transform(weather_test)\nweather_test_full.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Building Metadata"},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata = load_ashrae_energy_data('building_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata.drop(['floor_count', 'year_built'], axis=1, inplace=True)\nbuilding_metadata.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test Dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"test = load_ashrae_energy_data('test.csv')\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test['timestamp'] = pd.to_datetime(test['timestamp'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Merge Datasets"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data = test.merge(building_metadata, on=['building_id'], how='left').merge(weather_test_full, on=['timestamp', 'site_id'], how='left')\ntest_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del building_metadata, test, weather_test, weather_test_full","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create Batches"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_instances = len(test_data)\nprint(n_instances)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch_size = 100000\nn_batches = n_instances // batch_size\nn_batches","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batches = np.array_split(test_data, n_batches)\ndel test_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred =[]\nfor n, batch in enumerate(batches):\n    if n % 50 == 0:\n        print(\"batch number: \", n)\n    y_pred.extend(np.expm1(lgb_opt.predict(preprocess_pipeline.transform(batch))))\n\ndel batches","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_pred = np.array(y_pred)\nprint(y_pred.shape)\ny_pred.ravel()\nprint(y_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame(y_pred).describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = load_ashrae_energy_data('sample_submission.csv')\nsample_submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = sample_submission.copy()\ndel sample_submission\n\nsubmission['meter_reading'] = np.clip(y_pred, 0, a_max=None)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.6"}},"nbformat":4,"nbformat_minor":1}