{"cells":[{"metadata":{},"cell_type":"markdown","source":"Kaggle competition submission from TiimTiim. Started out from LGBM Baseline model https://www.kaggle.com/morituri/lgbm-baseline but is quite modified, with additional logic for missing values and 10 fold stratified models for each meter type."},{"metadata":{},"cell_type":"markdown","source":"# Imports"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":false},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport gc\nimport matplotlib.pyplot as plt\n\nimport seaborn as sns\nimport matplotlib.patches as patches\npd.set_option('max_columns', 150)\nfrom datetime import timedelta\nimport random\n\nfrom sklearn.model_selection import train_test_split, KFold, StratifiedKFold\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.experimental import enable_iterative_imputer \nfrom sklearn.impute import IterativeImputer","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training data"},{"metadata":{},"cell_type":"markdown","source":"Load the data reducing its size"},{"metadata":{"trusted":false},"cell_type":"code","source":"metadata_dtype = {'site_id':\"uint8\",'building_id':'uint16','square_feet':'float32','year_built':'float32','floor_count':\"float16\"}\nweather_dtype = {\"site_id\":\"uint8\",'air_temperature':\"float16\",'cloud_coverage':\"float16\",'dew_temperature':\"float16\",'precip_depth_1_hr':\"float16\",\n                 'sea_level_pressure':\"float32\",'wind_direction':\"float16\",'wind_speed':\"float16\"}\ntrain_dtype = {'meter':\"uint8\",'building_id':'uint16'}","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\n\nweather_train = pd.read_csv(\"../input/ashrae-energy-prediction/weather_train.csv\", parse_dates=['timestamp'], dtype=weather_dtype)\nmetadata = pd.read_csv(\"../input/ashrae-energy-prediction/building_metadata.csv\", dtype=metadata_dtype)\ntrain = pd.read_csv(\"../input/ashrae-energy-prediction/train.csv\", parse_dates=['timestamp'], dtype=train_dtype)\n\nprint('Size of train_df data', train.shape)\nprint('Size of weather_train_df data', weather_train.shape)\nprint('Size of building_meta_df data', metadata.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Improve data readability"},{"metadata":{"trusted":false},"cell_type":"code","source":"train['meter'].replace({0:\"Electricity\",1:\"ChilledWater\",2:\"Steam\",3:\"HotWater\"},inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Add missing measurements to weather"},{"metadata":{},"cell_type":"markdown","source":"There are some missing measurements alltogether as well. For example 2016.12.31 17:00. I simply interpolate between previous and next measurement. There are few days which have almost no weather data (2016.01.05 for example) and this will at the moment be left as is."},{"metadata":{"trusted":false},"cell_type":"code","source":"def add_missing_weather_times(df):\n    first = df['timestamp'].iloc[0]\n    last = df['timestamp'].iloc[-1]\n    \n    for site in df['site_id'].unique():\n        site_data = df[df['site_id'] == site].sort_values('timestamp')\n        site_timestamps = list(site_data['timestamp'])\n        site_num = len(site_timestamps)\n    \n        i = 0\n        while True:\n            time = first + timedelta(hours=i)\n            if time not in site_timestamps:\n                dic = {'site_id': site, 'timestamp': time}\n                    \n                # Add new rows with NAs, because these will be taken care of later\n                for col in ['air_temperature', 'cloud_coverage', 'dew_temperature', 'precip_depth_1_hr', 'sea_level_pressure', 'wind_direction', 'wind_speed']:\n                        dic[col] = np.nan\n\n                df.loc[df.index[-1] + 1] = dic\n\n            if time == last:\n                break\n\n            i += 1\n\n        del site_data, site_timestamps, site_num\n        gc.collect()\n                \n    df.sort_values(['site_id', 'timestamp'], inplace=True)\n    df.reset_index(inplace=True, drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"print(len(weather_train))\nadd_missing_weather_times(weather_train)\nprint(len(weather_train))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Engineer new features I\n\nDay, Month, Hour, Weekend"},{"metadata":{"trusted":false},"cell_type":"code","source":"def create_date_features(df):\n    df['Month'] = df['timestamp'].dt.month.astype(\"uint8\")\n    df['DayOfMonth'] = df['timestamp'].dt.day.astype(\"uint8\")\n    df['DayOfWeek'] = df['timestamp'].dt.dayofweek.astype(\"uint8\")\n    df['Hour'] = df['timestamp'].dt.hour.astype(\"uint8\")\n    df['Weekend'] = ((df['DayOfWeek'] == 6) | (df['DayOfWeek'] == 5)).astype('uint8')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"create_date_features(train)\ncreate_date_features(weather_train) # I'll also add these features to weather, which should help imputer","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"weather_train.drop(['Weekend', 'DayOfWeek', 'DayOfMonth', 'precip_depth_1_hr'], axis=1, inplace=True)\nweather_train","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Deal with NaNs"},{"metadata":{},"cell_type":"markdown","source":"If there is single NaN in weather data, it seems reasonable to linearly interpolate between previous and next measurement."},{"metadata":{"trusted":false},"cell_type":"code","source":"def replace_single_weather_nans(df):\n    for column in ['air_temperature', 'cloud_coverage', 'dew_temperature', 'sea_level_pressure', 'wind_speed']:\n        nan_index = df.index[df[column].isnull()]\n        \n        for i in range(0, len(nan_index)):\n            index = nan_index[i]\n            if index in [0, len(df)]:\n                continue\n            if  np.isfinite(df[column][index - 1]) and np.isfinite(df[column][index + 1]):\n                if df['site_id'][index-1] == df['site_id'][index] == df['site_id'][index + 1]:\n                    df[column][index] = (df[column][index-1] + df[column][index+1])/2","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"print(sum(weather_train.isna().sum()))\nreplace_single_weather_nans(weather_train)\nprint(sum(weather_train.isna().sum()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For other nan-s let's use sklearn's IterativeImputer with BayesianRidge estimator, to predict likely values for other nan-s."},{"metadata":{"trusted":false},"cell_type":"code","source":"imputer = IterativeImputer()\nimputer.fit(weather_train.drop('timestamp', axis=1))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def replace_other_weather_nans(df):\n    col_names = list(df.columns)\n    col_names.remove('timestamp')\n    # Must remove timestamp, since it doesn't work with imputer\n    weather_temp = pd.DataFrame(imputer.transform(df[col_names]), columns=col_names)\n    weather_temp.insert(1, 'timestamp', df['timestamp']) # Reinsert timestamps\n\n    return weather_temp.drop(['Month', 'Hour'], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"weather_train = replace_other_weather_nans(weather_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Map wind direction\n"},{"metadata":{"trusted":false},"cell_type":"code","source":"def map_wind_direction(df):\n    N_idx = (0 < df['wind_direction']) & ((315 < df['wind_direction']) | (df['wind_direction'] <= 45))\n    E_idx = (df['wind_direction'] > 45) & (df['wind_direction'] <= 135)\n    S_idx = (df['wind_direction'] > 135) & (df['wind_direction'] <= 225)\n    W_idx = (df['wind_direction'] > 225) & (df['wind_direction'] <= 315)\n    \n    df['wind_direction'][N_idx] = 1\n    df['wind_direction'][E_idx] = 2\n    df['wind_direction'][S_idx] = 3\n    df['wind_direction'][W_idx] = 4\n    \n    df['wind_direction'].astype('uint8')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"map_wind_direction(weather_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Drop features"},{"metadata":{},"cell_type":"markdown","source":"Drop some columns based on EDA"},{"metadata":{"trusted":false},"cell_type":"code","source":"# Dropping floor_count variable as it has 75% missing values\nmetadata.drop('floor_count',axis=1,inplace=True)\nmetadata.drop('year_built',axis=1,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Convert target to log scale"},{"metadata":{"trusted":false},"cell_type":"code","source":"train['meter_reading'] = np.log1p(train['meter_reading'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Preprocess metadata \n"},{"metadata":{"trusted":false},"cell_type":"code","source":"metadata['primary_use'].replace({\"Healthcare\":\"Other\",\"Parking\":\"Other\",\"Warehouse/storage\":\"Other\",\"Manufacturing/industrial\":\"Other\",\n                                \"Retail\":\"Other\",\"Services\":\"Other\",\"Technology/science\":\"Other\",\"Food sales and service\":\"Other\",\n                                \"Utility\":\"Other\",\"Religious worship\":\"Other\"},inplace=True)\nmetadata['square_feet'] = np.log1p(metadata['square_feet'])\nmetadata['square_feet'] = metadata['square_feet'].astype('float16') #Save space\n\n\n#metadata['year_built'].fillna(-999, inplace=True)\n#metadata['year_built'] = metadata['year_built'].astype('int16')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge data"},{"metadata":{"trusted":false},"cell_type":"code","source":"train = pd.merge(train,metadata,on='building_id',how='left')\nprint (\"Training Data+Metadata Shape {}\".format(train.shape))\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"train = pd.merge(train,weather_train,on=['site_id','timestamp'],how='left')\nprint (\"Training Data+Metadata+Weather Shape {}\".format(train.shape))\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"del weather_train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Drop some training data"},{"metadata":{"trusted":false},"cell_type":"code","source":"# Drop nonsense entries\n# As per the discussion in the following thread, https://www.kaggle.com/c/ashrae-energy-prediction/discussion/117083, there is some discrepancy in the meter_readings for different ste_id's and buildings. It makes sense to delete them\nidx_to_drop = list((train[(train['site_id'] == 0) & (train['timestamp'] < \"2016-05-21 00:00:00\")]).index)\nprint (len(idx_to_drop))\ntrain.drop(idx_to_drop,axis='rows',inplace=True)\n\n# dropping all the electricity meter readings that are 0, after considering them as anomalies.\nidx_to_drop = list(train[(train['meter'] == \"Electricity\") & (train['meter_reading'] == 0)].index)\nprint(len(idx_to_drop))\ntrain.drop(idx_to_drop,axis='rows',inplace=True)\n\n# Drop outliers from training data. Following https://www.kaggle.com/juanmah/ashrae-outliers\nidx_to_drop = list(train[(train.building_id == 1099) | \n                         (train.building_id == 799) | \n                         (train.building_id == 1088) | \n                         (train.building_id == 778) | \n                         (train.building_id == 1168) | \n                         (train.building_id == 1021)].index)\nprint(len(idx_to_drop))\ntrain.drop(idx_to_drop, axis='rows', inplace=True)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Encode features"},{"metadata":{"trusted":false},"cell_type":"code","source":"train.drop('timestamp',axis=1,inplace=True)\n\nle = LabelEncoder()\ntrain['meter']= le.fit_transform(train['meter']).astype(\"uint8\")\ntrain['primary_use']= le.fit_transform(train['primary_use']).astype(\"uint8\")\n\nprint (train.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Drop correlated variables"},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\n# Let's check the correlation between the variables and eliminate the one's that have high correlation\n# Threshold for removing correlated variables\nthreshold = 0.9\n\n# Absolute value correlation matrix\ncorr_matrix = train.corr().abs()\n# Upper triangle of correlations\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))\n\n# Select columns with correlations above threshold\nto_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n\ndel corr_matrix, upper\ngc.collect()\n\nprint('There are %d columns to remove.' % (len(to_drop)))\nprint (\"Following columns can be dropped {}\".format(to_drop))\n\ndef drop_correlated_features(df):\n    df.drop(to_drop,axis=1,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"drop_correlated_features(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\ny = train['meter_reading']\ntrain.drop('meter_reading',axis=1,inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\nx_trains = {meter: train[train['meter'] == meter] for meter in [0, 1, 2, 3]}\ny_trains = {meter: y[train['meter'] == meter] for meter in [0, 1, 2, 3]}","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"del y, train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Model\n\nI'm going to build lightGBM based model with separate model for each meter type and using 10 fold ensamble for each, in order to reduce overfitting."},{"metadata":{"trusted":false},"cell_type":"code","source":"import lightgbm as lgb","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Prediction function"},{"metadata":{"trusted":false},"cell_type":"code","source":"def predict(X, step=1000000):\n    results = {}\n    for meter in [0, 1, 2, 3]:\n        x = X[X.meter == meter]\n        predictions = None\n        for model in models[meter]:\n            preds = []\n            for i in range(0, len(x), step):\n                preds.extend(model.predict(x.iloc[i: min(i+step, len(x)), :], num_iteration=model.best_iteration))\n\n            # Average results\n            if predictions is None:\n                predictions = np.array(preds)/(len(models[meter]))\n            else:\n                predictions += np.array(preds)/(len(models[meter]))\n \n            print('... ', end ='')\n        predictions = np.expm1(predictions) # Back to kWh\n        print()\n        \n        # Create DFs\n        results[meter] = pd.DataFrame(x.index, columns=['row_id'])\n        results[meter]['meter_reading'] = predictions\n        results[meter]['meter_reading'].clip(lower=0,upper=None,inplace=True)\n        \n        del preds, predictions, x\n        gc.collect()\n        \n    # Merge results\n    result = pd.concat([*results.values()])\n    result.sort_values('row_id', inplace=True)\n    \n    return result","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"common_params = {'objective': 'regression',\n                'boosting_type': 'gbdt',\n                'bagging_seed': 11,\n                'metric': 'rmse',\n                'verbosity': -1,\n                'random_state': 47}\n\nparams = {0: {'feature_fraction': 0.75,\n          'bagging_fraction': 0.8,\n          'num_leaves': 300,\n          'max_depth': 15,\n          'learning_rate': 0.14,\n          'min_child_weight': 10,\n          'min_split_gain': 0.005,\n          'reg_alpha': 12.5,\n          'reg_lambda': 7.5,\n          **common_params\n         },\n         1: {'feature_fraction': 0.725,\n          'bagging_fraction': 0.75,\n          'num_leaves': 250, \n          'max_depth': 20,\n          'learning_rate': 0.17,\n          'min_child_weight': 1,\n          'min_split_gain': 0.005,\n          'reg_alpha': 3.,\n          'reg_lambda': 15.,\n          **common_params\n         },\n         2: {'feature_fraction': 0.75,\n          'bagging_fraction': 0.825,\n          'num_leaves': 300,\n          'max_depth': 25,\n          'learning_rate': 0.17,\n          'min_child_weight': 15,\n          'min_split_gain': 0.005,\n          'reg_alpha': 3.,\n          'reg_lambda': 5.,\n          **common_params\n         },\n         3: {'feature_fraction': 0.75,\n          'bagging_fraction': 0.750,\n          'num_leaves': 300,\n          'max_depth': 15,\n          'learning_rate': 0.16,\n          'min_child_weight': 15,\n          'min_split_gain': 0.01,\n          'reg_alpha': 1.,\n          'reg_lambda': 3.,\n          **common_params\n         }}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## KFold model"},{"metadata":{"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":false},"cell_type":"code","source":"%%time\ncategorical_cols = ['building_id','Month','meter','Hour','Weekend','primary_use','DayOfWeek','DayOfMonth', 'wind_direction']\n\nmodels = {0: [], 1: [], 2:[], 3:[]}\nfor meter in [0, 1, 2, 3]:\n    X = x_trains[meter]\n    y = y_trains[meter]\n    \n    X_np = np.array(X)\n    y_np = np.array(y)\n    \n    strats = np.array(pd.cut(y, 50, labels=list(range(50))))\n    \n    kf = StratifiedKFold(n_splits=10, random_state=200)\n    for train_i, test_i in kf.split(X_np, strats):\n        X_train_kf, X_test_kf = X.iloc[train_i], X.iloc[test_i]\n        y_train_kf, y_test_kf = y.iloc[train_i], y.iloc[test_i]\n\n        lgb_train = lgb.Dataset(X_train_kf, y_train_kf, categorical_feature=categorical_cols)\n        lgb_test = lgb.Dataset(X_test_kf, y_test_kf, categorical_feature=categorical_cols)\n\n        reg = lgb.train(params[meter], lgb_train, num_boost_round=500, valid_sets=[lgb_train, lgb_test], \n                    early_stopping_rounds=50, verbose_eval=500)\n        \n        print()\n        del X_train_kf, X_test_kf, y_train_kf, y_test_kf, lgb_train, lgb_test\n        gc.collect()\n    \n        models[meter].append(reg)\n    \n    del X, y, X_np\n    gc.collect()\n    \n    print('\\n------------------------\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Final predictions"},{"metadata":{},"cell_type":"markdown","source":"## Read and modify test data"},{"metadata":{"trusted":false},"cell_type":"code","source":"del x_trains, y_trains","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"test = pd.read_csv(\"../input/ashrae-energy-prediction/test.csv\", parse_dates=['timestamp'], usecols=['building_id','meter','timestamp'], dtype=train_dtype)\nweather_test = pd.read_csv(\"../input/ashrae-energy-prediction/weather_test.csv\", parse_dates=['timestamp'], dtype=weather_dtype)\n\ntest['meter'].replace({0:\"Electricity\",1:\"ChilledWater\",2:\"Steam\",3:\"HotWater\"},inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"add_missing_weather_times(weather_test)\ncreate_date_features(test)\ncreate_date_features(weather_test)\nweather_test.drop(['Weekend', 'DayOfWeek', 'DayOfMonth', 'precip_depth_1_hr'], axis=1, inplace=True)\nreplace_single_weather_nans(weather_test)\nweather_test = replace_other_weather_nans(weather_test)\nmap_wind_direction(weather_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\n# Merge data\ntest = pd.merge(test,metadata,on='building_id',how='left')\nprint (\"Training Data+Metadata Shape {}\".format(test.shape))\ngc.collect()\ntest = pd.merge(test,weather_test,on=['site_id','timestamp'],how='left')\nprint (\"Training Data+Metadata+Weather Shape {}\".format(test.shape))\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del metadata, weather_test\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"test.drop('timestamp',axis=1,inplace=True)\ntest['meter']= le.fit_transform(test['meter']).astype(\"uint8\")\ntest['primary_use']= le.fit_transform(test['primary_use']).astype(\"uint8\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"drop_correlated_features(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Predict"},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\npredictions = predict(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"predictions.to_csv(\"TiimTiim_submission_10.csv\",index=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.5"},"toc-autonumbering":true},"nbformat":4,"nbformat_minor":1}