{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Meter-specific LGBM for ASHRAE","metadata":{"execution":{"iopub.status.busy":"2021-06-18T12:45:41.771878Z","iopub.execute_input":"2021-06-18T12:45:41.772452Z","iopub.status.idle":"2021-06-18T12:45:41.775498Z","shell.execute_reply.started":"2021-06-18T12:45:41.772418Z","shell.execute_reply":"2021-06-18T12:45:41.774852Z"}}},{"cell_type":"markdown","source":"This is a modified version of [this other kernel](https://www.kaggle.com/michelezoccali/ashrae-energy-prediction-single-model) with meter-specific models rather than a single one, in an attempt to capture effects on consumption specific to each meter type.","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O\nimport os\nimport datetime\nimport warnings\nimport gc\n\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import KFold\nimport lightgbm as lgb\n\nfrom tqdm.notebook import tqdm\n\n# plotting\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:32.200743Z","iopub.execute_input":"2021-06-17T17:19:32.201188Z","iopub.status.idle":"2021-06-17T17:19:34.422475Z","shell.execute_reply.started":"2021-06-17T17:19:32.201094Z","shell.execute_reply":"2021-06-17T17:19:34.421658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = '../input/ashrae-energy-prediction'\n\nfor dirname, _, filenames in os.walk(path):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:34.423572Z","iopub.execute_input":"2021-06-17T17:19:34.423967Z","iopub.status.idle":"2021-06-17T17:19:34.430579Z","shell.execute_reply.started":"2021-06-17T17:19:34.423937Z","shell.execute_reply":"2021-06-17T17:19:34.429519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utilities","metadata":{}},{"cell_type":"markdown","source":"Memory reduction adapted from [this kernel.](https://www.kaggle.com/purist1024/ashrae-simple-data-cleanup-lb-1-08-no-leaks/notebook)","metadata":{}},{"cell_type":"code","source":"def reduce_mem(df):\n    result = df.copy()\n    for col in result.columns:\n        col_data = result[col]\n        dn = col_data.dtype.name\n        if not dn.startswith(\"datetime\"):  # avoid changing datetime\n            if dn == \"object\":  # only object features have low cardinality\n                result[col] = pd.to_numeric(col_data.astype(\"category\").cat.codes, downcast=\"unsigned\")\n            elif dn.startswith(\"int\") | dn.startswith(\"uint\"):\n                if col_data.min() >= 0:\n                    result[col] = pd.to_numeric(col_data, downcast=\"unsigned\")\n                else:\n                    result[col] = pd.to_numeric(col_data, downcast='integer')\n            else:\n                result[col] = pd.to_numeric(col_data, downcast='float')\n    return result","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:34.433039Z","iopub.execute_input":"2021-06-17T17:19:34.433479Z","iopub.status.idle":"2021-06-17T17:19:34.442979Z","shell.execute_reply.started":"2021-06-17T17:19:34.433436Z","shell.execute_reply":"2021-06-17T17:19:34.44187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Routine to add lag features to weather dataset, adapted from [this kernel](https://www.kaggle.com/corochann/ashrae-training-lgbm-by-meter-type/notebook).","metadata":{}},{"cell_type":"code","source":"def add_lag_features(weather_df, window=3):\n    group_df = weather_df.groupby('site_id')\n    cols = ['air_temperature', 'cloud_coverage', 'dew_temperature', 'precip_depth_1_hr']\n    rolled = group_df[cols].rolling(window=window, min_periods=0)\n    lag_mean = rolled.mean().reset_index().astype(np.float16)\n    lag_max = rolled.max().reset_index().astype(np.float16)\n    lag_min = rolled.min().reset_index().astype(np.float16)\n    for col in cols:\n        weather_df[f'{col}_min_lag{window}'] = lag_min[col]\n        weather_df[f'{col}_mean_lag{window}'] = lag_mean[col]\n        weather_df[f'{col}_max_lag{window}'] = lag_max[col]\n    return weather_df","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:34.445135Z","iopub.execute_input":"2021-06-17T17:19:34.445633Z","iopub.status.idle":"2021-06-17T17:19:34.456422Z","shell.execute_reply.started":"2021-06-17T17:19:34.445578Z","shell.execute_reply":"2021-06-17T17:19:34.455284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"def load_data(source='train'):\n    assert source in ['train','test']\n    df = pd.read_csv(f'{path}/{source}.csv', parse_dates=['timestamp'])\n    return reduce_mem(df)\n\ndef load_building():\n    df = pd.read_csv(f'{path}/building_metadata.csv').fillna(-1)\n    return reduce_mem(df)\n\ndef load_weather(source='train', fix_timezone=True, impute=True, add_lag=True):\n    assert source in ['train','test']\n    df = pd.read_csv(f'{path}/weather_{source}.csv', parse_dates=['timestamp'])\n    if fix_timezone:\n        offsets = [5,0,9,6,8,0,6,6,5,7,8,6,0,7,6,6]\n        offset_map = {site: offset for site, offset in enumerate(offsets)}\n        df.timestamp = df.timestamp - pd.to_timedelta(df.site_id.map(offset_map), unit='h')\n    if impute:\n        site_dfs = []\n        for site in df.site_id.unique():\n            if source == 'train':\n                new_idx = pd.date_range(start='2016-1-1', end='2016-12-31-23', freq='H')\n            else:\n                new_idx = pd.date_range(start='2017-1-1', end='2018-12-31-23', freq='H')\n            site_df = df[df.site_id == site].set_index('timestamp').reindex(new_idx)\n            site_df.site_id = site\n            for col in [c for c in site_df.columns if c != 'site_id']:\n                site_df[col] = site_df[col].interpolate(limit_direction='both', method='linear')\n                site_df[col] = site_df[col].fillna(df[col].median())\n            site_dfs.append(site_df)\n        df = pd.concat(site_dfs)\n        df['timestamp'] = df.index\n        df = df.reset_index(drop=True)\n        \n    if add_lag:\n        df = add_lag_features(df, window=3)\n    \n    return reduce_mem(df)\n\ndef merged_dfs(source='train', fix_timezone=True, impute=True, add_lag=True):\n    df = load_data(source=source).merge(load_building(), on='building_id', how='left')\n    df = df.merge(load_weather(source=source, fix_timezone=True, impute=True, add_lag=True),\n                 on=['site_id','timestamp'], how='left')\n    if source == 'train':\n        X = df.drop('meter_reading', axis=1)  \n        y = np.log1p(df.meter_reading)  # log-transform of target\n        return X, y\n    elif source == 'test':\n        return df","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:34.457868Z","iopub.execute_input":"2021-06-17T17:19:34.458308Z","iopub.status.idle":"2021-06-17T17:19:34.476532Z","shell.execute_reply.started":"2021-06-17T17:19:34.458263Z","shell.execute_reply":"2021-06-17T17:19:34.475526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nX_train, y_train = merged_dfs()\nX_train.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:19:34.47771Z","iopub.execute_input":"2021-06-17T17:19:34.478295Z","iopub.status.idle":"2021-06-17T17:20:06.397635Z","shell.execute_reply.started":"2021-06-17T17:19:34.478249Z","shell.execute_reply":"2021-06-17T17:20:06.396648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train.info()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:06.398735Z","iopub.execute_input":"2021-06-17T17:20:06.399005Z","iopub.status.idle":"2021-06-17T17:20:06.413281Z","shell.execute_reply.started":"2021-06-17T17:20:06.398979Z","shell.execute_reply":"2021-06-17T17:20:06.412291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This corresponds to more than 50% improvement in memory usage! Plus, we have already filled time gaps and missing values.","metadata":{}},{"cell_type":"markdown","source":"# Preprocessing","metadata":{}},{"cell_type":"markdown","source":"Let us remove the first 141 days of electrical meter readings at site 0, which are mostly zero or contain anomalous spikes.","metadata":{}},{"cell_type":"code","source":"def _delete_bad_sitezero(X, y):\n    cond = (X.timestamp > '2016-05-20') | (X.site_id != 0) | (X.meter != 0)\n    X = X[cond]\n    y = y.reindex_like(X)\n    return X.reset_index(drop=True), y.reset_index(drop=True)\n\ndef _extract_temporal(X):\n    X['hour'] = X.timestamp.dt.hour\n    X['weekday'] = X.timestamp.dt.weekday\n    # month and year cause overfit, could try other (holiday, business, etc.)\n    return reduce_mem(X)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:06.416002Z","iopub.execute_input":"2021-06-17T17:20:06.416366Z","iopub.status.idle":"2021-06-17T17:20:06.424506Z","shell.execute_reply.started":"2021-06-17T17:20:06.416317Z","shell.execute_reply":"2021-06-17T17:20:06.423619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# preprocessing\nX_train, y_train = _delete_bad_sitezero(X_train, y_train)\nX_train = _extract_temporal(X_train)\n\n# remove timestamp and other unimportant features\nto_drop = ['timestamp','sea_level_pressure','wind_direction','wind_speed']\nX_train.drop(to_drop, axis=1, inplace=True)\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:06.426027Z","iopub.execute_input":"2021-06-17T17:20:06.426326Z","iopub.status.idle":"2021-06-17T17:20:21.142716Z","shell.execute_reply.started":"2021-06-17T17:20:06.426299Z","shell.execute_reply":"2021-06-17T17:20:21.141665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling (one model per meter)","metadata":{}},{"cell_type":"markdown","source":"Let us write a small wrapper for LGBM, to pass to a CV routine.","metadata":{}},{"cell_type":"code","source":"def LGBM_wrapper(Xt, yt, Xv, yv, fold=-1):\n    \n    cat_features = ['building_id','site_id','primary_use','hour','weekday']\n    \n    dset = lgb.Dataset(Xt, label=yt, categorical_feature=cat_features)\n    dset_val = lgb.Dataset(Xv, label=yv, categorical_feature=cat_features)\n    \n    params = {\n        \"objective\": \"regression\",\n        \"boosting\": \"gbdt\",\n        \"num_leaves\": 40,\n        \"learning_rate\": 0.05,\n        \"feature_fraction\": 0.85,\n        \"reg_lambda\": 2,\n        \"metric\": \"rmse\",\n        \"verbose\": -1\n    }\n    \n    print(f'Fold {fold}')\n    \n    # filter some known warnings (open issue at https://github.com/microsoft/LightGBM/issues/3379)\n    with warnings.catch_warnings():\n        warnings.filterwarnings(\"ignore\", \"categorical_column in param dict is overridden\")\n        warnings.filterwarnings(\"ignore\", \"Overriding the parameters from Reference Dataset\")\n        model = lgb.train(params,\n                         train_set=dset,\n                         num_boost_round=1000,\n                         valid_sets=[dset, dset_val],\n                         verbose_eval=200,\n                         early_stopping_rounds=100,\n                         categorical_feature=cat_features)\n    \n    oof = model.predict(Xv, num_iteration=model.best_iteration)\n    score = np.sqrt(mean_squared_error(yv, oof))\n    print(f'Fold {fold} RMSLE: {score}\\n')\n    return model, oof, score","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:21.144175Z","iopub.execute_input":"2021-06-17T17:20:21.144511Z","iopub.status.idle":"2021-06-17T17:20:21.153888Z","shell.execute_reply.started":"2021-06-17T17:20:21.144478Z","shell.execute_reply":"2021-06-17T17:20:21.152594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's train a model for each meter, to maximise meter-specific learning (at the expense for interaction effects between meter types). \n\nFor each model we perform k-fold CV, without shuffling as this is a time series. An alternative would be to do a single train/validation split, possibly with a gap to mimic training/private split.","metadata":{}},{"cell_type":"code","source":"def CV_per_meter(n_splits=3):  \n    \n    oof_total = np.zeros(X_train.shape[0])\n\n    kf = KFold(n_splits=n_splits, shuffle=False)\n    models = {}\n    scores = {}\n\n    for meter in tqdm(X_train.meter.sort_values().unique()):\n\n        print(f'TRAINING MODEL FOR METER {meter}')\n        X_train_meter = X_train[X_train.meter==meter].drop(columns='meter')\n        y_train_meter = y_train.reindex_like(X_train_meter)\n\n        models[meter] = []\n        scores[meter] = []\n\n        for fold, (train_idx, val_idx) in enumerate(kf.split(X_train_meter), start=1):\n            Xt, yt = X_train_meter.iloc[train_idx], y_train_meter.iloc[train_idx]\n            Xv, yv = X_train_meter.iloc[val_idx], y_train_meter.iloc[val_idx]\n            model, oof, score = LGBM_wrapper(Xt, yt, Xv, yv, fold)\n\n            models[meter].append(model)\n            scores[meter].append(score)\n            oof_total[X_train_meter.index[val_idx]] = oof\n\n        rmsle_meter = np.sqrt(mean_squared_error(y_train_meter, oof_total[y_train_meter.index]))\n        \n        print(f'> Mean RMSLE for meter {meter}: {np.mean(scores[meter])}, std: {np.std(scores[meter])}')\n        print(f'> OOF RMSLE for meter {meter}: {rmsle_meter}\\n')\n    print('Training completed.')\n    return models, scores, oof_total","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:21.155453Z","iopub.execute_input":"2021-06-17T17:20:21.155878Z","iopub.status.idle":"2021-06-17T17:20:21.175851Z","shell.execute_reply.started":"2021-06-17T17:20:21.155832Z","shell.execute_reply":"2021-06-17T17:20:21.174187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_splits = 3\nmodels, scores, oof_total = CV_per_meter(n_splits=n_splits)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:20:21.177421Z","iopub.execute_input":"2021-06-17T17:20:21.177902Z","iopub.status.idle":"2021-06-17T17:32:54.730851Z","shell.execute_reply.started":"2021-06-17T17:20:21.177853Z","shell.execute_reply":"2021-06-17T17:32:54.730112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'TOTAL OOF RMSLE: {np.sqrt(mean_squared_error(y_train, oof_total))}')","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:32:54.737188Z","iopub.execute_input":"2021-06-17T17:32:54.737577Z","iopub.status.idle":"2021-06-17T17:32:54.903847Z","shell.execute_reply.started":"2021-06-17T17:32:54.737546Z","shell.execute_reply":"2021-06-17T17:32:54.903104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature importance","metadata":{}},{"cell_type":"markdown","source":"Let's see the average feature importance across models (per meter type). We can use this to retroactively drop further superfluous features during preprocessing.","metadata":{}},{"cell_type":"code","source":"def plot_importance(meter):    \n    importance = pd.DataFrame([model.feature_importance() for model in models[meter]],\n                              columns=X_train.drop(columns='meter').columns,\n                              index=[f'Fold {i}' for i in range(1, n_splits + 1)])\n    importance = importance.T\n    importance['Average importance'] = importance.mean(axis=1)\n    importance = importance.sort_values(by='Average importance', ascending=False)\n\n    plt.figure(figsize=(10,7))\n    sns.barplot(x='Average importance', y=importance.index, data=importance);","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:32:54.905016Z","iopub.execute_input":"2021-06-17T17:32:54.905451Z","iopub.status.idle":"2021-06-17T17:32:54.912645Z","shell.execute_reply.started":"2021-06-17T17:32:54.905403Z","shell.execute_reply":"2021-06-17T17:32:54.911832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot it for e.g. meter 3:","metadata":{}},{"cell_type":"code","source":"plot_importance(3)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:32:54.913835Z","iopub.execute_input":"2021-06-17T17:32:54.914278Z","iopub.status.idle":"2021-06-17T17:32:56.760822Z","shell.execute_reply.started":"2021-06-17T17:32:54.914242Z","shell.execute_reply":"2021-06-17T17:32:56.759684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del X_train, y_train\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:32:56.762262Z","iopub.execute_input":"2021-06-17T17:32:56.762583Z","iopub.status.idle":"2021-06-17T17:32:56.970924Z","shell.execute_reply.started":"2021-06-17T17:32:56.76255Z","shell.execute_reply":"2021-06-17T17:32:56.969922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test set, inference and submission","metadata":{}},{"cell_type":"code","source":"%%time\nX_test = merged_dfs('test')","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:32:56.97255Z","iopub.execute_input":"2021-06-17T17:32:56.973101Z","iopub.status.idle":"2021-06-17T17:34:10.472326Z","shell.execute_reply.started":"2021-06-17T17:32:56.973004Z","shell.execute_reply":"2021-06-17T17:34:10.470693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"row_ids = X_test.row_id # for submission file\nX_test = _extract_temporal(X_test)\nX_test.drop(columns=['row_id','timestamp']+to_drop, inplace=True)\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:34:10.474738Z","iopub.execute_input":"2021-06-17T17:34:10.475129Z","iopub.status.idle":"2021-06-17T17:34:36.952923Z","shell.execute_reply.started":"2021-06-17T17:34:10.475085Z","shell.execute_reply":"2021-06-17T17:34:36.951774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test.info()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:34:36.954413Z","iopub.execute_input":"2021-06-17T17:34:36.954745Z","iopub.status.idle":"2021-06-17T17:34:36.971421Z","shell.execute_reply.started":"2021-06-17T17:34:36.954711Z","shell.execute_reply":"2021-06-17T17:34:36.970129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For each meter, we compute the predictions on the test set for each CV model, then average the results. We split the computation in batches, to keep memory usage within the limits. Naturally, we transform the predictions back into linear space with the inverse of the log-transform.","metadata":{}},{"cell_type":"code","source":"def get_preds(n_iterations=5):   \n    \n    all_preds = np.zeros(X_test.shape[0])\n    \n    for meter in range(4):\n\n        X_test_meter = X_test[X_test.meter == meter].drop(columns='meter')\n        batch_size = len(X_test_meter) // n_iterations\n        preds = []\n\n        print(f'Getting predictions for meter {meter}...')\n        for i in tqdm(range(n_iterations)):\n            start = i * batch_size\n            fold_preds = [np.expm1(model.predict(X_test_meter.iloc[start:start + batch_size], \n                             num_iteration=model.best_iteration)) for model in models[meter]]\n            preds.extend(np.mean(fold_preds, axis=0))\n\n        all_preds[X_test_meter.index] = preds\n        \n    return all_preds\n\nall_preds = get_preds()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:34:36.972999Z","iopub.execute_input":"2021-06-17T17:34:36.973382Z","iopub.status.idle":"2021-06-17T17:46:25.451016Z","shell.execute_reply.started":"2021-06-17T17:34:36.973345Z","shell.execute_reply":"2021-06-17T17:46:25.449723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del X_test\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:46:25.452695Z","iopub.execute_input":"2021-06-17T17:46:25.453015Z","iopub.status.idle":"2021-06-17T17:46:25.653515Z","shell.execute_reply.started":"2021-06-17T17:46:25.452967Z","shell.execute_reply":"2021-06-17T17:46:25.652596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, save for submission and hope for the best.","metadata":{}},{"cell_type":"code","source":"submission = pd.DataFrame({'row_id':row_ids, 'meter_reading':np.clip(all_preds, 0, a_max=None)})\nsubmission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:46:25.670152Z","iopub.execute_input":"2021-06-17T17:46:25.670558Z","iopub.status.idle":"2021-06-17T17:48:58.183236Z","shell.execute_reply.started":"2021-06-17T17:46:25.670525Z","shell.execute_reply":"2021-06-17T17:48:58.182139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can check that the distribution of the predictions looks reasonable:","metadata":{}},{"cell_type":"code","source":"sns.displot(np.log1p(submission.meter_reading));","metadata":{"execution":{"iopub.status.busy":"2021-06-17T17:48:58.186233Z","iopub.execute_input":"2021-06-17T17:48:58.186699Z","iopub.status.idle":"2021-06-17T17:49:37.777978Z","shell.execute_reply.started":"2021-06-17T17:48:58.186646Z","shell.execute_reply":"2021-06-17T17:49:37.776927Z"},"trusted":true},"execution_count":null,"outputs":[]}]}