Large Scale Energy Anomaly Detection
\n\nIn this Notebook, I had implemented the anomaly detection model to predict whether the energy usage in a building is anomalous or not. This will help to save a lot of energy. In this notebook I had used the train_features.csv data to build a model. Look at my this notebook in which I had implement a model by just using the train.csv data.
\n\n\n
TABLE OF CONTENTS
\n \n* [1. Importing Libraries](#1)\n \n* [2. Exploratory Data Analysis](#2)\n \n* [3. Highly Correlated Features](#3)\n \n* [4. Model Development](#4) \n \n* [5. The End](#5) \n\n\n#IMPORTING LIBRARIES
","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sb\nfrom tqdm import tqdm, trange\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder, StandardScaler\nfrom sklearn import metrics\nfrom xgboost import XGBClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom imblearn.over_sampling import RandomOverSampler\n\nimport lightgbm as lgb\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.graph_objs as go\nimport plotly\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\nimport cufflinks as cf\ncf.set_config_file(offline=True)\n# Input data files are available i\n\nimport datetime\n\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom keras import layers\n\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:03.975127Z","iopub.execute_input":"2022-06-29T02:13:03.975522Z","iopub.status.idle":"2022-06-29T02:13:11.31023Z","shell.execute_reply.started":"2022-06-29T02:13:03.97542Z","shell.execute_reply":"2022-06-29T02:13:11.309234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n#Exploratory Data Analysis
\nWe have been provided with a whole lot of features. As per the Author of the competition data and cleaned and one can focused on the model building but till we are not sure with the features provided and what they represent how can we go forward.
","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('../input/energy-anomaly-detection/train_features.csv')\ntrain.shape","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:11.312168Z","iopub.execute_input":"2022-06-29T02:13:11.312546Z","iopub.status.idle":"2022-06-29T02:13:29.69757Z","shell.execute_reply.started":"2022-06-29T02:13:11.312509Z","shell.execute_reply":"2022-06-29T02:13:29.696989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:29.698529Z","iopub.execute_input":"2022-06-29T02:13:29.699261Z","iopub.status.idle":"2022-06-29T02:13:29.729328Z","shell.execute_reply.started":"2022-06-29T02:13:29.699229Z","shell.execute_reply":"2022-06-29T02:13:29.728436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for col in train.columns:\n k = train[col].isnull().sum()\n if k > 0:\n print(f'{col} : {k}')","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:29.731477Z","iopub.execute_input":"2022-06-29T02:13:29.731738Z","iopub.status.idle":"2022-06-29T02:13:31.378151Z","shell.execute_reply.started":"2022-06-29T02:13:29.731705Z","shell.execute_reply":"2022-06-29T02:13:31.377066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this data as well only meter_reading has null values.
There are no such entries where it is an anomalous case and meter_reading is null.
Very nice observation indeed has been shown in this discussion.
","metadata":{}},{"cell_type":"code","source":"for idx in tqdm(train['building_id'].unique()):\n if train[train['building_id']==idx]['primary_use'].nunique()>1:\n print(idx)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:31.499265Z","iopub.execute_input":"2022-06-29T02:13:31.499771Z","iopub.status.idle":"2022-06-29T02:13:35.865513Z","shell.execute_reply.started":"2022-06-29T02:13:31.499722Z","shell.execute_reply":"2022-06-29T02:13:35.864519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_by_use = train.groupby('primary_use').mean()['meter_reading']\nmean_by_use.plot.bar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:35.867062Z","iopub.execute_input":"2022-06-29T02:13:35.867701Z","iopub.status.idle":"2022-06-29T02:13:37.161797Z","shell.execute_reply.started":"2022-06-29T02:13:35.86765Z","shell.execute_reply":"2022-06-29T02:13:37.161123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each building provides only one type of facility/usecase and above chart shows the mean usage by use case throughout the year with minimum at the religious places and highest at the Education institutions.
\nNow let's explore what are the different features has been provided in this data as compare to the train.csv.
","metadata":{}},{"cell_type":"code","source":"ints = []\nobjects = []\nfloats = []\n\nfor col in train.columns:\n if train[col].dtype == int:\n ints.append(col)\n elif train[col].dtype == object:\n objects.append(col)\n else:\n floats.append(col)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:37.162862Z","iopub.execute_input":"2022-06-29T02:13:37.163651Z","iopub.status.idle":"2022-06-29T02:13:37.169302Z","shell.execute_reply.started":"2022-06-29T02:13:37.163612Z","shell.execute_reply":"2022-06-29T02:13:37.168683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Analysis of the columns with int datatype.
Columns from hour to is_holiday is nothing but the expanded form of the timestamp so, we can remove that feature.
year_built - To get the exact year_built add 1900 to the same cloud_coverage - It has all values from 0 - 9. But the value 255 in both year_built and cloud_coverage are replacement for the null values.No use of columns in the modelling with only one values for all the entries.
","metadata":{}},{"cell_type":"code","source":"train['date'] = pd.to_datetime(train['timestamp']).dt.date\ntrain['meterReadings_daily_std'] = train.groupby(['building_id','date'])['meter_reading'].transform('std')","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:42.94734Z","iopub.execute_input":"2022-06-29T02:13:42.94795Z","iopub.status.idle":"2022-06-29T02:13:44.598751Z","shell.execute_reply.started":"2022-06-29T02:13:42.947899Z","shell.execute_reply":"2022-06-29T02:13:44.597927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.drop(['timestamp', 'year', 'gte_meter'], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:13:44.600364Z","iopub.execute_input":"2022-06-29T02:13:44.600817Z","iopub.status.idle":"2022-06-29T02:13:45.50247Z","shell.execute_reply.started":"2022-06-29T02:13:44.600782Z","shell.execute_reply":"2022-06-29T02:13:45.50144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Analysis of the columns with object datatype.
It seems like these features are nothing but the concatenation of the building_id with the month, hour or other time realted features. In my opinion these are not going to help us as none of the values will overlap with the test data as the building id's are totally different for the test data set. So, the distribution of the data will vary.
Analysis of the columns with float datatype.
These features from it's name seems like they are derived by feature engineering from the time and the building_id except for some columns let's keep them as of now. It has been given in the data section that these features has been used by the winners so, may be these features proved to be helpful.
Highly Correlated Features
\nBy the name of the columns it is clear that they have been derived by combining two or more succint features. In such cases there is possibility of having high correlation between the developed feature and the original feature. And as much as I know highly correlated features doesn't help model.
","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(15,15))\nsb.heatmap(train.corr() > 0.95, annot=True, cbar=False)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:15:41.976054Z","iopub.execute_input":"2022-06-29T02:15:41.976414Z","iopub.status.idle":"2022-06-29T02:16:01.125911Z","shell.execute_reply.started":"2022-06-29T02:15:41.976363Z","shell.execute_reply":"2022-06-29T02:16:01.124948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.columns","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:16:01.127453Z","iopub.execute_input":"2022-06-29T02:16:01.128401Z","iopub.status.idle":"2022-06-29T02:16:01.1353Z","shell.execute_reply.started":"2022-06-29T02:16:01.128355Z","shell.execute_reply":"2022-06-29T02:16:01.134421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hcorr = [#'building_id',\n #'site_id',\n #'year_built',\n #'square_feet',\n 'primary_use',\n #'floor_count',\n #'precip_depth_1_hr', 'sea_level_pressure',\n #'wind_direction', 'wind_speed',\n 'gte_meter_hour',\n 'gte_meter_weekday',\n 'gte_meter_month',\n 'gte_meter_building_id',\n 'gte_meter_primary_use',\n 'gte_meter_site_id',\n 'gte_meter_building_id_hour',\n 'gte_meter_building_id_weekday',\n 'gte_meter_building_id_month',\n 'air_temperature_mean_lag73',\n 'air_temperature_max_lag73',\n 'air_temperature_min_lag73',\n 'air_temperature_mean_lag7',\n 'air_temperature_max_lag7',\n 'air_temperature_min_lag7']\n\ntrain = train.drop(hcorr, axis=1)\ntrain.shape","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:16:01.136717Z","iopub.execute_input":"2022-06-29T02:16:01.136973Z","iopub.status.idle":"2022-06-29T02:16:01.318338Z","shell.execute_reply.started":"2022-06-29T02:16:01.136942Z","shell.execute_reply":"2022-06-29T02:16:01.317235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15,15))\nsb.heatmap(train.corr() > 0.75, annot=True, cbar=False)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:16:01.319915Z","iopub.execute_input":"2022-06-29T02:16:01.320267Z","iopub.status.idle":"2022-06-29T02:16:10.824832Z","shell.execute_reply.started":"2022-06-29T02:16:01.320218Z","shell.execute_reply":"2022-06-29T02:16:10.823856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, there are no highly correlated features and we are good to go and build our anomaly detection model.
","metadata":{}},{"cell_type":"markdown","source":"\n#Model Development
\nTraining an XGBClassifier on millions of rows take hell lot of time that is why I have implemented a neural network model so, that the model can go through all the examples in feasible time. By using the RandomOverSampler I had balanced the negative and the positive exmaples so, that model can perform well for the negative as well as positive examples.
","metadata":{}},{"cell_type":"code","source":"test = pd.read_csv('../input/energy-anomaly-detection/test_features.csv', index_col=0)\ntest = impute_nulls(test)\ntest['cloud_coverage'] = test['cloud_coverage'].replace({255:10})\n\ntest['date'] = pd.to_datetime(test['timestamp']).dt.date\ntest['meterReadings_daily_std'] = test.groupby(['building_id','date'])['meter_reading'].transform('std')\n\ntest = test.drop(['timestamp', 'year', 'gte_meter'] + objects[2:] + hcorr, axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:16:10.829061Z","iopub.execute_input":"2022-06-29T02:16:10.829434Z","iopub.status.idle":"2022-06-29T02:19:22.275999Z","shell.execute_reply.started":"2022-06-29T02:16:10.8294Z","shell.execute_reply":"2022-06-29T02:19:22.274867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train['meter_reading'] = train['meter_reading']/train.groupby('building_id')['meter_reading'].transform('std')\n#test['meter_reading'] = test['meter_reading']/test.groupby('building_id')['meter_reading'].transform('std')","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:19:22.277775Z","iopub.execute_input":"2022-06-29T02:19:22.278156Z","iopub.status.idle":"2022-06-29T02:19:22.283209Z","shell.execute_reply.started":"2022-06-29T02:19:22.278092Z","shell.execute_reply":"2022-06-29T02:19:22.282246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train['meter_reading^2'] = train['meter_reading']**2\n#test['meter_reading^2'] = test['meter_reading']**2","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:19:22.284647Z","iopub.execute_input":"2022-06-29T02:19:22.285007Z","iopub.status.idle":"2022-06-29T02:19:22.296447Z","shell.execute_reply.started":"2022-06-29T02:19:22.284959Z","shell.execute_reply":"2022-06-29T02:19:22.295735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_feature = pd.read_csv('../input/energy-anomaly-detection/train_features.csv')\nfor shift_hours in tqdm([-1,1,-24,24,-7*24,7*24]):\n meter_reading_shift = train_feature[['building_id', 'timestamp', 'meter_reading']]\n meter_reading_shift['timestamp'] = pd.to_datetime(meter_reading_shift['timestamp']) + datetime.timedelta(hours=shift_hours)\n meter_reading_shift['timestamp'] = meter_reading_shift['timestamp'].astype('str')\n meter_reading_shift = meter_reading_shift.rename(columns={'meter_reading':'lag_value_'+str(shift_hours)})\n train_feature = train_feature.merge(meter_reading_shift, on=['building_id', 'timestamp'], how='left')\n \n train['lag_value_'+str(shift_hours)] = train_feature['lag_value_'+str(shift_hours)]-train_feature['meter_reading']\n \ntest_feature = pd.read_csv('../input/energy-anomaly-detection/test_features.csv')\nfor shift_hours in tqdm([-1,1,-24,24,-7*24,7*24]):\n meter_reading_shift = test_feature[['building_id', 'timestamp', 'meter_reading']]\n meter_reading_shift['timestamp'] = pd.to_datetime(meter_reading_shift['timestamp']) + datetime.timedelta(hours=shift_hours)\n meter_reading_shift['timestamp'] = meter_reading_shift['timestamp'].astype('str')\n meter_reading_shift = meter_reading_shift.rename(columns={'meter_reading':'lag_value_'+str(shift_hours)})\n test_feature = test_feature.merge(meter_reading_shift, on=['building_id', 'timestamp'], how='left')\n \n test['lag_value_'+str(shift_hours)] = test_feature['lag_value_'+str(shift_hours)]-test_feature['meter_reading']","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:19:22.29786Z","iopub.execute_input":"2022-06-29T02:19:22.298708Z","iopub.status.idle":"2022-06-29T02:21:37.393176Z","shell.execute_reply.started":"2022-06-29T02:19:22.298654Z","shell.execute_reply":"2022-06-29T02:21:37.391916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.describe()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:21:37.395202Z","iopub.execute_input":"2022-06-29T02:21:37.395746Z","iopub.status.idle":"2022-06-29T02:21:40.311324Z","shell.execute_reply.started":"2022-06-29T02:21:37.395695Z","shell.execute_reply":"2022-06-29T02:21:40.310192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#train['is_one'] = (train['meter_reading']==1).astype('int')\n#test['is_one'] = (test['meter_reading']==1).astype('int')","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:21:40.313046Z","iopub.execute_input":"2022-06-29T02:21:40.313357Z","iopub.status.idle":"2022-06-29T02:21:40.317767Z","shell.execute_reply.started":"2022-06-29T02:21:40.313311Z","shell.execute_reply":"2022-06-29T02:21:40.3165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.drop(['date','meterReadings_daily_std'],axis=1)\ntest = test.drop(['date','meterReadings_daily_std'],axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:33.097317Z","iopub.execute_input":"2022-06-29T02:27:33.097636Z","iopub.status.idle":"2022-06-29T02:27:33.601449Z","shell.execute_reply.started":"2022-06-29T02:27:33.097604Z","shell.execute_reply":"2022-06-29T02:27:33.600633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train.shape, test.shape)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:36.080744Z","iopub.execute_input":"2022-06-29T02:27:36.081565Z","iopub.status.idle":"2022-06-29T02:27:36.086685Z","shell.execute_reply.started":"2022-06-29T02:27:36.081509Z","shell.execute_reply":"2022-06-29T02:27:36.085922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"neg = train[train['anomaly'] == 0]\npos = train[train['anomaly'] == 1]\n\nprint(neg.shape, pos.shape)\nnegs1 = neg.sample(n = 37296, random_state=10)\nnegs2 = neg.sample(n = 37296, random_state=20)\ndf_eq = pd.concat([negs1, pos, negs2, pos], axis=0)\nprint(df_eq.shape)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:36.324985Z","iopub.execute_input":"2022-06-29T02:27:36.325636Z","iopub.status.idle":"2022-06-29T02:27:36.973455Z","shell.execute_reply.started":"2022-06-29T02:27:36.325594Z","shell.execute_reply":"2022-06-29T02:27:36.972524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = df_eq.drop(['anomaly'], axis=1)\ntarget = df_eq['anomaly']","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:36.975087Z","iopub.execute_input":"2022-06-29T02:27:36.975309Z","iopub.status.idle":"2022-06-29T02:27:37.050592Z","shell.execute_reply.started":"2022-06-29T02:27:36.975282Z","shell.execute_reply":"2022-06-29T02:27:37.049769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features.columns","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:37.052676Z","iopub.execute_input":"2022-06-29T02:27:37.05295Z","iopub.status.idle":"2022-06-29T02:27:37.059858Z","shell.execute_reply.started":"2022-06-29T02:27:37.052916Z","shell.execute_reply":"2022-06-29T02:27:37.058811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = features[train_feature['building_id']%5<4]\nX_val = features[train_feature['building_id']%5==4]\nY_train = target[train_feature['building_id']%5<4]\nY_val = target[train_feature['building_id']%5==4]\n\nprint(X_train.shape, X_val.shape)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:37.423631Z","iopub.execute_input":"2022-06-29T02:27:37.423932Z","iopub.status.idle":"2022-06-29T02:27:37.591878Z","shell.execute_reply.started":"2022-06-29T02:27:37.423899Z","shell.execute_reply":"2022-06-29T02:27:37.590889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = StandardScaler()\nX_train = scaler.fit_transform(X_train)\nX_val = scaler.transform(X_val)\ntest1 = scaler.transform(test)","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:37.806824Z","iopub.execute_input":"2022-06-29T02:27:37.80714Z","iopub.status.idle":"2022-06-29T02:27:39.057644Z","shell.execute_reply.started":"2022-06-29T02:27:37.807109Z","shell.execute_reply":"2022-06-29T02:27:39.05682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models = [XGBClassifier(), lgb.LGBMClassifier()]\n\nfor i in range(2):\n models[i].fit(X_train, Y_train)\n\n print(f'{models[i]} : ')\n print('Training Accuracy : ', metrics.roc_auc_score(Y_train, models[i].predict_proba(X_train)[:,1]))\n print('Validation Accuracy : ', metrics.roc_auc_score(Y_val, models[i].predict_proba(X_val)[:,1]))\n print()","metadata":{"execution":{"iopub.status.busy":"2022-06-29T02:27:39.059872Z","iopub.execute_input":"2022-06-29T02:27:39.06022Z","iopub.status.idle":"2022-06-29T02:28:23.574879Z","shell.execute_reply.started":"2022-06-29T02:27:39.060174Z","shell.execute_reply":"2022-06-29T02:28:23.573909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As compare to the model which I build by just using the train.csv I achieved around 95% validation accuracy and now 99% which means that the presence of the extra features has certainly helped to improve the performance of the model.