{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Strategy Evaluation\nThis competition provides an interesting grab-bag of sub-problems to solve: the data is messy and needs cleaned (site 0); the data sources are inconsistent (meter timestamps vs. weather timestamps); it asks you to predict several different vaguely related sets of values (meter type); it has hugely important categorical features (building_id); and we can find lots of external data to supplement the provided data. This makes it hard to know where to start and what areas are most important to tackle.\n\nLike everyone, I've tried a bunch of strategies and techniques designed to tackle all of the above issues as well the usual machine-learning pitfalls. However, as I built up layers of approaches, it became hard to keep track of which ones were actually useful and which ones simply didn't hurt. Since I'm an experimentalist at heart, I decided to build a framework which let me evaluate the actual effect of a bunch of different strategies. Since this is Kaggle, I decided to share the insights with the rest of the community."},{"metadata":{},"cell_type":"markdown","source":"## Summary\nAll of the code and the raw results are included below, and are well worth reading for additional insights and coding tricks. However, they can be summarized here. The percentage \"effect\" for each strategy is how much *worse* our results are when we don't use it.\n\n|Strategy|Effect|Details|\n|:-------|-----:|:------|\n|**Log-scaled target**|**45.7%**|Use log1p to rescale the meter readings, and expm1 to rescale the predictions back to kWh.|\n|**Site 0 data cleanup**|**17.8%**|Simply drop site 0 electrical data from the first 141 days.\n|**LGBM category-awareness**|**9.7%**|Declare which columns are categorical data.\n|**Per-meter models**|**3.9%**|Use separate models to predict for each meter type and then re-combine the predictions.\n|**Time-based features**|**3.6%**|Perform basic feature engineering, with hour-of-day and day-of-week.\n|**Timestamp alignment**|**0.4%**|Convert the GMT-based weather data to align with local-time-based readings.\n|**Interpolated weather data**|**0.4%**|Impute missing weather values by interpolating from time-adjacent readings.\n|**Missing value indicators**|**0.1%**|Add flags to inform the model of which values were imputed.\n|**0-valued meter cleanup**|**0.0%**|Drop all 0-valued electrical readings, assuming that there's always some electrical use, and that zeros indicate anomalies."},{"metadata":{},"cell_type":"markdown","source":"# Framework"},{"metadata":{"trusted":false},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os\nimport warnings\n\nfrom lightgbm import LGBMRegressor\nfrom sklearn.base import BaseEstimator, RegressorMixin, clone\nfrom sklearn.metrics import mean_squared_log_error\nfrom IPython.display import HTML\n\npd.set_option(\"max_columns\", 500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utilities\n`input_file` is a simple wrapper that hides the input directory location and automatically handles compressed files.\n\n`compress_dataframe` uses Pandas' built-in capabilities to downcast data into its smallest practical form."},{"metadata":{"trusted":false},"cell_type":"code","source":"def input_file(file):\n    path = f\"../input/ashrae-energy-prediction/{file}\"\n    if not os.path.exists(path): return path + \".gz\"\n    return path\n\ndef compress_dataframe(df):\n    result = df.copy()\n    for col in result.columns:\n        col_data = result[col]\n        dn = col_data.dtype.name\n        if dn == \"object\":\n            result[col] = pd.to_numeric(col_data.astype(\"category\").cat.codes, downcast=\"integer\")\n        elif dn == \"bool\":\n            result[col] = col_data.astype(\"int8\")\n        elif dn.startswith(\"int\") or (col_data.round() == col_data).all():\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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data loading and preprocessing\nCode to read and combine the standard input files, converting timestamps to number of hours since the beginning of 2016.\n\nOptionally, we perform cleanups on the weather data: adjusting timestamps to the appropriate time zone; using interpolation to fill missing values; and adding NaN presence indicators."},{"metadata":{"trusted":false},"cell_type":"code","source":"def read_train():\n    df = pd.read_csv(input_file(\"train.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    return compress_dataframe(df)\n\ndef read_building_metadata():\n    return compress_dataframe(pd.read_csv(\n        input_file(\"building_metadata.csv\")).fillna(-1)).set_index(\"building_id\")\n\nsite_GMT_offsets = [-5, 0, -7, -5, -8, 0, -5, -5, -5, -6, -7, -5, 0, -6, -5, -5]\n\ndef read_weather_train(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    df = pd.read_csv(input_file(\"weather_train.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    if fix_timestamps:\n        GMT_offset_map = {site: offset for site, offset in enumerate(site_GMT_offsets)}\n        df.timestamp = df.timestamp + df.site_id.map(GMT_offset_map)\n    if interpolate_na:\n        site_dfs = []\n        for site_id in df.site_id.unique():\n            # Make sure that we include all possible hours so that we can interpolate evenly\n            site_df = df[df.site_id == site_id].set_index(\"timestamp\").reindex(range(8784))\n            site_df.site_id = site_id\n            for col in [c for c in site_df.columns if c != \"site_id\"]:\n                if add_na_indicators: site_df[f\"had_{col}\"] = ~site_df[col].isna()\n                site_df[col] = site_df[col].interpolate(limit_direction='both', method='linear')\n                # Some sites are completely missing some columns, so use this fallback\n                site_df[col] = site_df[col].fillna(df[col].median())\n            site_dfs.append(site_df)\n        df = pd.concat(site_dfs).reset_index()  # make timestamp back into a regular column\n    elif add_na_indicators:\n        for col in df.columns:\n            if df[col].isna().any(): df[f\"had_{col}\"] = ~df[col].isna()\n    return compress_dataframe(df).set_index([\"site_id\", \"timestamp\"])\n\ndef combined_train_data(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    Xy = compress_dataframe(read_train().join(read_building_metadata(), on=\"building_id\").join(\n        read_weather_train(fix_timestamps, interpolate_na, add_na_indicators),\n        on=[\"site_id\", \"timestamp\"]).fillna(-1))\n    return Xy.drop(columns=[\"meter_reading\"]), Xy.meter_reading","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Drop potentially bad data\nThere are two different optional strategies here that can be used together. `_drop_electrical_zeros` removes every electrical meter reading with a value of 0; while `_drop_missing_site_0` gets rid of all electrical meter readings for site zero during the 141 days that its readings were mostly pegged to 0."},{"metadata":{"trusted":false},"cell_type":"code","source":"def _drop_electrical_zeros(X, y):\n    X = X[(y > 0) | (X.meter != 0)]\n    y = y.reindex(X.index)\n    return X, y\n\ndef _drop_missing_site_0(X, y):\n    X = X[(X.timestamp >= 3378) | (X.site_id != 0) | (X.meter != 0)]\n    y = y.reindex(X.index)\n    return X, y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Add basic time features\nWe add columns for \"hour of day\" and \"day of week\". Note that the source \"timestamp\" columns will be removed before training, so these are the only time-oriented features that will remain."},{"metadata":{"trusted":false},"cell_type":"code","source":"def _add_time_features(X):\n    return X.assign(tm_day_of_week=((X.timestamp // 24) % 7), tm_hour_of_day=(X.timestamp % 24))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cross-validation\nWe use a time-oriented split which deals in 12 equal time slices (which are almost, but not quite, months). Each fold includes a two adjacent slices for validation, 8 adjacent slices for training, and a 1 slice buffer on each side of the validation set. We use modulo arithmetic to wrap-around the year-end boundary. We thus have slices 4-11 vs 1-2; 6-12&1 vs 3-4; 8-12&1-3 vs 5-6; etc.\n\nWhile the CV score is always better than the LB score, relative differences between strategies tend to align well between CV and LB.\n\nBecause the data set is large, we use the `sample frac` parameter to train and evaluate on a specified fraction of the dataset. Again, this seems to track the full set fairly well with values of 0.05 and higher."},{"metadata":{"trusted":false},"cell_type":"code","source":"def np_sample(a, frac):\n    return a if frac == 1 else np.random.choice(a, int(len(a) * frac), replace=False)\n\ndef make_8121_splits(X, sample_frac):\n    np.random.seed(0)\n    time_sorted_idx = np.argsort(X.timestamp.values, kind='stable')\n    sections = np.array_split(time_sorted_idx, 12)\n    folds = []\n    for start_ix in range(0, 12, 2):\n        val_idxs = np.concatenate(sections[start_ix:start_ix + 2])  # no modulo necessary\n        train_idxs = np.concatenate(\n            [sections[ix % 12] for ix in range(start_ix + 3, start_ix + 11)])\n        folds.append((np_sample(train_idxs, sample_frac), np_sample(val_idxs, sample_frac)))\n    return folds\n\ndef make_cv_predictions(model, split, X, y, drop_electrical_zeros, verbose=True):\n    preds = []\n    for ix, (train_fold, val_fold) in enumerate(split):\n        Xt = X.iloc[train_fold]\n        yt = y.reindex_like(Xt)\n        if drop_electrical_zeros:\n            Xt, yt = _drop_electrical_zeros(Xt, yt)\n        Xv = X.iloc[val_fold]\n        yv = y.reindex_like(Xv)\n        if verbose: print(f\"Testing split {ix}: {len(Xt)} train rows & {len(Xv)} val rows\")\n        model.fit(Xt, yt)\n        preds.append(pd.DataFrame(dict(target=yv, prediction=model.predict(Xv)), index=yv.index))\n    result = pd.concat(preds).sort_index()\n    return result.target, result.prediction","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Per-meter models\nThe `CatSplitRegressor` wrapper class hides the process of splitting the training and validation data up by unique values of a specified column (i.e. \"meter\"), using a different sub-model for each value, and then re-integrating the resulting predictions. "},{"metadata":{"trusted":false},"cell_type":"code","source":"class CatSplitRegressor(BaseEstimator, RegressorMixin):\n    def __init__(self, model, col):\n        self.model = model\n        self.col = col\n\n    def fit(self, X, y):\n        self.fitted = {}\n        importances = []\n        for val in X[self.col].unique():\n            X1 = X[X[self.col] == val].drop(columns=[self.col])\n            self.fitted[val] = clone(self.model).fit(X1, y.reindex_like(X1))\n            importances.append(self.fitted[val].feature_importances_)\n            del X1\n        fi = np.average(importances, axis=0)\n        col_index = list(X.columns).index(self.col)\n        self.feature_importances_ = [*fi[:col_index], 0, *fi[col_index:]]\n        return self\n\n    def predict(self, X):\n        result = np.zeros(len(X))\n        for val in X[self.col].unique():\n            ix = np.nonzero((X[self.col] == val).to_numpy())\n            predictions = self.fitted[val].predict(X.iloc[ix].drop(columns=[self.col]))\n            result[ix] = predictions\n        return result","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## LGBM category awareness\nThe LGBM uses some very effective strategies to produce better splits when it knows that a feature is categorical. Since we have a reasonable number of integer-valued features that are, in fact, unordered or cyclical categoricals, we benefit from explicitly declaring them to LGBM.\n\nUnfortunately, the scikit wrapper from LGBM doesn't handle these declarations gracefully, and prints a warning that we've overridden its default guesses. We include extra code here to silence those warnings."},{"metadata":{"trusted":false},"cell_type":"code","source":"categorical_columns = [\n    \"building_id\", \"meter\", \"site_id\", \"primary_use\", \"had_air_temperature\", \"had_cloud_coverage\",\n    \"had_dew_temperature\", \"had_precip_depth_1_hr\", \"had_sea_level_pressure\", \"had_wind_direction\",\n    \"had_wind_speed\", \"tm_day_of_week\", \"tm_hour_of_day\"\n]\n\nclass LGBMWrapper(BaseEstimator, RegressorMixin):\n    def __init__(self, categorical_feature=None, **params):\n        self.model = LGBMRegressor(**params)\n        self.categorical_feature = categorical_feature\n\n    def fit(self, X, y):\n        with warnings.catch_warnings():\n            cats = None if self.categorical_feature is None else list(\n                X.columns.intersection(self.categorical_feature))\n            warnings.filterwarnings(\"ignore\",\n                                    \"categorical_feature in Dataset is overridden\".lower())\n            self.model.fit(X, y, **({} if cats is None else {\"categorical_feature\": cats}))\n            self.feature_importances_ = self.model.feature_importances_\n            return self\n\n    def predict(self, X):\n        return self.model.predict(X)\n\n    def get_params(self, deep=True):\n        return {**self.model.get_params(deep), \"categorical_feature\": self.categorical_feature}\n\n    def set_params(self, **params):\n        ctf = params.pop(\"categorical_feature\", None)\n        if ctf is not None: self.categorical_feature = ctf\n        self.model.set_params(params)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Catch-all experimental framework\nThis loads the data, runs cross-validation, reports the CV error, and displays the most important features as declared by the model.\n\nThe extra complexity comes from allowing us to switch various strategies on and off:\n* fix_timestamps: adjust for time-zone errors in weather data\n* interpolate_na: use simple linear interpolation to provide fully-populated weather data\n* add_na_indicators: add boolean columns which let the model know that values were imputed\n* drop_electrical_zeros: drops all electrical meter readings which are zero, *but only from the training set*. We still evaluate our performance on such readings.\n* drop_missing_site_0: drops all electical meter readings for the first 141 days for site 0. These are known to be bogus, so we neither train *nor* evaluate on them\n* log_target: convert the target variable to a logarithmic scale, and convert the predictions back.\n* LGBM_cataware: declare to LGBM the columns which are known to be categorical\n* add_time_features: add basic day of week and hour of day features.\n* per_meter_models: train separate models for each distinct meter type\n\nBy switching *off* one of these strategies at a time, we can get an indication of how much worse we are without them, and thus how much they helped out the baseline \"all strategies on\" evaluation."},{"metadata":{"trusted":false},"cell_type":"code","source":"def run_experiment(fix_timestamps=True, interpolate_na=True, add_na_indicators=True,\n                   drop_electrical_zeros=True, drop_missing_site_0=True, log_target=True,\n                   LGBM_cataware=True, add_time_features=True, per_meter_models=True,\n                   sample_frac=0.2, verbose=False, baseline='return'):\n    if verbose: print(f\"Loading data\")\n    X, y = combined_train_data(fix_timestamps, interpolate_na, add_na_indicators)\n    if verbose: X.info(memory_usage=\"deep\", verbose=False)\n\n    if add_time_features:\n        if verbose: print(f\"Adding time features\")\n        X = compress_dataframe(_add_time_features(X))\n        if verbose: X.info(memory_usage=\"deep\", verbose=False)\n\n    if drop_missing_site_0:\n        if verbose: print(f\"Dropping missing site 0 data\")\n        # We do this here rather than without our CV, because we are *certain* that there is no\n        # equivalent data in the test set.\n        X, y = _drop_missing_site_0(X, y)\n        if verbose: X.info(memory_usage=\"deep\", verbose=False)\n\n    model = LGBMWrapper(random_state=0, n_jobs=-1,\n                        categorical_feature=None if not LGBM_cataware else categorical_columns)\n    if per_meter_models:\n        model = CatSplitRegressor(model, \"meter\")\n\n    if verbose: print(f\"Creating splits\")\n    splits = make_8121_splits(X, sample_frac)\n\n    X = X.drop(columns=\"timestamp\")  # Raw timestamp doesn't carry over to test data\n\n    if log_target:\n        y = np.log1p(y)\n    if verbose: print(\"Making predictions\")\n    sampled_y, sampled_prediction = make_cv_predictions(model, splits, X, y, drop_electrical_zeros,\n                                                        verbose=verbose)\n    if log_target:\n        sampled_y = np.expm1(sampled_y)\n        sampled_prediction = np.expm1(sampled_prediction)\n    rmsle = np.sqrt(mean_squared_log_error(sampled_y, np.clip(sampled_prediction, 0, None)))\n    if baseline == 'return':\n        print(f\"RMSLE = {rmsle:.4f} (baseline)\")\n    else:\n        baseline_change_pct = ((rmsle/baseline) - 1) * 100\n        print(f\"RMSLE = {rmsle:.4f} ({baseline_change_pct:.1f}% over baseline)\")\n\n    if verbose: print(\"Finding important features\")\n    X_feature_importance = X.sample(frac=sample_frac, random_state=0)\n    y_feature_importance = y.reindex(X_feature_importance.index)\n    model.fit(X_feature_importance, y_feature_importance)\n    importances = pd.Series(model.feature_importances_, index=X.columns).rename(\"Importance\")\n    display(importances.sort_values(ascending=False).to_frame().T)\n\n    if baseline == 'return':\n        return rmsle","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The experiments"},{"metadata":{},"cell_type":"markdown","source":"## All strategies enabled\nThis is the baseline against which we evaluate each individual strategy. All subsequent experiments will *disable* a single strategy and report how much *worse* we perform compared to this evaluation."},{"metadata":{"trusted":false},"cell_type":"code","source":"baseline_score = run_experiment()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No log-scaled target\nGiven that we are evaluated on log-error, it's not surprising that we do better when operating in log-space. The naive linear-space approach is our biggest loser at 45% worse."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, log_target=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Keep bad site 0 readings\nThose bad readings for a single site really do drag everything else down -- by a full 17%. Without pre-emptively throwing them away, we'll never achieve greatness."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, drop_missing_site_0=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No LGBM category awareness\nThe category-aware capabilities of LGBM are unreasonably effective, netting us a 9% improvement. This is a lot of extra power for almost no cost."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, LGBM_cataware=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No per-meter models\nThe different types of meters are pretty much completely unrelated to each other, so it makes sense to split them up and better take advantage of true commonalities between similar meters. At 4% improvement, it's not mind-bogglingly impressive, but it's still a no-brainer."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, per_meter_models=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No time features\nIt's hard to imagine *not* using these basic time features -- hour of day and day of week. We all know that power usage will be lower during non-peak hours. However, the actual effect is fairly small. Oh well, c'est la vie."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, add_time_features=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No adjusted weather timestamps\nThis is a surprising result. We know that the timestamps for the weather data aren't aligned with the timestamps for the meter readings, and can be off by up to 8 hours. However, it seems that we lose less than half a percent in performance if we simply ignore this appalling condition."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, fix_timestamps=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No interpolated weather\nThere's a lot of missing weather data, but it seems that it really doesn't much matter whether we fix it. This seems odd."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, interpolate_na=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## No NA indicators\nGiven the interpolation that we did, there isn't likely to be much difference between an imputed value and the true value. Thus, flagging the imputed values doesn't really give the model much useful information. We can probably just drop the indicators to save a bit of memory."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, add_na_indicators=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Keep electrical meter zero readings\nThis is a misleading result. In fact, the point of dropping all readings of 0 for electrical meters is to provide an alternate method of dealing with the bad zeros in the first part of the site 0 data. However, unlike with the \"drop bad site 0 readings\" strategy, we dropped these values from training but not evaluation. The fact that we do better by dropping them than by keeping them indicates that zero-readings for electricity really are misleading outliers and that dropping all of them is a valid alternative to specially handling the ones from site 0."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, drop_electrical_zeros=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Disable all strategies\nThis is another baseline, showing what would happen if we just used the raw data and didn't try to improve anything. Not surprisingly, it does *really* badly."},{"metadata":{"trusted":false},"cell_type":"code","source":"run_experiment(baseline=baseline_score, fix_timestamps=False, interpolate_na=False,\n               add_na_indicators=False, drop_electrical_zeros=False, drop_missing_site_0=False,\n               log_target=False, LGBM_cataware=False, add_time_features=False,\n               per_meter_models=False)","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.7"}},"nbformat":4,"nbformat_minor":1}