{"cells":[{"metadata":{},"cell_type":"markdown","source":"## ASHRAE - Great Energy Predictor III\n\nThis Notebook is divided into couple of sections which we are going to discuss futher. This will give us an edge for understanding the dataset better and perform the feature engineering later\n\n<b>Section I. Import dataset\n\nSection II. Descriptive Stats\n\nSection III. Train Test Comparison\n\nSection IV. Normalization and Merge\n\nSection V. Exploration of Null values\n\nSection VI. Univariate Distribution\n\nSection VII. Bivariate Distribution\n\nSection VIII.Understanding seasonality\n\nSection IX. Feature Normalization\n\nSection X. Model Fitting\n\nSection XI. Model Evaluation & Feature Importance\n\n</b>"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","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# Let's see the discribution of target in train data\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nfrom plotly.offline import init_notebook_mode, plot, iplot\ninit_notebook_mode(connected = True)\nimport plotly.graph_objs as go\nimport warnings\nwarnings.filterwarnings('ignore')\nplt.style.use('ggplot')\n\n%matplotlib inline\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.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Give an overall definition of the dataframe\ndef define_df(df):\n    tup = df.shape\n    print(\"Number of rows in dataframe: {0} & cols: {1}\".format(tup[0], tup[1]))\n    # Get number of null values in dataframe\n    print(\"-----------Number of Null values-----------\")\n    print(df.isna().sum())\n    #Descriptive stats\n    print(\"-----------Descriptive Stats-----------\")\n    print(df.describe())\n    #Total memory consumption\n    print(\"-----------Data Types + Memory Consumption-----------\")\n    print(df.info())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"define_df(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building_metadata = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/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":"define_df(building_metadata)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/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":"define_df(weather_train)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's use the usual code to minimize the size of dataframe\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        # Custom implementation for this dataset\n        if col_type != object and col_type != np.datetime64 and col != 'timestamp':\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int': #Encode with the most relevant datatype.g\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\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":{"trusted":true},"cell_type":"code","source":"from time import time","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time\ntrain = reduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time\nbuilding_metadata = reduce_mem_usage(building_metadata)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time\nweather_train = reduce_mem_usage(weather_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For a particular train row to match with its building metadata we have to understand the nuances of weather conditions at a particular time. It's entirely possible that due to some extreme weather the meters stopped working which can be captured by timestamp in weather date.\n\nProcessing date -- Get day, week, month and year of the energy consumption. This will help us in understanding if there is any seasonality within the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"def process_date(df):\n    # Thanks to the solution provided at https://stackoverflow.com/questions/25146121/extracting-just-month-and-year-separately-from-pandas-datetime-column\n    # This made my code way shorter\n    lst = ['month', 'day', 'hour', 'dayofweek']\n    df['timestamp'] = pd.to_datetime(df.timestamp)\n    df = df.join(pd.concat((getattr(df.timestamp.dt, i).rename(i) for i in lst), axis=1))\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Now merge the data\ntemp = pd.merge(building_metadata, weather_train, on = 'site_id', how = 'inner')\ntrain = pd.merge(train, temp, on=['building_id', 'timestamp'], how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# reduce memory usage\n%time\ntrain = reduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del weather_train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import gc\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = process_date(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# reduce memory usage\n%time\ntrain = reduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The 99.1 percentile of target i.e. meter_reading is at 5449. Thus, the target follows a pareto distribution. I will plot 2 distribution of meter reading which gives us an idea of how it looks like."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(15,5))\nsns.distplot(train.meter_reading, ax=ax[0])\nsns.distplot(train[train.meter_reading < 5450].meter_reading, ax=ax[1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are zeros in the target which have multiple theories. Let's try to understand how many buildings in each site id have monthly meter reading as ~0."},{"metadata":{"trusted":true},"cell_type":"code","source":"def process_mr(reading):\n    if reading < 1:\n        return 1\n    elif reading >= 1 and reading < 10:\n        return 2\n    elif reading >= 10 and reading < 100:\n        return 3\n    elif reading >= 100 and reading < 1000:\n        return 4\n    elif reading >= 1000 and reading < 5450:\n        return 5\n    else:\n        return 6\ntemp = train.groupby(['site_id', 'building_id', 'month']).agg({'meter_reading':np.mean}).reset_index()\ntemp['groups'] = temp.meter_reading.apply(lambda x: process_mr(x))\ntemp = temp.groupby(['site_id', 'month', 'groups']).agg({'building_id': 'count'}).reset_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Thanks to for some of the great plot at this kernel [https://www.kaggle.com/nroman/eda-for-ashrae]. I have used reused some of the codebase over in my kernel too."},{"metadata":{"trusted":true},"cell_type":"code","source":"for site in temp.site_id.unique():\n    t = temp[temp.site_id == site]\n    fig, ax = plt.subplots(1,2, figsize=(16,4))\n    sns.heatmap(data=t.pivot_table(index='groups', columns='month', values='building_id'), ax=ax[0], cmap=\"YlGnBu\")\n    ax[0].set_title('Month vs Group - site_id: {}'.format(site))\n    train[train['site_id'] == site][['timestamp', 'meter_reading']].set_index('timestamp').resample('H').mean()['meter_reading'].plot(ax=ax[1], alpha=0.8, label='By hour', color='tab:blue').set_ylabel('Mean meter reading', fontsize=13);\n    train[train['site_id'] == site][['timestamp', 'meter_reading']].set_index('timestamp').resample('D').mean()['meter_reading'].plot(ax=ax[1], alpha=1, label='By day', color='tab:orange').set_xlabel('');\n    ax[1].set_title('Train Month vs Meter reading - site_id: {}'.format(site))\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is quite fascinating with some useful insights. \n1. Group 6 i.e. with meter reading > 5450  is present in only few of the site id namely: 1,6,10,13,14.\n2. There is been a lot of talk recently of how site id 0 having meter reading ~0 for first few months for most of the buildings and this is evident from the visualization where almost 90% of buildings have < 1 meter reading.\n3. The groups are not equally distributed which means all site ids don't have equal distribution for the meter reading.\n\nI am going to ignore the readings for point site id 0 for first 4 months. This may be a instrumental error which we have to look more closely into."},{"metadata":{},"cell_type":"markdown","source":"#### Let's understand the types of meter and their respective reading. As defined in data description,we have 4 types: {0: electricity, 1: chilledwater, 2: steam, 3: hotwater}. Let's plot an overall distribution for it."},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/test.csv\")\nweather_test = pd.read_csv(\"/kaggle/input/ashrae-energy-prediction/weather_test.csv\")\n%time\ntest =reduce_mem_usage(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.drop(columns=['row_id'], inplace=True)\n# Now merge the data\ntemp = pd.merge(building_metadata, weather_test, on = 'site_id', how = 'inner')\ntest = pd.merge(test, temp, on=['building_id', 'timestamp'], how='left')\n# reduce memory usage\n%time\ntest = reduce_mem_usage(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test = process_date(test)\ntest.drop(columns=['timestamp'], inplace=True)\n# reduce memory usage\n%time\ntest = reduce_mem_usage(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(train.meter.value_counts()).reset_index()\ntemp2 = pd.DataFrame(test.meter.value_counts()).reset_index()\ntemp.columns = ['meter', 'Counts']\ntemp2.columns = ['meter', 'Counts']\ntrace1 = go.Bar(x = temp.meter, y = temp.Counts, name='Train')\ntrace2 = go.Bar(x = temp2.meter, y = temp2.Counts, name='Test')\ndata = [trace1, trace2]\nlayout = go.Layout(\n    title = \"Meter Type Distribution\",\n    xaxis = dict(\n        title = \"Meter Type\",\n        tickfont = dict(\n        color = 'rgb(107,107,107)'\n        )\n    ),\n    yaxis = dict(\n        title = \"Counts\",\n        titlefont = dict(\n        color = 'rgb(107,107,107)'\n        ),\n        tickfont=dict(\n        color = 'rgb(107,107,107)'\n        )\n    )\n)\nfig = go.Figure(data = data, layout = layout)\niplot(fig,image_height=12,image_width=15, filename = 'Meter Type Dist')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(train.primary_use.value_counts()).reset_index()\ntemp2 = pd.DataFrame(test.primary_use.value_counts()).reset_index()\ntemp.columns = ['primary_use', 'Counts']\ntemp2.columns = ['primary_use', 'Counts']\ntrace1 = go.Bar(x = temp.primary_use, y = temp.Counts, name='Train')\ntrace2 = go.Bar(x = temp2.primary_use, y = temp2.Counts, name='Test')\ndata = [trace1, trace2]\nlayout = go.Layout(\n    title = \"Primary Use Distribution\",\n    xaxis = dict(\n        title = \"Primary Use\",\n        tickfont = dict(\n        color = 'rgb(107,107,107)'\n        )\n    ),\n    yaxis = dict(\n        title = \"Counts\",\n        titlefont = dict(\n        color = 'rgb(107,107,107)'\n        ),\n        tickfont=dict(\n        color = 'rgb(107,107,107)'\n        )\n    )\n)\nfig = go.Figure(data = data, layout = layout)\niplot(fig,image_height=12,image_width=15, filename = 'Primary Use Dist')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(train.year_built.value_counts()).reset_index()\ntemp2 = pd.DataFrame(test.year_built.value_counts()).reset_index()\ntemp.columns = ['year_built', 'Counts']\ntemp2.columns = ['year_built', 'Counts']\ntrace1 = go.Bar(x = temp.year_built, y = temp.Counts, name='Train')\ntrace2 = go.Bar(x = temp2.year_built, y = temp2.Counts, name='Test')\ndata = [trace1, trace2]\nlayout = go.Layout(\n    title = \"Year Built Distribution\",\n    xaxis = dict(\n        title = \"Year Built\",\n        tickfont = dict(\n        color = 'rgb(107,107,107)'\n        )\n    ),\n    yaxis = dict(\n        title = \"Counts\",\n        titlefont = dict(\n        color = 'rgb(107,107,107)'\n        ),\n        tickfont=dict(\n        color = 'rgb(107,107,107)'\n        )\n    )\n)\nfig = go.Figure(data = data, layout = layout)\niplot(fig,image_height=12,image_width=15, filename = 'Year Built Dist')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(train.floor_count.value_counts()).reset_index()\ntemp2 = pd.DataFrame(test.floor_count.value_counts()).reset_index()\ntemp.columns = ['floor_count', 'Counts']\ntemp2.columns = ['floor_count', 'Counts']\ntrace1 = go.Bar(x = temp.floor_count, y = temp.Counts, name='Train')\ntrace2 = go.Bar(x = temp2.floor_count, y = temp2.Counts, name='Test')\ndata = [trace1, trace2]\nlayout = go.Layout(\n    title = \"floor_count Distribution\",\n    xaxis = dict(\n        title = \"floor_count\",\n        tickfont = dict(\n        color = 'rgb(107,107,107)'\n        )\n    ),\n    yaxis = dict(\n        title = \"Counts\",\n        titlefont = dict(\n        color = 'rgb(107,107,107)'\n        ),\n        tickfont=dict(\n        color = 'rgb(107,107,107)'\n        )\n    )\n)\nfig = go.Figure(data = data, layout = layout)\niplot(fig,image_height=12,image_width=15, filename = 'Floor Count Dist')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 'year_built', 'floor_count'\nfor col in ['square_feet', 'air_temperature',\n       'cloud_coverage', 'dew_temperature', 'precip_depth_1_hr',\n       'sea_level_pressure', 'wind_direction', 'wind_speed']:\n    plt.figure(figsize=(10,5))\n    sns.distplot(train[~train[col].isna()][col],  kde_kws={\"lw\": 3, \"label\": 'Train'})\n    sns.distplot(test[~test[col].isna()][col],  kde_kws={\"lw\": 3, \"label\": 'Test'})\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Electricity meter type has the maximum readings. Also, this follows an exponentially decreasing relationship with 0>1>2>3\n\nLet's try to understand their relationship with the target. Once again I am going to plot a general distribution of target but w.r.t Meter Type"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(4,2, figsize=(15, 20))\nmeter_dict = {0: \"electricity\", 1: \"chilledwater\", 2: \"steam\", 3: \"hotwater\"}\nmeter = [0,1,2,3]\nfor i, label in enumerate(meter):\n    sns.distplot(train[train.meter == label].meter_reading, ax=ax[int(i)][0], kde_kws={\"lw\": 3, \"label\": meter_dict[label]})\n    sns.distplot(train[(train.meter == label) & (train.meter_reading < 5450)].meter_reading, ax=ax[int(i)][1], kde_kws={\"lw\": 3, \"label\": meter_dict[label]})\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This gives us some better idea about target distribution. Most of the zeros are coming in steam meter type while smallest one is for hot water. "},{"metadata":{},"cell_type":"markdown","source":"Identify if there is seasonality according to week or month, \n\nAlso make a time plot with days distribution along with hours values in it. Perform the same for month and days along with year and month"},{"metadata":{"trusted":true},"cell_type":"code","source":"def seasonality(df, col):\n    Z = df.groupby(col).agg({'meter_reading': sum}).reset_index()\n    fig, ax = plt.subplots(ncols=1, sharey=True, figsize = (10,5))\n    sns.barplot(data=Z, x=col, y='meter_reading', ax = ax, palette=\"BrBG\")\n    plt.title('Overall Meter reading over the course of an year')\n    plt.show()\n    \n    # Next plot if according to different meter types\n    Z = df.groupby([col, 'meter']).agg({'meter_reading': sum}).reset_index()\n    \n    fig, ax = plt.subplots(ncols=4, sharey=True, figsize = (20,5))\n    sns.barplot(data=Z[Z.meter == 0], x=col, y='meter_reading', ax = ax[0], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 1], x=col, y='meter_reading', ax = ax[1], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 2], x=col, y='meter_reading', ax = ax[2], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 3], x=col, y='meter_reading', ax = ax[3], palette=\"BrBG\")\n    plt.show()\n    \n    # Let's truncate the values to smaller values\n    Z = df[df.meter_reading<5450].groupby(col).agg({'meter_reading': sum}).reset_index()\n    fig, ax = plt.subplots(ncols=1, sharey=True, figsize = (10,5))\n    sns.barplot(data=Z, x=col, y='meter_reading', ax = ax, palette=\"BrBG\")\n    plt.title('Overall Meter reading over the course of an year')\n    plt.show()\n    \n    # Next plot if according to different meter types\n    Z = df[df.meter_reading<5450].groupby([col, 'meter']).agg({'meter_reading': sum}).reset_index()\n    \n    fig, ax = plt.subplots(ncols=4, sharey=True, figsize = (20,5))\n    sns.barplot(data=Z[Z.meter == 0], x=col, y='meter_reading', ax = ax[0], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 1], x=col, y='meter_reading', ax = ax[1], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 2], x=col, y='meter_reading', ax = ax[2], palette=\"BrBG\")\n    sns.barplot(data=Z[Z.meter == 3], x=col, y='meter_reading', ax = ax[3], palette=\"BrBG\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seasonality(train,'month')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seasonality(train, 'day')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seasonality(train, 'hour')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seasonality(train, 'dayofweek')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So far we have seen the relationship of one time variable w.r.t. target. Let's spice things up and see the distribution of 2 time variables w.r.t. target"},{"metadata":{"trusted":true},"cell_type":"code","source":"month_day_agg_sum = train[['month', 'day', 'meter_reading']].groupby(['month', 'day']).meter_reading.sum().reset_index()\nmonth_day_agg_sum_limited = train[['month', 'day', 'meter_reading']][train.meter_reading<5450].groupby(['month', 'day']).meter_reading.sum().reset_index()\nfig, ax = plt.subplots(1,2, figsize=(16,8))\nsns.heatmap(data=month_day_agg_sum.pivot_table(index='day', columns='month', values='meter_reading'), ax=ax[0], cmap=\"YlGnBu\", label='Month vs day')\nsns.heatmap(data=month_day_agg_sum_limited.pivot_table(index='day', columns='month', values='meter_reading'), ax=ax[1], cmap=\"YlGnBu\", label='Month vs dat--limited')\nax[0].set_title('Month vs Day')\nax[1].set_title('Month vs Day--limited')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"month_day_agg_sum = train[['day', 'hour', 'meter_reading']].groupby(['day', 'hour']).meter_reading.sum().reset_index()\nmonth_day_agg_sum_limited = train[['day', 'hour', 'meter_reading']][train.meter_reading<5450].groupby(['day', 'hour']).meter_reading.sum().reset_index()\nfig, ax = plt.subplots(2,1, figsize=(16,12))\nsns.heatmap(data=month_day_agg_sum.pivot_table(index='hour', columns='day', values='meter_reading'), ax=ax[0], cmap=\"YlGnBu\")\nsns.heatmap(data=month_day_agg_sum_limited.pivot_table(index='hour', columns='day', values='meter_reading'), ax=ax[1], cmap=\"YlGnBu\")\nax[0].set_title('Day vs hour')\nax[1].set_title('Day vs hour--limited')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Few interesting things that we captured:\n    1. The meter reading is in general higher for steam (approximately 10 times). If we restrict the values to smaller extent then steam type changes its distribution. The overall distribution is then dominated by electricity type.\n    2. If we try to locate the highest consumption days in month then it usually starts from 6th and ends in 12th month with approximately even distribution over the days. \n    3. The day distribution makes sense (for limited values) since most of the readings are in working hours compared to non working hours. \n    4. If we try to understand the seasonality according to meter types per month then, it makes sense that chilled water is mostly used during the time of summer while steam and hot water is mostly used in the time of winter.\n    "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Now delete the test related data in order to understand feature importance otherwise the kernels will be overloaded.\ndel test, weather_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.corr()['meter_reading']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### The weather condition for the time at which reading was taken won't be helpful. We have to consider other approaches. For an instance, we have to take mean, max, min values of different weather conditions over the course of month [for each site], week of year [for each site]. This will help out model in understanding the weather condition much better which helps in explaining the consumption of energy.\n#### Remove site id 0 for first few month since they are anomalous events."},{"metadata":{"trusted":true},"cell_type":"code","source":"# This function helps in aggregating the variables over different statistical functions\ndef agg_numeric(df, group_var = [], to_group = [], df_name = 'dummy'):\n    \n    if len(to_group) == 0 or len(group_var) == 0:\n        raise ValueError('Please check grouping and to be grouped variables again!!!')\n        \n    cols = group_var + to_group\n    \n    agg = df[cols].groupby(group_var).agg([np.min, np.mean, np.max]).reset_index()\n    \n    columns = group_var \n    \n    for var in agg.columns.levels[0]:\n        if var not in group_var:\n            for stat in agg.columns.levels[1][:-1]:\n                 columns.append('%s_%s_%s' % (df_name, var, stat))\n                    \n    agg.columns = columns\n    return agg","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Removing site id 0 for first few months\ntrain = train[~((train.site_id == 0) & (train.month <= 5))]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = agg_numeric(train, group_var=['site_id', 'month'], \n                   to_group = ['air_temperature', 'cloud_coverage', 'wind_speed'],\n                   df_name = 'weather')\ntrain = pd.merge(train, temp, on = ['site_id', 'month'], how = 'inner')\nreduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Another feature engineering -- Age of building\ntrain['age'] = 2016 - train['year_built']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.corr()['meter_reading']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def weekend(dayofweek):\n    if dayofweek in [5, 6]:\n        return 1\n    else:\n        return 0\ntemp['dayofweek'] = train.dayofweek.apply(lambda x: weekend(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.drop(columns = ['precip_depth_1_hr', 'timestamp', 'year_built', 'wind_direction', 'day', 'hour'], inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# building features over every month\ntemp = agg_numeric(train, group_var=['building_id', 'month'], \n                   to_group = ['meter_reading'],\n                   df_name = 'building')\ntrain = pd.merge(train, temp, on = ['building_id', 'month'], how = 'inner')\nreduce_mem_usage(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.corr()['meter_reading']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Analyse the usefulness of variables using feature importance methodology."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Taken from XGBoost documentation\n\ndef gradient(predt, dtrain) -> np.ndarray:\n    '''Compute the gradient squared log error.'''\n    y = dtrain.get_label()\n    return (np.log1p(predt) - np.log1p(y)) / (predt + 1)\n\ndef hessian(predt, dtrain) -> np.ndarray:\n    '''Compute the hessian for squared log error.'''\n    y = dtrain.get_label()\n    return ((-np.log1p(predt) + np.log1p(y) + 1) /\n            np.power(predt + 1, 2))\n\ndef squared_log(predt,\n                dtrain) -> tuple([np.ndarray, np.ndarray]):\n    '''Squared Log Error objective. A simplified version for RMSLE used as\n    objective function.\n    '''\n    grad = gradient(predt, dtrain)\n    hess = hessian(predt, dtrain)\n    return grad, hess","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import math\ndef rmsle(predt, dtrain):\n    ''' Root mean squared log error metric.'''\n    y = dtrain.get_label()\n    elements = np.power(np.log1p(y) - np.log1p(predt), 2)\n    return 'PyRMSLE', float(np.sqrt(np.sum(elements) / len(y))), False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import lightgbm as lgb\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold\n # lightbgm\nparams = {\n     'max_depth' : 8,\n     'boosting_type': 'gbdt',\n     'objective': 'regression',\n     'metric': {'rmse'},\n     'subsample': 0.2,\n     'learning_rate': 0.1,\n     'feature_fraction': 0.8,\n     'bagging_fraction': 0.9,\n     'alpha': 0.1,\n     'lambda': 0.1\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = train.columns.tolist()\ncols.remove(\"meter_reading\"); #cols.remove(\"year_built\")\ntrain_label = np.log1p(train.meter_reading)\ntrain = train[cols]\n# Create a label encoder -- for object type\nlabel_encoder = LabelEncoder()\n\nfor i, col in enumerate(cols):\n    if train[col].dtype == 'object':\n        # Map the categorical features to integer\n        train[col] = label_encoder.fit_transform(np.array(train[col].astype(str)).reshape((-1,)))\n\n# Define categorical cols as well\ncat_cols = [0,1,2,3,11]\nfolds = 3\nseed = 666\nkf = KFold(n_splits=folds, shuffle=True, random_state=seed)\nfeat_importance=None\ncount = 0\nfor train_index, val_index in kf.split(train):\n    train_X = train.iloc[train_index]\n    val_X = train.iloc[val_index]\n    train_y = train_label.iloc[train_index]\n    val_y = train_label.iloc[val_index]\n    lgb_train = lgb.Dataset(train_X, train_y, categorical_feature=cat_cols)\n    lgb_eval = lgb.Dataset(val_X, val_y, categorical_feature=cat_cols)\n    gbm = lgb.train(params,\n                     lgb_train,\n                     num_boost_round=50, #300,\n                     valid_sets=(lgb_train, lgb_eval),\n                     early_stopping_rounds= 5,#100,\n                    fobj = squared_log, \n                    feval = rmsle,\n                     verbose_eval=2) #100)\n    if count == 0:\n        feat_importance = pd.DataFrame(zip(gbm.feature_importance(), gbm.feature_name()), columns=['Value','Feature'])\n    else:\n        feat_importance2 = pd.DataFrame(zip(gbm.feature_importance(), gbm.feature_name()), columns=['Value{0}'.format(count),'Feature'])\n        feat_importance=pd.merge(feat_importance, feat_importance2, on='Feature', how='inner')\n    count+=1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"feat_importance['Value'] = (feat_importance.Value+feat_importance.Value1+feat_importance.Value2) / 3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfi_df =feat_importance.sort_values('Value', ascending=False)\n\ntrace = go.Bar(x = fi_df.Feature, y = fi_df.Value)\n\ndata = [trace]\nlayout = go.Layout(\n    title = \"Feature importance of LGBM Model\",\n    xaxis = dict(\n        title = \"Columns\",\n        tickfont = dict(\n        color = 'rgb(107,107,107)'\n        )\n    ),\n    yaxis = dict(\n        title = \"Feature Importance\",\n        titlefont = dict(\n        color = 'rgb(107,107,107)'\n        ),\n        tickfont=dict(\n        color = 'rgb(107,107,107)'\n        )\n    )\n)\nfig = go.Figure(data = data, layout = layout)\niplot(fig,image_height=12,image_width=15, filename = 'feature importance status')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#### Next Stop:\n1. Lag variables for the same building will also help in identifying the potential energy consumption but that's not possible for target since we are doing projection for next 2 years. We can do lags on weather conditions though [This lag has to be applied by combining both train and test]. \n2. Holiday/weekend based features. -- Done\n2. Building history in terms of energy consumption -- Done\n3. Remove highly correlated variables or apply dimensionality reduction.\n5. Implement customized function for RMSLE and use it within the algorithm -- Done.\n6. Follow the corochann kernel to make a model fit according to different meter type since energy consumption is on different scales for them."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}