{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Really simple cleanup\nIn this kernel, I provide basic routines to identify \"bad\" rows that can simply be dropped to improve overall performance.\n\nI identify four sorts of bad rows, with significant overlap:\n* Unjustified runs of zero readings: We identify sequences of more than 48 hours with zero readings which do not occur during the typical seasons for the designated meter type\n* Zero readings for electical meters: There's no reason for a building to ever have zero electrical usage, so we simply throw them all away.\n* The first 141 days of electricity for site 0: Most of these would be covered by the previous sets, but there are a few stray non-zero values that we ignore because they don't fit the overall pattern.\n* Abnormally high readings from building 1099: These values are just absurdly high and don't fit an established pattern. Leaderboard probes show that we do indeed benefit by dropping the outliers.\n\n"},{"metadata":{},"cell_type":"markdown","source":"## Results\nAfter deleting these bad rows, adding basic features, and fitting/predicting with a proven-useful model, I submit the predictions and get a LB score of 1.08. Note that I have made use of no data leaks, so this is a respectable (but not record-breaking) score."},{"metadata":{},"cell_type":"markdown","source":"# Data cleanup routines\nPresented here are fast, fully-debugged routines for identifying bad rows. You should be able to incorporate them into your code whether or not you use the rest of the code in this kerel. You simply need to make sure that your 'timestamp' column has been converted into hours since Jan 1, 2016.\n\nAlternatively, you can simply incorporate the list of indices that we write out later in the kernel."},{"metadata":{},"cell_type":"markdown","source":"`find_bad_zeros` identifies rows with zero-readings that have the following characteristics:\n\n* all electrical readings with zero values.\n* 48+ hour runs of steam and hotwater zero-readings *except* for those which are entirely contained within what we identify as typical \"core summer months\".\n* 48+ hour runs of chilledwater zero-readings *except* for those which occur simultaneously at the start and end of the year (i.e. \"core winter months\").\n\nThe exact time periods for summer and winter were determined empirically by looking at hundreds of charts for meters which performed particularly poorly. (See [\"The worst meters\"](https://www.kaggle.com/purist1024/ashrae-the-worst-meters) for the code which identifies these poor performers.)"},{"metadata":{"trusted":false},"cell_type":"code","source":"def make_is_bad_zero(Xy_subset, min_interval=48, summer_start=3000, summer_end=7500):\n    \"\"\"Helper routine for 'find_bad_zeros'.\n    \n    This operates upon a single dataframe produced by 'groupby'. We expect an \n    additional column 'meter_id' which is a duplicate of 'meter' because groupby \n    eliminates the original one.\"\"\"\n    meter = Xy_subset.meter_id.iloc[0]\n    is_zero = Xy_subset.meter_reading == 0\n    if meter == 0:\n        # Electrical meters should never be zero. Keep all zero-readings in this table so that\n        # they will all be dropped in the train set.\n        return is_zero\n\n    transitions = (is_zero != is_zero.shift(1))\n    all_sequence_ids = transitions.cumsum()\n    ids = all_sequence_ids[is_zero].rename(\"ids\")\n    if meter in [2, 3]:\n        # It's normal for steam and hotwater to be turned off during the summer\n        keep = set(ids[(Xy_subset.timestamp < summer_start) |\n                       (Xy_subset.timestamp > summer_end)].unique())\n        is_bad = ids.isin(keep) & (ids.map(ids.value_counts()) >= min_interval)\n    elif meter == 1:\n        time_ids = ids.to_frame().join(Xy_subset.timestamp).set_index(\"timestamp\").ids\n        is_bad = ids.map(ids.value_counts()) >= min_interval\n\n        # Cold water may be turned off during the winter\n        jan_id = time_ids.get(0, False)\n        dec_id = time_ids.get(8283, False)\n        if (jan_id and dec_id and jan_id == time_ids.get(500, False) and\n                dec_id == time_ids.get(8783, False)):\n            is_bad = is_bad & (~(ids.isin(set([jan_id, dec_id]))))\n    else:\n        raise Exception(f\"Unexpected meter type: {meter}\")\n\n    result = is_zero.copy()\n    result.update(is_bad)\n    return result\n\ndef find_bad_zeros(X, y):\n    \"\"\"Returns an Index object containing only the rows which should be deleted.\"\"\"\n    Xy = X.assign(meter_reading=y, meter_id=X.meter)\n    is_bad_zero = Xy.groupby([\"building_id\", \"meter\"]).apply(make_is_bad_zero)\n    return is_bad_zero[is_bad_zero].index.droplevel([0, 1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`find_bad_sitezero` identifies the \"known-bad\" electrical readings from the first 141 days of the data for site 0 (i.e. UCF)."},{"metadata":{"trusted":false},"cell_type":"code","source":"def find_bad_sitezero(X):\n    \"\"\"Returns indices of bad rows from the early days of Site 0 (UCF).\"\"\"\n    return X[(X.timestamp < 3378) & (X.site_id == 0) & (X.meter == 0)].index","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`find_bad_building1099` identifies the most absurdly high readings from building 1099. These are orders of magnitude higher than all data, and have been emperically seen in LB probes to be harmful outliers."},{"metadata":{"trusted":false},"cell_type":"code","source":"def find_bad_building1099(X, y):\n    \"\"\"Returns indices of bad rows (with absurdly high readings) from building 1099.\"\"\"\n    return X[(X.building_id == 1099) & (X.meter == 2) & (y > 3e4)].index","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, `find_bad_rows` combines all of the above together to allow you to do a one-line cleanup of your data."},{"metadata":{"trusted":false},"cell_type":"code","source":"def find_bad_rows(X, y):\n    return find_bad_zeros(X, y).union(find_bad_sitezero(X)).union(find_bad_building1099(X, y))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Framework\nThe following 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":{"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\n\npd.set_option(\"max_columns\", 500)\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 _add_time_features(X):\n    return X.assign(tm_day_of_week=((X.timestamp // 24) % 7), tm_hour_of_day=(X.timestamp % 24))\n\nclass 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\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":"The following functions are simple variants of the ones above, but deal with loading in the test set. They could easily be refactored to share code with those functions, but we keep them separate for this demonstration."},{"metadata":{"trusted":false},"cell_type":"code","source":"def read_test():\n    df = pd.read_csv(input_file(\"test.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    return compress_dataframe(df).set_index(\"row_id\")\n\ndef read_weather_test(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    df = pd.read_csv(input_file(\"weather_test.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, 26304))\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_test_data(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    X = compress_dataframe(read_test().join(read_building_metadata(), on=\"building_id\").join(\n        read_weather_test(fix_timestamps, interpolate_na, add_na_indicators),\n        on=[\"site_id\", \"timestamp\"]).fillna(-1))\n    return X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fit, predict, and submit\n"},{"metadata":{},"cell_type":"markdown","source":"First, read in the data, and identify the bad rows using the provided functions. Write out the bad rows in an index-per-line format for fast and easy re-use."},{"metadata":{"trusted":false},"cell_type":"code","source":"X, y = combined_train_data()\n\nbad_rows = find_bad_rows(X, y)\npd.Series(bad_rows.sort_values()).to_csv(\"rows_to_drop.csv\", header=False, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Drop the bad rows that we identified above, and then train the model using our favorite features and regressors. See [Strategy evaluation: What helps and by how much?](https://www.kaggle.com/purist1024/strategy-evaluation-what-helps-and-by-how-much) for information on the specific strategies."},{"metadata":{"trusted":false},"cell_type":"code","source":"X = X.drop(index=bad_rows)\ny = y.reindex_like(X)\n\n# Additional preprocessing\nX = compress_dataframe(_add_time_features(X))\nX = X.drop(columns=\"timestamp\")  # Raw timestamp doesn't help when prediction\ny = np.log1p(y)\n\nmodel = CatSplitRegressor(\n    LGBMWrapper(random_state=0, n_jobs=-1, categorical_feature=categorical_columns), \"meter\")\n\nmodel.fit(X, y)\ndel X, y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load the test set and predict meter readings. We must, of course, use exponentiation to convert our predictions back from log-scale to the desired kWh values. We also clip to a minimum of zero, since we know that there will be no negative readings."},{"metadata":{"trusted":false},"cell_type":"code","source":"X = combined_test_data()\nX = compress_dataframe(_add_time_features(X))\nX = X.drop(columns=\"timestamp\")  # Raw timestamp doesn't help when prediction\n\npredictions = pd.DataFrame({\n    \"row_id\": X.index,\n    \"meter_reading\": np.clip(np.expm1(model.predict(X)), 0, None)\n})\n\ndel X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, write the predictions out for submission. After that, it's Miller Time (tm)."},{"metadata":{"trusted":false},"cell_type":"code","source":"predictions.to_csv(\"submission.csv\", index=False, float_format=\"%.4f\")","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}