{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd \nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport gc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#standard data usage reduction, making sure filetype isn't unnescessarily large\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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/train.csv'))\nbuildings = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/building_metadata.csv'))\nweather_train = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/weather_train.csv'))\ntest_df = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/test.csv'))\nweather_test = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/weather_test.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#convert 'timestamp' to datetime format\ntrain_df[\"timestamp\"] = pd.to_datetime(train_df[\"timestamp\"], format='%Y-%m-%d %H:%M:%S')\nweather_train[\"timestamp\"] = pd.to_datetime(weather_train[\"timestamp\"], format='%Y-%m-%d %H:%M:%S' )\ntest_df[\"timestamp\"] = pd.to_datetime(test_df[\"timestamp\"], format='%Y-%m-%d %H:%M:%S' )\nweather_test[\"timestamp\"] = pd.to_datetime(weather_test[\"timestamp\"], format='%Y-%m-%d %H:%M:%S' )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"merged_train_df = pd.merge(train_df,buildings,on='building_id',how='left')\nmerged_train_df = pd.merge(merged_train_df,weather_train,on=['site_id','timestamp'],how='left')\nmerged_test_df = pd.merge(test_df,buildings,on='building_id',how='left')\nmerged_test_df = pd.merge(merged_test_df,weather_test,on=['site_id','timestamp'],how='left')\ndel train_df, test_df, weather_train,weather_test, buildings\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = merged_train_df\ntest_df = merged_test_df\ndel merged_train_df,merged_test_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We now have all of our data in train_df and test_df"},{"metadata":{},"cell_type":"markdown","source":"# EDA"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A look at target variable"},{"metadata":{},"cell_type":"markdown","source":"### Target variable distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['meter_reading'].describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"10^5 gap between 3rd quartile and maximum, data is heaviliy skewed"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nplt.subplot(211)\nplt.hist(train_df['meter_reading'],30)\nplt.title('distribution of meter_reading in train data')\nplt.xlabel('value of meter reading')\nplt.ylabel('count')\nplt.subplot(212)\nplt.hist(train_df['meter_reading'],30)\nplt.title('logscale distribution of meter_reading in train data')\nplt.xlabel('value of meter reading')\nplt.ylabel('log(count)')\nplt.yscale('log')\nplt.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Many very low values, let's zoom in on the lower values"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nplt.subplot(211)\nplt.hist(train_df[train_df['meter_reading']<=268].meter_reading.dropna(),30)\nplt.title('distribution of meter_reading in train data')\nplt.xlabel('value of meter reading')\nplt.ylabel('count')\nplt.subplot(212)\nplt.hist(train_df[train_df['meter_reading']<=268].meter_reading.dropna(),30)\nplt.title('logscale distribution of meter_reading in train data, focus on quartiles 1 to 3')\nplt.xlabel('value of meter reading')\nplt.ylabel('log(count)')\nplt.yscale('log')\nplt.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Still a high bar near zero, let's count actual zero reading in the training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"number_of_zero_readings = train_df[train_df['meter_reading']==0].meter_reading.count()\nnumber_of_readings = train_df['meter_reading'].count()\nprint('there are ',100*number_of_zero_readings/number_of_readings,'% zero readings in train_data')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Comparing meter types"},{"metadata":{},"cell_type":"markdown","source":"> There are four different meter types and are displayed below:\n>     0: electricity\n>     1: chilledwater\n>     2: steam\n>     3: hotwater"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(y='meter',data = train_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(20, 5))\n\nlow_reading = train_df['meter_reading']<=30\nnonnul_reading = train_df['meter_reading']!=0\n\nelectricity_reading = train_df['meter']==0\nsns.distplot(train_df[low_reading & electricity_reading & nonnul_reading].meter_reading, color = 'blue',ax =ax)\n\nchillwater_reading = train_df['meter']==1\nsns.distplot(train_df[low_reading & chillwater_reading & nonnul_reading].meter_reading, color = 'green',ax =ax)\n\nsteam_reading = train_df['meter']==2\nsns.distplot(train_df[low_reading & steam_reading & nonnul_reading].meter_reading, color = 'red',ax =ax)\n\nhotwater_reading = train_df['meter']==3\nsns.distplot(train_df[low_reading & hotwater_reading & nonnul_reading].meter_reading, color = 'orange',ax =ax)\n\nplt.xlabel('meter_reading', fontsize=16)\nplt.legend(['elec', 'chill water','steam','hotwater'])\nplt.show()    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Notice that while electricity anbd chillwater are continuous, **steam and hotwater measurements are by step of approx 5**"},{"metadata":{},"cell_type":"markdown","source":"## Weather data"},{"metadata":{},"cell_type":"markdown","source":"### looking at weatehr data correlation"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train = reduce_mem_usage(pd.read_csv('../input/ashrae-energy-prediction/weather_train.csv'))\nplt.figure(figsize = (6, 6))\nsns.heatmap(weather_train.groupby('site_id').count(), annot = True)\n# No precipitations data for site 1, 5 and 12\n# No sea-level pressure for site 5\ndel weather_train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Distribution + comparing test & train data"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_distribution(train_df,test_df,column):\n    fig, ax = plt.subplots(figsize=(20, 5))\n    sns.distplot(train_df[column].dropna(), color = 'red',ax =ax).set_title(column,fontsize = 10)\n    sns.distplot(test_df[column].dropna(), color = 'blue',ax =ax).set_title(column,fontsize = 10)\n    plt.xlabel(column, fontsize=16)\n    plt.legend(['training data', 'test data'])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_distribution(train_df,test_df,'cloud_coverage')\nplot_distribution(train_df,test_df,'air_temperature')\nplot_distribution(train_df,test_df,'dew_temperature')\nplot_distribution(train_df,test_df,'precip_depth_1_hr')\nplot_distribution(train_df,test_df,'sea_level_pressure')\nplot_distribution(train_df,test_df,'wind_direction')\nplot_distribution(train_df,test_df,'wind_speed')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distributions are similiar between test and learn data, except for precipitation -> handle with care"},{"metadata":{},"cell_type":"markdown","source":"### Is it all in the northern hemisphere ?\nOtherwise, might need to move data by 6 months for given site_id's"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_df = train_df[['timestamp','building_id','site_id','air_temperature']].copy()\ntemp_df['date'] = temp_df.timestamp.dt.date\ntemp_df = temp_df.groupby(['site_id','date'],as_index=False).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,ax=plt.subplots(4,4,figsize=(25,30))\nplt.subplots_adjust(hspace=0.5)\nfor k in train_df.site_id.unique():\n    temp_df[temp_df['site_id']==k].plot(x='date', y='air_temperature', ax=ax[(k//4),(k%4)], title=f'site {k}', figure=fig)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can observe a few things :\n* all peak in the middle of the year -> northen hemisphere\n* there are some striking similarities between a few sites : 14,15,11,7,6 ahve the same freezing february day\nToDo : cluster the sites based on the distance between their air temperature to add a feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"del temp_df\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Looking for seasonality\nidea : try to cluster buildings into climate types could improve model"},{"metadata":{"trusted":true},"cell_type":"code","source":"#group per day of the week\nweekday_group_df = train_df[['timestamp','meter_reading','primary_use','meter']].copy()\nweekday_group_df['weekday'] = weekday_group_df['timestamp'].dt.weekday\nweekday_group_df = weekday_group_df.groupby(['primary_use','meter','weekday'],group_keys=False,as_index=False).sum()\nweekday_group_df = reduce_mem_usage(weekday_group_df)\nweekday_group_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weekday_group_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"primary_use_list = train_df['primary_use'].unique()\nlen(primary_use_list)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(16, 4,figsize=(25,50))\nplt.subplots_adjust(hspace=0.5)\ni=0\nfor use in primary_use_list:\n    use_filter = weekday_group_df['primary_use']==use\n    for k in range(0,4):\n        meter_filter = weekday_group_df['meter']==k\n        if weekday_group_df[use_filter & meter_filter].shape[0] != 0:\n            weekday_group_df[use_filter & meter_filter].plot(x='weekday' , y='meter_reading' , ax=ax[i,k] , figure=fig , title=f'{use} - meter type {k}')\n    i=i+1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can observe :\n* Strong weekly trends\n* Different for each primary use and meter type\n-> primary_use is a strongly descriptive variable"},{"metadata":{"trusted":true},"cell_type":"code","source":"del weekday_group_df,i,fig,ax\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Per primary use"},{"metadata":{"trusted":true},"cell_type":"code","source":"#group per day of the year\ndaily_group_df = train_df[['timestamp','meter_reading','primary_use','meter','site_id']].copy()\ndaily_group_df['date'] = daily_group_df['timestamp'].dt.date\ndaily_group_df = daily_group_df.groupby(['site_id','primary_use','meter','date'],group_keys=False,as_index=False).sum()\ndaily_group_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#let'ts look at yearly trends per meter type\nplt.figure(figsize=(20,8))\nax = sns.lineplot(x='date', y=\"meter_reading\", data=daily_group_df,hue = 'meter').set_title('Energy consumption per meter type')\nplt.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Electricity is the main source of energy, followed by steam. they are strongly seasonal"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_per_primary_use(df,primary_use):\n    plt.figure(figsize=(20,8))\n    ax = sns.lineplot(x='date', y=\"meter_reading\", data=df[df['primary_use']==primary_use],hue = 'meter').set_title('Energy consumption per type for '+primary_use+' buildings')\n    plt.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"primary_use_list = train_df['primary_use'].unique()\nprimary_use_list","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for primary_use in primary_use_list:\n    plot_per_primary_use(daily_group_df,primary_use)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Per site ID (geography)"},{"metadata":{"trusted":true},"cell_type":"code","source":"colour_palette_16=sns.color_palette(n_colors=16)\ncolour_palette_4=sns.color_palette(n_colors=4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,8))\nax = sns.lineplot(x='date', y=\"meter_reading\", data=daily_group_df,hue = 'site_id',legend= 'full',palette = colour_palette_16).set_title('Energy consumption per site id')\nplt.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(4, 4,figsize=(25,30))\nplt.subplots_adjust(hspace=0.5)\nfor i in range(0,16):\n    daily_group_df[(daily_group_df.site_id==i) & daily_group_df.meter==0].plot(y='meter_reading',x='date',title=f'site:{i}',figure=fig,ax=ax[(i//4),(i%4)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_per_site_id(df,site_id):\n    plt.figure(figsize=(20,5))\n    ax = sns.lineplot(x='date', y=\"meter_reading\", data=df[df['site_id']==site_id],hue = 'meter').set_title('Energy consumption per type for '+str(site_id)+' site id')\n    plt.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for n in range(0,15):\n    plot_per_site_id(daily_group_df,site_id=n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del daily_train_df\ngc()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# NA values processing"},{"metadata":{},"cell_type":"markdown","source":"#### Where do we have missing values ?"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#% of missing values\n100*train_df.isna().sum()/train_df.shape[0]"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#let's check test data\n100*test_df.isna().sum()/test_df.shape[0]"},{"metadata":{},"cell_type":"markdown","source":"Similar missing values for train and test"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#for now let's drop them\nNA_cols = ['year_built','floor_count','cloud_coverage','precip_depth_1_hr']\ntest_df.drop(NA_cols,axis=1,inplace=True)\ntrain_df.drop(NA_cols,axis=1,inplace=True)"},{"metadata":{},"cell_type":"markdown","source":"#### for the other ones, let's fill with mean value"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"few_NA_cols = ['air_temperature','dew_temperature','sea_level_pressure','wind_direction','wind_speed']"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#fill brutally with column mean, need to find better segmentation\nfor col in few_NA_cols:\n    train_df[col].fillna(train_df[col].mean,inplace=True)\n    test_df[col].fillna(test_df[col].mean,inplace=True)"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#% of missing values in train\n100*train_df.isna().sum()/train_df.shape[0]"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#let's check test data\n100*test_df.isna().sum()/test_df.shape[0]"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#I need to group the data by building for interpolation to make sense\n#ToDo : make dict with all NA n% over 40%\n#maxed_bluidings ={}\nfor col in few_NA_cols:\n    max_NA = 0\n    for building_k in train_df['building_id'].unique():\n        only_building_k = train_df['building_id']==building_k\n        temp = max(100*train_df[only_building_k].isna().sum()/train_df[only_building_k].shape[0])\n        if temp > 40:\n            max_NA=temp\n            maxed_building = building_k\n            maxed_buildings[building_k]=max_NA\n    print('in column ',col)\n    print('Maximum % of missing values for a same building is : ',max_NA)\n    print('Achieved by :',maxed_building)"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"def interpolate_NA_values(df,column):\n    temp_df = df.groupby('building_id')[column].interpolate()\n    #for building_k in range(0,df['building_id'].unique):\n    #    only_building_k = df['building_id']==building_k\n    #    df[only_building_k].column.interpolate()\n    return temp_df"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"train_df = interpolate_NA_values(train_df,'air_temperature')"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"100*train_df.isna().sum()/train_df.shape[0]"},{"metadata":{},"cell_type":"markdown","source":"#### Weather - done above"},{"metadata":{"trusted":true},"cell_type":"code","source":"# For the following fillings, I was thinking to add another separation into Morning / Afternoon / Evening / Night, in case there could be cycles of the data\n#weather_train['sea_level_pressure'] = weather_train.groupby('site_id').transform(lambda x: x.fillna(x.mean()))['sea_level_pressure']\n#weather_train['air_temperature'] = weather_train.groupby('site_id').transform(lambda x: x.fillna(x.mean()))['air_temperature']\n#weather_train['dew_temperature'] = weather_train.groupby('site_id').transform(lambda x: x.fillna(x.mean()))['dew_temperature']\n#weather_train['wind_direction'] = weather_train.groupby('site_id').transform(lambda x: x.fillna(x.mean()))['wind_direction']\n#weather_train['wind_speed'] = weather_train.groupby('site_id').transform(lambda x: x.fillna(x.mean()))['wind_speed']\n\n# MUST BE LAST ACTION \n# Sea-level pressure : WW mean = 1013.25 (https://en.wikipedia.org/wiki/Atmospheric_pressure)\n#weather_train.fillna(1013.25, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Features generation"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Year were the building was built. Arbitrary thresholds for now.\n#train_df['is_old'] = train_df['year_built'] <= 1990\n#buildings['is_very_old'] = buildings['year_built'] <= 1950\n#buildings.drop('year_built', axis = 1, inplace = True) #@greg why drop it ?","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Number of floors. I've tried to distinguish buildings with relatively big and small floors.\n#buildings_floor = buildings.dropna()[['square_feet', 'floor_count']]\n#buildings_floor['avg_area'] = buildings_floor['square_feet'] / buildings_floor['floor_count']\n#avg_area = buildings_floor['avg_area'].mean()\n#avg_area","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#buildings['big_floors'] = (buildings['square_feet'] / buildings['floor_count']) > 3 * avg_area\n#buildings['small_floors'] = (buildings['square_feet'] / buildings['floor_count']) < avg_area / 3\n#buildings['few_floors'] = buildings['floor_count'] == 1\n#buildings['many_floors'] = buildings['floor_count'] >= 8\n#buildings.drop('floor_count', axis = 1, inplace = True)","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}