{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Introduction\nThis kernel is created with respect to the Kaggle competition [ASHRAE - Great Energy Predictor III](https://www.kaggle.com/c/ashrae-energy-prediction). In this competition we need to build a models to predict metered building energy usage in the following areas: chilled water, electric, hot water, and steam meters.\n
"},{"metadata":{},"cell_type":"markdown","source":"In this kernel, we have done an **Exploratory Data Analysis(EDA)** of the entire dataset and have also created a model with the help of **LightGBM Regressor**. We have also performed data-frame memory usage reduction, data pre-processing(Dealing with the NaN values) and feature Engineering before training the model."},{"metadata":{},"cell_type":"markdown","source":"## Load Packages and Data"},{"metadata":{"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)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport gc\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import explained_variance_score,mean_absolute_error,mean_squared_error,r2_score\nfrom sklearn.ensemble import RandomForestRegressor\nfrom catboost import CatBoostRegressor\nfrom sklearn.model_selection import KFold, StratifiedKFold\nimport lightgbm as lgb\nfrom tqdm import tqdm_notebook as tqdm\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":{"trusted":true},"cell_type":"code","source":"BASE_PATH = '/kaggle/input/ashrae-energy-prediction/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def Pre_process_data(df,col):\n '''\n Input: Data-frame and Column name.\n Operation: Fills the nan values with the minimum value in their respective column.\n and if the column is \"dew_temperature\" or \"air_temperature\" it uses forward \n fill and backward fill.\n Output: Returns the pre-processed data-frame.\n '''\n #df['primary_use'] = df['primary_use'].astype(\"category\").cat.codes\n print(\"Name of column with NaN: \"+str(col))\n print(df[col].value_counts(dropna=False, normalize=True).head())\n if col=='dew_temperature' or col=='air_temperature':\n df[col] = df[col].ffill(axis = 0) \n df[col] = df[col].bfill(axis = 0)\n else:\n df[col] = df[col].fillna(df[col].min())\n #df.drop(['building_id'], axis=1, inplace=True)\n return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Based on this great kernel https://www.kaggle.com/arjanso/reducing-dataframe-memory-size-by-65\ndef reduce_mem_usage(df):\n '''\n Input - data-frame.\n Operation - Reduce memory usage of the data-frame.\n '''\n start_mem_usg = df.memory_usage().sum() / 1024**2 \n print(\"Memory usage of properties dataframe is :\",start_mem_usg,\" MB\")\n #NAlist = [] # Keeps track of columns that have missing values filled in. \n for col in df.columns:\n if df[col].dtype != object: # Exclude strings \n # Print current column type\n print(\"******************************\")\n print(\"Column: \",col)\n print(\"dtype before: \",df[col].dtype) \n # make variables for Int, max and min\n IsInt = False\n mx = df[col].max()\n mn = df[col].min()\n print(\"min for this col: \",mn)\n print(\"max for this col: \",mx)\n # Integer does not support NA, therefore, NA needs to be filled\n if not np.isfinite(df[col]).all(): \n #NAlist.append(col)\n df = Pre_process_data(df,col)\n \n # test if column can be converted to an integer\n asint = df[col].fillna(0).astype(np.int64)\n result = (df[col] - asint)\n result = result.sum()\n if result > -0.01 and result < 0.01:\n IsInt = True \n # Make Integer/unsigned Integer datatypes\n if IsInt:\n if mn >= 0:\n if mx < 255:\n df[col] = df[col].astype(np.uint8)\n elif mx < 65535:\n df[col] = df[col].astype(np.uint16)\n elif mx < 4294967295:\n df[col] = df[col].astype(np.uint32)\n else:\n df[col] = df[col].astype(np.uint64)\n else:\n if mn > np.iinfo(np.int8).min and mx < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif mn > np.iinfo(np.int16).min and mx < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif mn > np.iinfo(np.int32).min and mx < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif mn > np.iinfo(np.int64).min and mx < np.iinfo(np.int64).max:\n df[col] = df[col].astype(np.int64) \n # Make float datatypes 32 bit\n else:\n df[col] = df[col].astype(np.float32)\n \n # Print new column type\n print(\"dtype after: \",df[col].dtype)\n print(\"******************************\")\n # Print final result\n print(\"___MEMORY USAGE AFTER COMPLETION:___\")\n mem_usg = df.memory_usage().sum() / 1024**2 \n print(\"Memory usage is: \",mem_usg,\" MB\")\n print(\"This is \",100*mem_usg/start_mem_usg,\"% of the initial size\")\n return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_train = pd.read_csv(BASE_PATH+'weather_train.csv')\nweather_train = weather_train.set_index('timestamp')\nweather_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_test = pd.read_csv(BASE_PATH+'weather_test.csv')\nweather_test = weather_test.set_index('timestamp')\nweather_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## EDA for weather data-set"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(17,15)})\nfig, ax = plt.subplots(2, 1)\nweather_train['air_temperature'].plot(marker='.',color='r', alpha=0.1,ax=ax[0])\nweather_test['air_temperature'].plot(marker='.', alpha=0.1,ax=ax[1])\nax[0].set_ylabel('Air Temperature')\nax[1].set_ylabel('Air Temperature')\n#weather_test['air_temperature'].plot(marker='.', alpha=0.1);\n#weather_train['air_temperature'].plot(marker='.', alpha=0.3);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(17,15)})\nfig, ax = plt.subplots(2, 1)\nweather_train['cloud_coverage'].plot(marker='o',color='r',alpha=0.3,linestyle='-',ax=ax[0])\nweather_test['cloud_coverage'].plot(marker='o',linestyle='-', alpha=0.3,ax=ax[1])\nax[0].set_ylabel('Cloud Coverage')\nax[1].set_ylabel('Cloud Coverage')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(17,15)})\nfig, ax = plt.subplots(2, 1)\nweather_train['dew_temperature'].plot(marker='.',color='r',alpha=0.3,ax=ax[0])\nweather_test['dew_temperature'].plot(marker='.', alpha=0.3,ax=ax[1])\nax[0].set_ylabel('Dew Temperature')\nax[1].set_ylabel('Dew Temperature')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(16,9)})\nweather_train['air_temperature'].plot(marker='o',color='r',alpha=0.3,label='air_temperature')\nweather_train['dew_temperature'].plot(marker='.',color='y', alpha=0.1,label='dew_temperature')\nplt.legend()\n#ax[0].set_ylabel('Dew Temperature')\n#ax[1].set_ylabel('Dew Temperature')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(16,9)})\nweather_test['air_temperature'].plot(marker='o',color='r',alpha=0.3,label='air_temperature')\nweather_test['dew_temperature'].plot(marker='.',color='y', alpha=0.1,label='dew_temperature')\nplt.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.set(rc={'figure.figsize':(17,15)})\nk1 = sns.lineplot(weather_train['wind_direction'],weather_train['wind_speed'],label='2016-2017')\nk2 = sns.lineplot(weather_test['wind_direction'],weather_test['wind_speed'],label='2017-2018')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(BASE_PATH+'train.csv')\n#train = train.set_index('timestamp')\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building = pd.read_csv(BASE_PATH+'building_metadata.csv')\nbuilding.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## EDA for our Training Data-frame"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"train_final = train.merge( building, left_on = \"building_id\",right_on = \"building_id\", how = \"left\")\ntrain_final = reduce_mem_usage(train_final)\n#train_final.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"train_final = train_final.merge( weather_train, on=['site_id', 'timestamp'], how='left')\ntrain_final = reduce_mem_usage(train_final)\ntrain_final.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,20))\nsns.countplot(y='primary_use', data=train_final)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plt.figure(figsize=(20,20))\n#N = 10000\n#k=sns.barplot(train_final['meter_reading'][0:N], train_final['primary_use'][0:N])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_final = train_final.set_index('timestamp')\nsns.set(rc={'figure.figsize':(20,15)})\ntrain_final['meter_reading'].plot(marker='o',c='r',alpha=0.3,label='2016-2017')\nplt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> **I would come up with more EDA representing the training data frame. Plotting some of the graphs is taking a lot of time sometimes even one hour and more. I would come up with a solution to plot those graphs and add it to the kernel**"},{"metadata":{"trusted":true},"cell_type":"code","source":"del train,\ndel k1\ndel k2\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_final = train_final.reset_index()\n'''train_final.drop(['timestamp'], \n axis=1, inplace=True)'''\ntrain_final.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del building, weather_train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering\nWe have added more features with respect to features that are already present in the data frame and turned categorical values into integer values.\nFeature Engineering done:\n* Transforming the \"Primary_use\" column into an unique integer column.\n* Adding only the month from the \"timestamp\" column.
\n Logic: Weather change takes place between an interval of months, therefore energy consumption also differs month-wise.\n* wind_load - The generic formula for wind load is F = A x P x Cd where F is the force or wind load, A is the projected area of the object, P is the wind pressure, and Cd is the drag coefficient.\n* temp_diff - Difference between dew temperature and air_temperature.\n* Floor_Area - Square area of the building multiplied with number of floors of the building."},{"metadata":{"trusted":true},"cell_type":"code","source":"def Feature_Engineering(df):\n '''\n Input: Data-frame\n Operation: Adding features to the data-frame\n Output: Inproved Data frame.\n '''\n df['primary_use'] = df['primary_use'].astype(\"category\").cat.codes\n df[\"timestamp\"] = pd.to_datetime(df[\"timestamp\"])\n df[\"month\"] = df[\"timestamp\"].dt.month.astype(np.uint8)\n df.drop(['building_id','timestamp'], axis=1, inplace=True)\n df['wind_load'] = df['sea_level_pressure']*(df['wind_speed']**2)\n df['windspeed_mean'] = df['wind_speed'] / df.groupby(['wind_direction'])['wind_speed'].transform('mean')\n df['windspeed_diff_std'] = df['wind_speed'] - df.groupby(['wind_direction'])['wind_speed'].transform('std')\n df['temp_diff'] = df['dew_temperature'] - df['air_temperature']\n df['Floor_Area'] = df['floor_count']*df['square_feet']\n \n return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_final = Feature_Engineering(train_final)\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_final.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Correlation matrix(Heat map)"},{"metadata":{"trusted":true},"cell_type":"code","source":"f,ax = plt.subplots(figsize=(18,18))\nsns.heatmap(train_final.corr(), annot=True, linewidths=.5, fmt= '.2f',ax=ax,cmap=\"YlGnBu\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lesson Learned\nAccording to my knowledge, all types of data pre-processing, and data engineering are always done with the input labels. The target label is kept as it is.
In this case also I assumed the same thing, but no matter how many times I fit the model and validated its performance there was a huge RMSE(Root Mean Squared Error) and MAE(Mean Absolute Error). I tried to fit the model with different regressor algorithms but the same thing was repeated again and again.
Then I went through some of the notebooks in the Notebooks panel of the competition and found that the users are transforming the target variable in this case to obtain better accuracy. It was a completely new thing for me to learn that the target variable could also be processed according to our requirement to get better results.
According to my opinion the main reason for this problem to occur was because the ML model was not able to find proper pattern in the data to generate the target data, but once you transform the target data the model is able to find a pattern in the input data corresponding to the target data and hence the errors are reduced while predicting the target variable.\n
In this case we have transformed the target value corresponding to the equation **log(1 + x)**."},{"metadata":{"trusted":true},"cell_type":"code","source":"y = np.log1p(train_final[\"meter_reading\"])\ntrain_final.drop(\"meter_reading\", axis=1, inplace=True)\nX = train_final","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#x_train, x_val, y_train, y_val = train_test_split(X, y,test_size=0.3,random_state=42)\ndel train_final\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model Fit\nI have tried to fit and validate the data with many regressor algorithms, but LightGBM and CatBoost regressor turned out to be the most perfect to solve this problem. With CatBoost regressor we achieved a **Variance score** of **0.711**, **MAE** of **0.734**, **MSE** of **1.33**, **R-Score** of **0.711** over 2000 iteration. I have commented out the portion where the CatBoost regressor is used for model fit and evaluation.
I am thinking of creating another kernel where I would evaluate and show the performance of Catboost with a higher number of iterations.
With LightGBM regressor we achieved a **Variance score** of **0.87**, **MAE** of **0.448**, **MSE** of **0.594**, **R-Score** of **0.87** over 2000 iteration.
\nThe performance of LightGBM could be seen in the output of the cell down below."},{"metadata":{"trusted":true},"cell_type":"code","source":"folds = 4\nseed = 42\nnum_epochs =300\nnum_batch_size = 1080\nkfold = KFold(n_splits = folds, shuffle = True, random_state = seed)\nmodels = []\nfor train_index, val_index in kfold.split(X,y):\n X_train, X_val = X.loc[train_index,:], X.loc[val_index,:] \n y_train, y_val = y[train_index], y[val_index]\n print({'train size':len(X_train), 'eval size':len(X_val)})\n \n lgb_train = lgb.Dataset(X_train, y_train)\n lgb_eval = lgb.Dataset(X_val, y_val)\n \n clf = lgb.LGBMRegressor(metric='rmse',\n learning_rate=0.4,\n #feature_fraction= 0.9,\n n_estimators= 600,\n subsample=0.3, # batches of 25% of the data\n subsample_freq=1,\n #lgb_train,\n #num_boost_round=2000,\n #valid_sets=(lgb_train, lgb_eval),\n #early_stopping_rounds=20,\n verbose_eval = 500\n )\n clf.fit(X_train, y_train,\n eval_set=[(X_val, y_val)],\n #early_stopping_rounds=50,\n verbose=200)\n #clf = CatBoostRegressor(iterations=2000,depth= 9,random_seed = 23,\n # task_type = \"GPU\")\n #clf.fit(X_train,y_train)\n y_pred = clf.predict(X_val)\n print(\"\\nVariance_Score\\t:\"+str(explained_variance_score(y_val,y_pred)))\n print(\"Mean_Absolute_Error\\t:\"+str(mean_absolute_error(y_val,y_pred)))\n print(\"Mean_Squared_Error\\t:\"+str(mean_squared_error(y_val,y_pred)))\n print(\"R2-Score\\t:\"+str(r2_score(y_val,y_pred)))\n models.append(clf)\n del X_train, X_val, y_train, y_val, clf, lgb_train, lgb_eval\n gc.collect()\n\n print (20*'---')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15,15))\nplt.bar(range(len(models[0].feature_importances_)), models[0].feature_importances_)\nplt.title(\"Feature Importance\")\nplt.xticks(np.arange(len(X.columns)),X.columns, rotation=90)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"building = pd.read_csv(BASE_PATH+'building_metadata.csv')\ntest = pd.read_csv(BASE_PATH+'test.csv')\n#test = test.set_index('timestamp')\n#test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"test_final = test.merge( building, left_on = \"building_id\",right_on = \"building_id\", how = \"left\")\ntest_final = reduce_mem_usage(test_final)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"test_final = test_final.merge( weather_test, on=['site_id', 'timestamp'], how='left')\ntest_final = reduce_mem_usage(test_final)\ntest_final.drop(['row_id'], axis=1, inplace=True)\n#test_final.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del building,test,X,y,weather_test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_final = Feature_Engineering(test_final)\ntest_final.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission\nPreparing for submission. The submission code is inpired from the kernel [Starter EDA and Feature selection ASHRAE3](https://www.kaggle.com/hmendonca/starter-eda-and-feature-selection-ashrae3)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# split test data into batches\nset_size = len(test_final)\niterations = 100\nbatch_size = set_size // iterations\n\nprint(set_size, iterations, batch_size)\nassert set_size == iterations * batch_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meter_reading = []\nfor i in tqdm(range(iterations)):\n pos = i*batch_size\n fold_preds = [np.expm1(model.predict(test_final[test_final.columns].iloc[pos : pos+batch_size])) for model in models]\n meter_reading.extend(np.mean(fold_preds, axis=0))\n\n#print(len(meter_reading))\nassert len(meter_reading) == set_size","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**numpy.clip()** function is used to Clip (limit) the values in an array.\n\nGiven an interval, values outside the interval are clipped to the interval edges. For example, if an interval of [0, 1] is specified, values smaller than 0 become 0, and values larger than 1 become 1."},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(BASE_PATH+'sample_submission.csv')\nsubmission['meter_reading'] = np.clip(meter_reading, a_min=0, a_max=None) # clip min at zero\ndel meter_reading, test_final\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission2.csv', index=False)\nsubmission.head(9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## References\nThese are the following kernels which helped me a lot while understanding the problem and also I think you must explore them to know different types of EDA that could be done with the data and different approaches to solving the prediction problem.\n* [Starter EDA and Feature selection ASHRAE3](https://www.kaggle.com/hmendonca/starter-eda-and-feature-selection-ashrae3)\n* [Simple LightGBM LB 1.24](https://www.kaggle.com/isaienkov/simple-lightgbm-lb-1-24)\n* [Starter ⚡ Great Energy Predictor ](https://www.kaggle.com/jesucristo/starter-great-energy-predictor)\n* [ASHRAE - Energy prediction](https://www.kaggle.com/allunia/ashrae-energy-prediction)"},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"## Thank you\n### Thank you for staying with me throughout the kernel. Please **up-vote** if you think this kernel was informative or if you like it.\n### Also comment down below about how much did you like this kernel, or what all improvements that you think could be added to this kernel."}],"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}