{"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":"# Single LGBM for ASHRAE","metadata":{}},{"cell_type":"markdown","source":"This is a simple kernel which does the following:\n- **preprocessing**: during data import, minimise memory usage (typecasting) given the size of the dataset and the instance's RAM, correct weather data (timezone alignment, imputation) and log-transform the target.\n- **cleaning**: drops bad readings in site 0.\n- **FE**: basic temporal and lag features, to avoid overfitting.\n- **modelling**: trains a single LGBM with (unshuffled) k-fold CV (plus a random forest for reference).\n- **inference**: bagging of model's versions trained on each CV fold.","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\nfrom sklearn.ensemble import RandomForestRegressor\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-23T06:49:07.763052Z","iopub.execute_input":"2021-06-23T06:49:07.763494Z","iopub.status.idle":"2021-06-23T06:49:10.233203Z","shell.execute_reply.started":"2021-06-23T06:49:07.763399Z","shell.execute_reply":"2021-06-23T06:49:10.232274Z"},"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-23T06:49:10.234664Z","iopub.execute_input":"2021-06-23T06:49:10.234972Z","iopub.status.idle":"2021-06-23T06:49:10.240856Z","shell.execute_reply.started":"2021-06-23T06:49:10.234942Z","shell.execute_reply":"2021-06-23T06:49:10.239935Z"},"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\"):\n            if dn == \"object\":  # only object feature has 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-23T06:49:10.312536Z","iopub.execute_input":"2021-06-23T06:49:10.312934Z","iopub.status.idle":"2021-06-23T06:49:10.321523Z","shell.execute_reply.started":"2021-06-23T06:49:10.3129Z","shell.execute_reply":"2021-06-23T06:49:10.320418Z"},"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-23T06:49:10.323057Z","iopub.execute_input":"2021-06-23T06:49:10.323535Z","iopub.status.idle":"2021-06-23T06:49:10.335487Z","shell.execute_reply.started":"2021-06-23T06:49:10.323504Z","shell.execute_reply":"2021-06-23T06:49:10.334416Z"},"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-23T07:03:00.91485Z","iopub.execute_input":"2021-06-23T07:03:00.915347Z","iopub.status.idle":"2021-06-23T07:03:00.933342Z","shell.execute_reply.started":"2021-06-23T07:03:00.915308Z","shell.execute_reply":"2021-06-23T07:03:00.932276Z"},"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-23T06:49:10.355554Z","iopub.execute_input":"2021-06-23T06:49:10.355833Z","iopub.status.idle":"2021-06-23T06:49:43.242185Z","shell.execute_reply.started":"2021-06-23T06:49:10.355807Z","shell.execute_reply":"2021-06-23T06:49:43.241357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train.info()","metadata":{"execution":{"iopub.status.busy":"2021-06-23T06:49:43.243329Z","iopub.execute_input":"2021-06-23T06:49:43.243612Z","iopub.status.idle":"2021-06-23T06:49:43.261707Z","shell.execute_reply.started":"2021-06-23T06:49:43.243585Z","shell.execute_reply":"2021-06-23T06:49:43.260658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This corresponds to roughly a 50% improvement in memory usage, despite having added lag features.","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. This is the type of outlier which causes the most trouble and is comparatively easier to remove. We also extract some basic temporal features.","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-22T12:29:04.856943Z","iopub.execute_input":"2021-06-22T12:29:04.857598Z","iopub.status.idle":"2021-06-22T12:29:04.871255Z","shell.execute_reply.started":"2021-06-22T12:29:04.857557Z","shell.execute_reply":"2021-06-22T12:29:04.869752Z"},"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-22T12:29:04.873169Z","iopub.execute_input":"2021-06-22T12:29:04.873856Z","iopub.status.idle":"2021-06-22T12:29:23.062129Z","shell.execute_reply.started":"2021-06-22T12:29:04.873804Z","shell.execute_reply":"2021-06-22T12:29:23.060813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"markdown","source":"Before training the LGBM, let us train a \"baseline\" random forest for comparison. We write small wrappers for each, to be passed to a CV routine.","metadata":{}},{"cell_type":"code","source":"def RF_wrapper(Xt, yt, Xv, yv, fold=-1):\n    \n    model = RandomForestRegressor(n_jobs=-1, n_estimators=40,\n                              max_samples=200000, max_features=0.5,\n                              min_samples_leaf=5, oob_score=False).fit(Xt, yt)\n    print(f'Training fold {fold}...')\n    \n    score_train = np.sqrt(mean_squared_error(model.predict(Xt), yt))\n    oof = model.predict(Xv)\n    score = np.sqrt(mean_squared_error(oof, yv))\n    print(f'Fold {fold}: training RMSLE: {score_train},   validation RMSLE: {score}\\n')\n    return model, oof, score","metadata":{"execution":{"iopub.status.busy":"2021-06-22T13:03:02.158129Z","iopub.execute_input":"2021-06-22T13:03:02.158987Z","iopub.status.idle":"2021-06-22T13:03:02.166907Z","shell.execute_reply.started":"2021-06-22T13:03:02.158912Z","shell.execute_reply":"2021-06-22T13:03:02.165706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def LGBM_wrapper(Xt, yt, Xv, yv, fold=-1):\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\": 500,\n        \"learning_rate\": 0.04,\n        \"feature_fraction\": 0.7,\n        \"subsample\": 0.4,\n        \"metric\": \"rmse\",\n        \"seed\": 42,\n        \"n_jobs\": -1,\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} validation RMSLE: {score}\\n')\n    return model, oof, score","metadata":{"execution":{"iopub.status.busy":"2021-06-22T12:54:29.992932Z","iopub.execute_input":"2021-06-22T12:54:29.993357Z","iopub.status.idle":"2021-06-22T12:54:30.003739Z","shell.execute_reply.started":"2021-06-22T12:54:29.993321Z","shell.execute_reply":"2021-06-22T12:54:30.002475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let us 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. Otherwise, one could try something like Time-series split CV.","metadata":{}},{"cell_type":"code","source":"def perform_CV(wrapper, n_splits=3):\n    \n    kf = KFold(n_splits=n_splits, shuffle=False)\n\n    models = []\n    scores = []\n    oof_total = np.zeros(X_train.shape[0])\n\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X_train), start=1):\n        Xt, yt = X_train.iloc[train_idx], y_train[train_idx]\n        Xv, yv = X_train.iloc[val_idx], y_train[val_idx]\n        model, oof, score = wrapper(Xt, yt, Xv, yv, fold)\n\n        models.append(model)\n        scores.append(score)\n        oof_total[val_idx] = oof\n\n    print('Training completed.')\n    print(f'> Mean RMSLE across folds: {np.mean(scores)}, std: {np.std(scores)}')\n    print(f'> OOF RMSLE: {np.sqrt(mean_squared_error(y_train, oof_total))}')\n    return models, scores, oof_total","metadata":{"execution":{"iopub.status.busy":"2021-06-22T12:59:11.439303Z","iopub.execute_input":"2021-06-22T12:59:11.439806Z","iopub.status.idle":"2021-06-22T12:59:11.449281Z","shell.execute_reply.started":"2021-06-22T12:59:11.439764Z","shell.execute_reply":"2021-06-22T12:59:11.448129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's train the random forest.","metadata":{}},{"cell_type":"code","source":"n_splits = 3\n\n_, _, _ = perform_CV(RF_wrapper, n_splits=n_splits)","metadata":{"execution":{"iopub.status.busy":"2021-06-22T13:13:02.910385Z","iopub.execute_input":"2021-06-22T13:13:02.910868Z","iopub.status.idle":"2021-06-22T13:14:15.926446Z","shell.execute_reply.started":"2021-06-22T13:13:02.910829Z","shell.execute_reply":"2021-06-22T13:14:15.925183Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's repeat this for LGBM.","metadata":{}},{"cell_type":"code","source":"cat_features = ['building_id','meter','site_id','primary_use','hour','weekday']\n\nmodels, scores, oof_total = perform_CV(LGBM_wrapper, n_splits=n_splits)","metadata":{"execution":{"iopub.status.busy":"2021-06-22T13:04:53.958163Z","iopub.execute_input":"2021-06-22T13:04:53.958876Z","iopub.status.idle":"2021-06-22T13:09:36.611307Z","shell.execute_reply.started":"2021-06-22T13:04:53.95881Z","shell.execute_reply":"2021-06-22T13:09:36.609806Z"},"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. We can use this to retroactively drop further superfluous features during preprocessing.","metadata":{}},{"cell_type":"code","source":"importance = pd.DataFrame([model.feature_importance() for model in models],\n                          columns=X_train.columns,\n                          index=[f'Fold {i}' for i in range(1, n_splits + 1)])\nimportance = importance.T\nimportance['Average importance'] = importance.mean(axis=1)\nimportance = importance.sort_values(by='Average importance', ascending=False)\n\nplt.figure(figsize=(10,7))\nsns.barplot(x='Average importance', y=importance.index, data=importance);","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del X_train, y_train\ngc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test set, inference and submission","metadata":{}},{"cell_type":"markdown","source":"Let us load the test data and apply the same transformations as for the training set.","metadata":{}},{"cell_type":"code","source":"%%time\nX_test = merged_dfs('test')","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test.info()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's compute the predictions on the test set for each model, then average the results. This would improve the stability of the predictions, assuming the errors from each model are independent. Indeed, this gives a considerable boost in performance on PB.\n\nWe 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":"n_iterations = 20\nbatch_size = len(X_test) // n_iterations\n\npreds = []\nfor i in tqdm(range(n_iterations)):\n    start = i * batch_size\n    fold_preds = [np.expm1(model.predict(X_test.iloc[start:start + batch_size], \n                                         num_iteration=model.best_iteration)) for model in models]\n    preds.extend(np.mean(fold_preds, axis=0))\n\ndel X_test\ngc.collect()","metadata":{"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(preds, 0, a_max=None)})\nsubmission.to_csv('submission.csv', index=False)","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's it. Do upvote this kernel if you found it of any use! 🖖\n\nBelow are some links to kernels I either found interesting or directly borrowed from.","metadata":{}},{"cell_type":"markdown","source":"# References","metadata":{}},{"cell_type":"markdown","source":"Instructive kernels:\n- https://www.kaggle.com/purist1024/ashrae-simple-data-cleanup-lb-1-08-no-leaks\n- https://www.kaggle.com/gunesevitan/ashrae-lightgbm-1-048-no-leak\n\nInteresting discussions:\n- https://www.kaggle.com/kyakovlev/ashrae-cv-options/comments\n- https://www.kaggle.com/c/ashrae-energy-prediction/discussion/122471\n\nAnd summaries:\n- https://www.kaggle.com/c/ashrae-energy-prediction/discussion/125017\n- https://www.kaggle.com/c/ashrae-energy-prediction/discussion/112872#651685","metadata":{}}]}