{"cells":[{"metadata":{},"cell_type":"markdown","source":"Reference: https://www.kaggle.com/hmendonca/shapley-values-for-feature-selection-ashrae by Henrique Mendonça"},{"metadata":{"trusted":false},"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)\n\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.\n\nfrom pathlib import Path\nimport os\nimport os, gc\nimport random\nimport datetime\n\nfrom tqdm import tqdm_notebook as tqdm\n\n# matplotlib and seaborn for plotting\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold, StratifiedKFold\nfrom sklearn import preprocessing\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"## Function to reduce the DF size\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":{},"cell_type":"markdown","source":"## Read Datasets"},{"metadata":{"trusted":false},"cell_type":"code","source":"def dataset_reader():\n    list=['weather_test.csv'\n          ,'weather_train.csv'\n          ,'test.csv'\n          ,'train.csv'\n          ,'building_metadata.csv']\n    input = Path('/kaggle/input/ashrae-energy-prediction')\n    #list= [c for c in os.listdir(input)]\n    \n    wtest = pd.read_csv(input/list[0],parse_dates=['timestamp'])\n    wtrain = pd.read_csv(input/list[1],parse_dates=['timestamp'])\n    test = pd.read_csv(input/list[2],parse_dates=['timestamp'])\n    train = pd.read_csv(input/list[3],parse_dates=['timestamp'])\n    bmdata = pd.read_csv(input/list[4])\n\n    train['is_train'] = 1\n    test['is_train'] = 0\n    \n    # Concatenate and Merge\n    full = pd.concat([train,test],sort=True,ignore_index = True)\n    mean_mr = train.groupby('building_id').meter_reading.mean().reset_index()\n    bmdata = bmdata.merge(mean_mr,on='building_id',how='left')\n    wfull = pd.concat([wtrain,wtest],sort=True,ignore_index = True)\n    return full,wfull,bmdata\n\nfull,wfull,bmdata = dataset_reader()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Building metadata"},{"metadata":{},"cell_type":"markdown","source":"### primary usage"},{"metadata":{"trusted":false},"cell_type":"code","source":"bmdata.groupby('primary_use').meter_reading.mean().sort_values().plot(kind='bar')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"bmdata.groupby('primary_use').meter_reading.mean().sort_values().reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"primary_use = {0: 'Religious worship',\n  1: 'Warehouse/storage',\n  2: 'Technology/science',\n  3: 'Other',\n  4: 'Retail',\n  5: 'Parking',\n  6: 'Lodging/residential',\n  7: 'Manufacturing/industrial',\n  8: 'Public services',\n  9: 'Food sales and service',\n  10: 'Entertainment/public assembly',\n  11: 'Utility',\n  12: 'Office',\n  13: 'Healthcare',\n  14: 'Services',\n  15: 'Education'}\ninv_map = {v: k for k, v in primary_use.items()}\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def primary_use_encoding(x):\n    \n    for use in inv_map.keys():\n        if use == x:\n            return inv_map[use]\n        \nbmdata['pu_label'] = bmdata['primary_use'].apply(primary_use_encoding)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"###  square_feet (log-scale)"},{"metadata":{"trusted":false},"cell_type":"code","source":"bmdata['log_sqf'] = np.log(bmdata.square_feet)\n#from scipy import stats\n#bmdata['bcx_sqf'] = stats.boxcox(bmdata.square_feet)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"'''\nfull['diff_pp'] = full.loc[(~full.precip_depth_1_hr.isnull())&(full.dew_temperature.isnull() | full.air_temperature.isnull())].precip_depth_1_hr.apply(pp_encoding)\nfull['diff_cd'] = full.loc[(~full.cloud_coverage.isnull())&(full.dew_temperature.isnull() | full.air_temperature.isnull())].cloud_coverage.apply(cd_encoding)\ndef pp_encoding(x):\n    \n    if np.isnan(x):\n        return x\n    \n    else:\n        for diff in diff_precip['mean'].keys():\n\n            if diff == x:\n                return round(float(np.random.normal(diff_precip['mean'][diff], diff_precip['std'][diff], 1)),2)\n        \n    return x\ndef cd_encoding(x):\n    \n    if np.isnan(x):\n        return x\n    else:\n        for diff in diff_cloud['mean'].keys():\n\n            if diff == x:\n                return round(float(np.random.normal(diff_cloud['mean'][diff], diff_cloud['std'][diff], 1)),2)\n    return x\n    \n    full['est_dew_p'] = full.air_temperature.loc[~full.diff_pp.isnull()] - full.diff_pp\nfull['est_air_p'] = full.dew_temperature.loc[~full.diff_pp.isnull()] + full.diff_pp\nfull['est_dew_c'] = full.air_temperature.loc[~full.diff_cd.isnull()] - full.diff_cd\nfull['est_air_c'] = full.dew_temperature.loc[~full.diff_cd.isnull()] + full.diff_cd\n\n#precipitation first\nfull.air_temperature.fillna(full.est_air_p, inplace=True)\nfull.dew_temperature.fillna(full.est_dew_p, inplace=True)\n#cloud next\nfull.air_temperature.fillna(full.est_air_c, inplace=True)\nfull.dew_temperature.fillna(full.est_dew_c, inplace=True)\nfull['diff_pp'] = full.loc[(~full.precip_depth_1_hr.isnull())\n                           &(full.dew_temperature.isnull() | full.air_temperature.isnull())].precip_depth_1_hr.apply(pp_encoding)\nfull['diff_cd'] = full.loc[(~full.cloud_coverage.isnull())\n                           &(full.dew_temperature.isnull() | full.air_temperature.isnull())].cloud_coverage.apply(cd_encoding)\n\nfull['est_dew_p'] = full.air_temperature.loc[~full.diff_pp.isnull()] - full.diff_pp\nfull['est_air_p'] = full.dew_temperature.loc[~full.diff_pp.isnull()] + full.diff_pp\nfull['est_dew_c'] = full.air_temperature.loc[~full.diff_cd.isnull()] - full.diff_cd\nfull['est_air_c'] = full.dew_temperature.loc[~full.diff_cd.isnull()] + full.diff_cd\n\n#precipitation first\nfull.air_temperature.fillna(full.est_air_p, inplace=True)\nfull.dew_temperature.fillna(full.est_dew_p, inplace=True)\n#cloud next\nfull.air_temperature.fillna(full.est_air_c, inplace=True)\nfull.dew_temperature.fillna(full.est_dew_c, inplace=True)\n'''\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### floor_count \n(mean filling// if NaN , 1)"},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"bmdata.groupby('primary_use').floor_count.mean().apply(np.ceil).plot(kind='bar')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"floor_avg = bmdata.groupby('primary_use').floor_count.mean().apply(np.ceil).fillna(1).to_dict()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"floor_avg = bmdata.groupby('primary_use').floor_count.mean().apply(np.ceil).fillna(1).to_dict()\ndef floor_encoding(x):\n    \n    if pd.isna(x):\n        return np.nan\n    else:\n        for floor in floor_avg.keys():\n            if floor in x:\n                return floor_avg[floor]\n    return np.nan\n\nnew = bmdata.loc[bmdata.floor_count.isnull()].primary_use.apply(floor_encoding)\nbmdata['floor_count'].fillna(new, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"from sklearn import preprocessing\nencoder = LabelEncoder()\nbmdata['primary_use'] = encoder.fit_transform(bmdata['primary_use'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### year_built"},{"metadata":{"trusted":false},"cell_type":"code","source":"bmdata['year_built']=bmdata['year_built'].fillna(bmdata['year_built'].mean())\nbmdata['year_built']=bmdata['year_built']-1900\nbmdata.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## weather data"},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"#from tqdm import tqdm\n#lists = ['air_temperature','dew_temperature','cloud_coverage','sea_level_pressure','wind_direction','wind_speed','precip_depth_1_hr']\n#size = full.building_id.nunique()\n#for li in lists:\n    #print(li)\n    #for i in tqdm(range(size)):\n        #full[li].update(full.loc[full.building_id==i][li].interpolate(method='pchip',limit_direction='both'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* air_temperature\n* dew_temperature\n* cloud_coverage\n* sea_level_pressure wind_direction\n* wind_speed\n* precip_depth_1_hr \n\nTentative estimation\n($forward\\ fill \\rightarrow backward\\ fill$)"},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"from tqdm import tqdm\nlists = ['air_temperature','dew_temperature','cloud_coverage','sea_level_pressure','wind_direction','wind_speed','precip_depth_1_hr']\nsize = wfull.site_id.nunique()\nfor li in lists:\n    print(li)\n    for i in tqdm(range(size)):\n        wfull[li].update(wfull.loc[wfull.site_id==i][li].interpolate(method='ffill'))\n        wfull[li].update(wfull.loc[wfull.site_id==i][li].interpolate(method='bfill'))\nwfull.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"wfull.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Observe : Some sites are having etirely missing feature."},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"wfull.groupby('site_id')[['air_temperature', 'cloud_coverage', 'dew_temperature',\n       'precip_depth_1_hr', 'sea_level_pressure', 'site_id', 'timestamp',\n       'wind_direction', 'wind_speed']].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"wfull['Month']= wfull.timestamp.dt.month\nwfull['Day']= wfull.timestamp.dt.day\nwfull['Hour'] = wfull.timestamp.dt.hour\nwfull['Weekday'] = wfull.timestamp.dt.weekday","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Mean filling w.r.t Month"},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"for i in tqdm(range(12)):\n    wfull.update(wfull[wfull.Month==i+1].fillna(wfull[wfull.Month==i+1].mean()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Add Relative Humidity Approx\n\n$$100*\\frac{(\\exp\\frac{17.27*dew}{dew+237.3})}{(\\exp\\frac{17.27*air}{air+237.3})}$$"},{"metadata":{"trusted":false},"cell_type":"code","source":"dt = wfull.dew_temperature\nat = wfull.air_temperature\nws = wfull.wind_speed\n#Relative Humidity\nwfull['RH'] = 100*(0.6108*np.exp((17.27*dt)/(dt+237.3)))/(0.6108*np.exp((17.27*at)/(at+237.3)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### rescale some features\n\n$sea\\ level\\ pressure \\rightarrow sea\\ level\\ pressure - 1000$\n\n$wind\\ direction \\rightarrow wind\\ direction \\mod 360$"},{"metadata":{"trusted":false},"cell_type":"code","source":"wfull.sea_level_pressure = wfull.sea_level_pressure - 1000\nwfull.wind_direction = wfull.wind_direction%360\nwfull.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"wfull.groupby('Month')[['air_temperature','dew_temperature','cloud_coverage','sea_level_pressure','wind_direction','wind_speed','precip_depth_1_hr']].mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"import gc\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Merge all datasets"},{"metadata":{"trusted":false},"cell_type":"code","source":"full = reduce_mem_usage(full)\nwfull = reduce_mem_usage(wfull)\nbmdata = reduce_mem_usage(bmdata)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"del bmdata['meter_reading']\nfull = full.merge(bmdata, on='building_id', how='left')\nfull = full.merge(wfull, on=['site_id', 'timestamp'], how='left')\ndel bmdata\ndel wfull\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"full.sample(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"train = full.loc[(full.is_train==1)&(full.building_id<15)]\ntest = full[full.is_train==0]\nprint('done')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"train[\"meter_reading\"]=np.log1p(train[\"meter_reading\"])\nprint('done here')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"target= 'meter_reading'\ndo_not_use = ['meter_reading'\n                 ,'is_train'\n                ,'row_id'\n                ,'square_feet'\n                ,'timestamp'\n                ,'primary_use'\n                ,'random']\n\nfeature_columns = [c for c in full.columns if c not in do_not_use ]\nfeature_columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"gc.collect()\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"## Training(LGBM)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"import lightgbm as lgb\nfolds = 4\nseed = 777\nmodels=[]\nfeature_importance = pd.DataFrame()\nkf = StratifiedKFold(n_splits=folds, shuffle=True, random_state=seed)\nfor train_idx, valid_idx in tqdm(kf.split(train,train['building_id']),total=folds):\n    print(f'Training and predicting for target {target}')\n    Xtr = train[feature_columns].iloc[train_idx]\n    Xv = train[feature_columns].iloc[valid_idx]\n    ytr = train[target].iloc[train_idx].values\n    yv = train[target].iloc[valid_idx].values\n    print('Train_size: ',Xtr.shape[0],'Validation_size: ', ytr.shape[0])\n    \n    dtrain = lgb.Dataset(Xtr, label=ytr)\n    dvalid = lgb.Dataset(Xv, label=yv)\n    \n    params = {\n            'boosting_type': 'gbdt',\n            'objective': 'regression',\n            'metric': {'rmse'},\n            'learning_rate': 0.5,\n            'feature_fraction': 0.8,\n            'bagging_fraction': 0.8,\n            'bagging_freq' : 5\n            }\n    model = lgb.train(params,\n                dtrain,\n                num_boost_round=2000,\n                valid_sets=(dtrain, dvalid),\n               early_stopping_rounds=20,\n               verbose_eval = 20)\n    \n    \n    #feature importance\n    #f_imp = pd.DataFrame()\n    #f_imp['feature'] = feature_columns\n    #f_imp[\"importance\"] = model.feature_importances_\n    #f_imp[\"fold\"] = nfold\n    #nfold += 1\n    #feature_importance = pd.concat([feature_importance, f_imp],axis=0,ignore_index=True)\n    models.append(model)\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Overfitting :D"},{"metadata":{"trusted":false},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfeature_imp = pd.DataFrame(sorted(zip(models[0].feature_importance(), models[0].feature_name()),reverse = True), columns=['Value','Feature'])\nplt.figure(figsize=(10, 5))\nsns.barplot(x=\"Value\", y=\"Feature\", data=feature_imp.sort_values(by=\"Value\", ascending=False))\nplt.title('LightGBM Features (avg over folds)')\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"# split test data into batches\nset_size = len(test)\niterations = 50\nbatch_size = set_size // iterations\n\nprint(set_size, iterations, batch_size)\nassert set_size == iterations * batch_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"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[feature_columns].iloc[pos : pos+batch_size])) for model in models]\n    meter_reading.extend(np.mean(fold_preds, axis=0))\n\nprint(len(meter_reading))\nassert len(meter_reading) == set_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"submission = pd.read_csv('/kaggle/input/ashrae-energy-prediction/sample_submission.csv')\nsubmission['meter_reading'] = np.clip(meter_reading, a_min=0, a_max=None)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)\nsubmission.head(9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### SHAP Explainer(Practice)"},{"metadata":{"trusted":false},"cell_type":"code","source":"#explainer = shap.TreeExplainer(models[0])\n#shap_values = explainer.shap_values(train[feature_columns])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#shap.force_plot(explainer.expected_value,shap_values[0,:] ,train[feature_columns].iloc[0,:], matplotlib=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#shap.summary_plot(shap_values, train[feature_columns], plot_type=\"bar\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n","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.6.6"}},"nbformat":4,"nbformat_minor":1}