{"cells":[{"metadata":{},"cell_type":"markdown","source":"# The worst meters\nIn this experiment, we evaluate each building/meter pair separately and chart the errors so that we can see which meters contribute most to our global error. The goal is to demonstrate where we should focus the most effort -- either to clean up bad meter-reading data or to produce models which can better deal with those particular readings.\n\nI believe that the harm done by the worst of the meters is two-fold: the impossibility of predicting them directly raises the error, but I suspect that it also degrades the overall model so that it performs less well even on the \"good\" meters. Further experimentation is required to figure out how to best mitigate this harm.\n\nNote that we already eliminate the first 141 days of electrical readings for site 0. We already know that those data points are bad, and that we best deal with them by simply dropping them. This experiment is intended to find out what *further* cleanup is required after this basic process."},{"metadata":{},"cell_type":"markdown","source":"## Summary\nAs expected, a small number of meters account for a moderate fraction of our overall error -- with just 4 meters accounting for 1%. Of particular note are building 1072, with the highest error contribution; and building 954, with 3 different meters in our top ten wall of shame. Interestingly, the much maligned buildings 1099 and 778 come in at #2 and #12, respectively, rather than headlining the chart. A surprising number of offenders were electrical meters (which we'd expect to be more predictable), starting with building 799 at #3.\n\nThe RMSLE and global contribution of the 10 worst meters are shown here:\n\n|Building|Meter|RMSLE|Contribution|\n|-------:|----:|----:|-----------:|\n|1072|2|6.685|0.141|\n|1099|2|5.801|0.122|\n|799|0|5.715|0.120|\n|954|1|5.650|0.119|\n|803|0|5.839|0.119|\n|1303|2|5.442|0.114|\n|1021|3|5.199|0.109|\n|954|2|4.967|0.104|\n|802|0|4.862|0.102|\n|954|0|4.783|0.101|\n"},{"metadata":{},"cell_type":"markdown","source":"# Framework\nThe framework code is taken from my previous kernel: [Strategy evaluation: What helps and by how much?](https://www.kaggle.com/purist1024/strategy-evaluation-what-helps-and-by-how-much). It is described in more detail there and so, in order to get to the point, we incorporate it here without the descriptions."},{"metadata":{"jupyter":{"source_hidden":true},"_kg_hide-input":true,"_kg_hide-output":true,"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)\n\n\ndef 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\n\ndef 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\n\ndef _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\n\ndef _add_time_features(X):\n    return X.assign(tm_day_of_week=((X.timestamp // 24) % 7), tm_hour_of_day=(X.timestamp % 24))\n\ndef 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\n\ncategorical_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":"# Code for per-building/meter evaluation\nOur `run_experiment` function evaluates model performance for each building & meter via a 6-way time-based cross-validation. It produces a chart giving the RMSLE for each model, and a \"contribution\" number which is the RMSLE that would be achieved over the whole dataset if we perfectly predicted every other building and meter. This makes it a crude approximation of the maximum improvement we could get if we fixed all of the problems with the given meter.\n\n`fractions_report` describes the number of meters required to account for specific fractions of the overall error.\n\n`plot_contributions` shows all of the contribution values in 4 easy charts. Note that matplotlib is perfectly willing to simply drop any lines which it considers \"too thin\", which is why we had to adjust the figure size to absurd values."},{"metadata":{"trusted":false},"cell_type":"code","source":"def run_experiment(n_estimators, sample_frac=1):\n    X, y = combined_train_data()\n\n    # Reduce evaluation cost by subsampling the data\n    X = X.sample(frac=sample_frac).sort_index()\n    y = y.reindex(X.index)\n\n    # Additional preprocessing\n    X, y = _drop_missing_site_0(X, y)\n    X = compress_dataframe(_add_time_features(X))\n    y = np.log1p(y)\n\n    model = LGBMWrapper(random_state=0, n_jobs=-1, n_estimators=n_estimators,\n                        categorical_feature=categorical_columns)\n    contribution_chart = pd.DataFrame(columns=[\"building_id\", \"meter\", \"RMSLE\", \"contribution\"])\n\n    for building in sorted(X.building_id.unique()):\n        for meter in range(4):\n            X_subset = X[(X.building_id == building) & (X.meter == meter)]\n            y_subset = y.reindex(X_subset.index)\n            if len(y_subset) == 0: continue\n            splits = make_8121_splits(X_subset, 1)  # We already subsampled, so no need to resample\n            X_subset = X_subset.drop(columns=\"timestamp\")\n\n            cv_y, cv_prediction = make_cv_predictions(model, splits, X_subset, y_subset,\n                                                      drop_electrical_zeros=False, verbose=False)\n            sle = np.square(cv_y - cv_prediction).sum()\n            rmsle = np.sqrt(sle / len(cv_y))\n            contribution = np.sqrt(sle / len(y))\n            contribution_chart.loc[len(contribution_chart)] = (building, meter, rmsle, contribution)\n    return contribution_chart","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def fractions_report(chart):\n    contribution = chart.contribution\n    cum_sum = (contribution.sort_values(ascending=False) / contribution.sum()).cumsum()\n    for frac in [0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 1]:\n        count = cum_sum[cum_sum < frac].count()+1  # Add one to account for boundary issues\n        frac_frac = count / len(cum_sum)\n        print(f\"{count: 5} meters ({frac_frac:6.1%} of total) account for {frac:4.0%} of error.\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"%matplotlib inline\nimport matplotlib.pyplot as plt\n\ndef plot_contributions(chart):\n    for meter in range(4):\n        subchart = chart[chart.meter == meter].sort_values(\"building_id\")\n        ax = plt.figure(0, (24, 8)).add_subplot(111)\n        ax.set(xlim=(0, chart.building_id.max() + 1), ylim=(0, chart.contribution.max() * 1.05))\n        ax.bar(subchart.building_id, subchart.contribution, width=1, label=f\"meter {meter}\")\n        ax.legend()\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluation\nWe evaluate the error contribution of each building/meter pair, and chart the results. The resulting \"contribution chart\" is also saved to a file so that you can easily view it in its entirety. We run a \"light\" (10 estimators) version of LGBM, both for evaluation speed, and because our experiments have shown that more estimators at the per-building level just cause overfitting.\n\nNote that the total contributions sum to a much higher total than the actual observed error. This is presumably because the global model *does* benefit from cross-building correlations. This doesn't negate the fact that some meter readings contain a much lower signal-to-noise ratio, and that we could benefit by cleaning them up or even eliminating them."},{"metadata":{"trusted":false},"cell_type":"code","source":"contribution_chart = run_experiment(10, 1)\ncontribution_chart.sort_values(\"contribution\", ascending=False).to_csv(\"contribution_chart.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Contributions by number of meters"},{"metadata":{"trusted":false},"cell_type":"code","source":"fractions_report(contribution_chart)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Contributions by meter type"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(contribution_chart.groupby(\"meter\").contribution.sum().to_frame().T)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Top contributions by building/meter"},{"metadata":{"trusted":false},"cell_type":"code","source":"display(contribution_chart.sort_values(\"contribution\", ascending=False).reset_index(drop=True).head(25))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## All contributions by meter"},{"metadata":{"trusted":false},"cell_type":"code","source":"plot_contributions(contribution_chart)","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}