{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Overview\n\nCheck the page: https://www.kaggle.com/c/ashrae-energy-prediction/overview","metadata":{}},{"cell_type":"markdown","source":"# Dependency","metadata":{}},{"cell_type":"code","source":"#matplotlab inline\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utility functions","metadata":{}},{"cell_type":"code","source":"# Data loading functions\ndef loadTrainData():\n    building_metadata = pd.read_csv(\"../input/ashrae-energy-prediction/building_metadata.csv\")\n    train = pd.read_csv(\"../input/ashrae-energy-prediction/train.csv\")\n    weather_train = pd.read_csv(\"../input/ashrae-energy-prediction/weather_train.csv\")  \n    return building_metadata,train,weather_train\n\ndef loadTestData():\n    test = pd.read_csv(\"../input/ashrae-energy-prediction/test.csv\")\n    weather_test = pd.read_csv(\"../input/ashrae-energy-prediction/weather_test.csv\")\n    return test,weather_test\n\ndef summary(dfs):\n    for df in dfs:\n        print(df[1]+':\\n')\n        print(df[0].info(null_counts=True))\n        print('\\n\\n')\n        \ndef convertType(df): \n    for col in df.columns:\n        colType = str(df[col].dtypes)\n        if colType[:3]=='int':\n            if df[col].max()<=np.iinfo(np.int8).max:\n                df[col] = df[col].astype(np.int8)\n            elif df[col].max()<=np.iinfo(np.int16).max:\n                df[col] = df[col].astype(np.int16)\n            elif df[col].max()<=np.iinfo(np.int32).max:\n                df[col] = df[col].astype(np.int32)\n        elif colType[:5]=='float':\n            if df[col].max()<=np.finfo(np.float16).max:\n                df[col] = df[col].astype(np.float16)\n            elif df[col].max()<=np.finfo(np.float32).max:\n                df[col] = df[col].astype(np.float32)\n    return df\n\n\n\n# Visualize Trend (Moving Average) and Seasonality (Serial Correlation)\ndef plotTrendSeasonality(axes,df,building_id=0,meter=0,span=20,plotRaw=True):\n    if plotRaw:\n        s = df[(df['building_id']==building_id)&(df['meter']==meter)]['meter_reading']\n        axes[0].plot(s,label='Raw')\n\n    # Exponentially-weighted moving average (EWMA)\n    rolling = s.ewm(span=span)\n    axes[0].plot(rolling.mean(),label='EWMA')\n\n    # Autocorrelation\n    if not s.empty:\n        axes[1].acorr(s-rolling.mean(),usevlines=True, normed=True, maxlags=500)\n        \n        \n           \n# Detect the constant meter reading for select meter and building\ndef getSpanIdx(df,building_id,meter,maxLength=24):\n    tempData = df[(df['building_id']==building_id)&(df['meter']==meter)]['meter_reading']\n    x= tempData.values\n    idx = np.where(x[:-1]!=x[1:])[0]\n    idx = np.concatenate((idx,[len(x)-1]))\n    sp = np.concatenate(([idx[0]+1],np.diff(idx)))\n    if any(sp>maxLength):\n        return tempData.index[np.concatenate([range(x-y+1,x+1) for (x,y) in zip(idx[sp>maxLength],sp[sp>maxLength])])]\n    else:\n        return None\n\ndef plotConstantSpan(ax,df,building_id=0,meter=0,maxLength=24,plotRaw=True):\n    if plotRaw:\n        ax.plot(df[(df['building_id']==building_id)&(df['meter']==meter)]['meter_reading'],label='Raw')\n    idx = getSpanIdx(df=df,building_id=building_id,meter=meter,maxLength=maxLength)\n    if idx is not None:\n        ax.scatter(x=df_train.loc[idx].index,y=df_train.loc[idx]['meter_reading'],c='r',label='Constant Span')\n        \n        \n        \n# Visualize random building on select/random meter and site\ndef plotRandBuild(df_train,df_building_metadata,site_id=None, meter=None, maxLength=24):\n    if site_id is None:\n        rng = np.random.default_rng()\n        site_id = rng.integers(0,df_building_metadata.site_id.max()+1)\n    if meter is None:\n        rng = np.random.default_rng()\n        meter = rng.integers(0,4)\n    building_id = df_building_metadata[df_building_metadata['site_id']==site_id]['building_id'].sample(n=1).values[0]\n    \n    fig,axes = plt.subplots(1,2,figsize=[30,5])\n    axes[0].set_title('Site_id = {}, Building_id = {}, Meter = {}'.format(site_id,building_id,meter))\n    axes[1].set_title('Autocorrelation of residual')\n    axes[1].set_xlabel('Lags')\n    plotTrendSeasonality(axes,df=df_train,building_id=building_id,meter=meter,plotRaw=True)\n    plotConstantSpan(ax=axes[0],df=df_train,building_id=building_id,meter=meter,plotRaw=False)\n    axes[0].legend()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Import training data","metadata":{}},{"cell_type":"code","source":"df_building_metadata, df_train, df_weather_train = loadTrainData()\nsummary([(df_building_metadata,'building_metadata'), (df_train,'train'), (df_weather_train,'weather_train')])","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Reduce size for numeric columns.","metadata":{}},{"cell_type":"code","source":"df_building_metadata = convertType(df_building_metadata)\ndf_train = convertType(df_train) \ndf_weather_train = convertType(df_weather_train)\nsummary([(df_building_metadata,'building_metadata'), (df_train,'train'), (df_weather_train,'weather_train')])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Merge data","metadata":{}},{"cell_type":"code","source":"df_train_merge = df_train.merge(df_building_metadata, how='left',on='building_id').merge(df_weather_train,how='left',on=['site_id','timestamp'])\nsummary([(df_train_merge,'building_metadata')])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA\n\nWe summarize the comments at the beginning of each section and then show the stats and visualizations.","metadata":{}},{"cell_type":"markdown","source":"## Building data\n\n* Missing values in year_built and floor_count\n* Inf in year_built\n\n\n* No obvious linear dependency, except a bit between floor_count and square_feet\n\n\n* The top primary_use is from Education, as well as Office, Entertainment/public assembly, Lodging/residential, and Public services\n* Public services are mainly on site 3 and 8\n* Site 7 and 11 are only for Education","metadata":{}},{"cell_type":"code","source":"print(df_building_metadata.describe())\nsns.heatmap(df_building_metadata.corr(),vmin=-1,vmax=1,linewidths=.5,annot=True)\nsns.pairplot(df_building_metadata)\nsns.catplot(y='primary_use',kind='count',data=df_building_metadata.sort_values('primary_use'),color='c')\nsns.catplot(y='primary_use',kind='count',data=df_building_metadata.sort_values('primary_use'),color='c',col='site_id',height=5,aspect=.3)\nsns.catplot(x='site_id',y='square_feet',data=df_building_metadata,kind='bar')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Weather training data\n\nNote: timeseries visualizations per site are commented out as the figure takes too much memory\n\n* Missing values in all columns, except for site_id. cloud_coverage and precip_depth_1_hr contain a lot of missing values\n* Most of precip_depth_1_hr are 0\n* Timeseries data on site 15 are about 300 less than on other sites\n\n\n* aire_temp and dew_temp are highly correlated\n* cloud_coverage is of int type\n\n\n* Abnormal temp on site 13\n* Missing cloud coverage on site 7,11\n* Abnormal cloud coverage on site 1,5,6,9,10,14\n* Missing precip data on site 1,5,12 (or no precip at these locations)\n* Missing pressure data on site 5\n* Abnormal pressure behavior on site 1,12\n* A horizontal margin of wind speed appares on most sites\n* The pattern of wind speed on site 12 is different from others\n\n","metadata":{}},{"cell_type":"code","source":"print(df_weather_train.describe())\nprint('\\n')\nprint('Timeseries data count per site: \\n')\nprint(df_weather_train.site_id.value_counts(sort=False))\nsns.heatmap(df_weather_train.corr(),vmin=-1,vmax=1,linewidths=.5,annot=True)\nsns.pairplot(df_weather_train)\nsns.catplot(x='site_id',y='cloud_coverage',data=df_weather_train,kind='box')\ng = sns.FacetGrid(df_weather_train,subplot_kws=dict(projection='polar'),despine=False,height=10,sharex=False, sharey=False)\ng.map(sns.countplot,'wind_direction',alpha=.5)\ng = sns.FacetGrid(df_weather_train,col='site_id',subplot_kws=dict(projection='polar'),despine=False,col_wrap=4)\ng.map(sns.countplot,'wind_direction',alpha=.5)\n# sns.relplot(x='timestamp',y='air_temperature',data=df_weather_train,col='site_id')\n# sns.relplot(x='timestamp',y='precip_depth_1_hr',data=df_weather_train,col='site_id')\n# sns.relplot(x='timestamp',y='sea_level_pressure',data=df_weather_train,col='site_id')\n# sns.relplot(x='timestamp',y='wind_speed',data=df_weather_train,col='site_id')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Meter training data\n\nGiven each meter type (0:electricity, 1:chilledwater, 2:steam, 3:hotwater), we will explore the data at different levels as follows:\n* meter at all locations\n* meter on each site\n* meter in each building\n\nNote that we may not drill down all the analysis to the building level given that many buildings. Instead, we will only explore the meters in certain buildings to identify outliers and missing/ill-condition data. For example, the timeseries visualization per meter.\n\nKey results:\n* Meter 0 has the most observations, while the others have much less\n* No apparent correlation between overall meter_reading and other features\n* Meter 0 has low mean, max, and range\n* Meter 2 has the highest mean, deviation, max, and range\n* Meter 0 has some linear dependency to square_feet and floor_count\n* Meter 1 has small linear dependency to square_feet and year_build\n* Religious worship has extremely low consumption in meter 0. \n* Religious and Warehouse/storage have no reading for meter 1-3. Maybe such equipment/meter are not installed for some reason.\n* With respect to year_built, patterns from meters 0-2 somehow have similar behavior, especially after 1950. This may be driven by the facts on Site 3,8,15.\n* Site 7 has extremely high mean reading in meter 0 from Education\n* Site 13 has extremely high mean reading in meter 2 from Education\n* Site 0 has extremely high mean reading in meter 1 from Education, Lodging/residential and Other\n* Site 6 has extremely high mean reading in meter 1 from Entertainment/public_assembly\n* Site 15 has extremely high mean reading in meter 3 from Education\n* Site 7 has high mean reading in meter 3 from Education\n* Site 10 has high mean reading in meter 3 from Entertainment/public_assembly\n* The overall daily mean meter_reading for meter 1 grows significantly in Q3, which is mainly driven by site 6\n* The overall daily mean meter_reading for meter 2 is extremely high in H1 and has a spike in Q4, which is mainly driven by site 13\n* Meter 0 on site 0 seems not working until mid of Q2\n* Meter 0 on sites 8,9, and 14 may contain outliers for certain days\n* Meter 0,2, and 3 on site 15 have consecutive missing values during the same time period in Q1\n* Meter 1 on site 13 may contain outliers in Q3 and Q4\n* Meter 3 on sites 1 may contain outliers for certain days\n* Meters from some buildings have conseutive reading for more than a day which is abnormal\n\n\n","metadata":{}},{"cell_type":"code","source":"# Overall for all meters\nprint(df_train_merge.describe())\nprint('\\n')\nprint('Timeseries data count per meter: \\n')\nprint(df_train_merge.meter.value_counts(sort=False))\nsns.heatmap(df_train_merge.corr(),vmin=-1,vmax=1,linewidths=.5,annot=True)\n\n# Overall for each meter\nsns.catplot(x='meter',y='meter_reading',data=df_train_merge,kind='box',showfliers=False)\n\nplt.figure(figsize=(50,10))\nfor i in range(4):\n    plt.subplot(1,4,i+1)\n    sns.heatmap(df_train_merge[df_train_merge.meter==i].corr(),vmin=-1,vmax=1,linewidths=.5,annot=True)\n    \nsns.catplot(x='meter_reading',y='primary_use',data=df_train_merge,col='meter',kind='bar',sharex=False)\nsns.relplot(x='year_built',y='meter_reading',data=df_train_merge,hue='meter',kind='line',ci=None,aspect=5)\nsns.catplot(y='meter_reading',x='site_id',data=df_train_merge,col='meter',kind='bar',sharey=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Each meter per site\nsns.catplot(x='meter_reading',y='primary_use',data=df_train_merge,col='meter',row='site_id',kind='bar',sharex=False)\nsns.relplot(x='year_built',y='meter_reading',data=df_train_merge,hue='meter',kind='line',ci=None,aspect=5,row='site_id',facet_kws=dict(sharey=False))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Timeseries for each meter\ndf_train_merge['timestamp'] = pd.to_datetime(df_train_merge['timestamp'])\n\nplt.figure(figsize=[30,5])\nplt.subplot(131)\nsns.lineplot(y='meter_reading',x='timestamp', data=df_train_merge.groupby([df_train_merge['timestamp'].dt.month,'meter']).mean().reset_index(),ci=None,hue='meter').set_xlabel('Month')\nplt.subplot(132)\nsns.lineplot(y='meter_reading',x='timestamp', data=df_train_merge.groupby([df_train_merge['timestamp'].dt.week,'meter']).mean().reset_index(),ci=None,hue='meter').set_xlabel('Week')\nplt.subplot(133)\nsns.lineplot(y='meter_reading',x='timestamp', data=df_train_merge.groupby([df_train_merge['timestamp'].dt.dayofyear,'meter']).mean().reset_index(),ci=None,hue='meter').set_xlabel('Day of Year')\n\nsns.relplot(y='meter_reading',x='timestamp', data=df_train_merge.groupby([df_train_merge['timestamp'].dt.dayofyear,'meter']).mean().reset_index(),\n            ci=None,row='meter',kind='line',facet_kws=dict(sharey=False),aspect=4)\nsns.relplot(y='meter_reading',x='timestamp', data=df_train_merge.groupby([df_train_merge['timestamp'].dt.dayofyear,'meter','site_id']).mean().reset_index(),\n            ci=None,col='meter',row='site_id',kind='line',facet_kws=dict(sharey=False),aspect=4)\n\nsns.relplot(x='timestamp',y='meter_reading',data=df_train_merge.set_index(['site_id','meter']).loc[(9,0)],kind='line',ci=None,aspect=4)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Trend, seasonality and constant span for random meters\nplotRandBuild(df_train,df_building_metadata)\nplotRandBuild(df_train,df_building_metadata,meter=0)\nplotRandBuild(df_train,df_building_metadata,site_id=6,meter=1)\nplotRandBuild(df_train,df_building_metadata,site_id=13,meter=2)\nplotRandBuild(df_train,df_building_metadata,site_id=15)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}