{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# ASHRAE - Great Energy Predictor III\n", "\n", "\n", "**IMPORTANT NOTE!!** The inline graphics make the content frame over 16Mb of data - it may take a while to load and render fully, please have patience (I hope your browser renders this bit quickly at least :)\n", "\n", "\n", "## Heatmap EDA\n", "\n", "This notebook shows all the meter readings over time in a 2D heatmap (one per building), as a simple exploration of the time series structure of the train/test sets.\n", "\n", "Each plot shows a year of data: a pixel represents one hour, with a row of pixels as one week (24 * 7 = 168 hours), with ~52 rows, one per week of the year.\n", "\n", "The plots are saved to png files and can be downloaded from the **Output Files** section (left), for closer inspection.\n", "\n", "This gives a good overview of the kind of interactions between hourly, weekly and long term seasonality: many patterns are present, strong seasonal variation, weekly cycles, spikes, and intricate patterns that would be hard to spot on a simple line plot.\n", "\n", "The plot for every single building is shown at the end, sorted by site, and act as a map of the entire training set, illustrating the kinds of *conditionality* needed in the models, e.g. quiet Thursdays, busy summers, shorter usage-hours at weekends..."]}, {"cell_type": "code", "execution_count": 1, "metadata": {"_kg_hide-input": true}, "outputs": [], "source": ["%matplotlib inline\n", "import pandas as pd\n", "import numpy as np\n", "import gc, os, sys\n", "import matplotlib.pyplot as plt\n", "import cv2\n", "from IPython.display import Image, display, HTML\n", "import calendar\n", "\n", "np.seterr(divide='ignore', invalid='ignore')\n", "\n", "DTYPE = {\n", " 'building_id': ('int16'),\n", " 'meter': ('int8'),\n", " 'meter_reading': ('float32')\n", "}\n", "\n", "INPUT = '../input/ashrae-energy-prediction'"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": ["buildings = pd.read_csv(f'{INPUT}/building_metadata.csv', index_col='building_id')\n", "buildings.shape"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": ["train = pd.read_csv(f'{INPUT}/train.csv', parse_dates=['timestamp'], dtype=DTYPE)\n", "train.shape"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": ["weather = pd.read_csv(f'{INPUT}/weather_train.csv', parse_dates=['timestamp'])\n", "weather.shape"]}, {"cell_type": "markdown", "metadata": {}, "source": ["One-liners \"R\" us - first heatmap is count of building types at each site."]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": ["buildings.groupby(['primary_use','site_id']).size().unstack().fillna(0).style.background_gradient(axis=None)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create (x,y) coordinates for every row, a row is one week, in pandas, `dayofweek` 0 is *Monday*. (Note that all of train is in 2016.)"]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": ["train.timestamp.dt.year.value_counts()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["#### Note: pandas reports first days of training set as weekofyear=53 (of previous year, 2015); use mod 53 as an ugly fix to make those times 'week 0', then week of year is 0..52 inclusive"]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": ["train.timestamp.dt.weekofyear.min(), train.timestamp.dt.weekofyear.max()"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": ["def add_xy(df):\n", " dt = df.timestamp.dt\n", " df['x'] = ((dt.dayofweek * 24) + dt.hour).astype('int16')\n", " df['y'] = (dt.weekofyear % 53).astype('int8')"]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": ["add_xy(train)\n", "add_xy(weather)\n", "train.x.max(), train.y.max()"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": ["WIDTH = train.x.max() + 1\n", "HEIGHT = train.y.max() + 1\n", "WIDTH, HEIGHT"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Grayscale plots\n", "\n", "One per building / meter."]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": ["PLOTS_GRAY = 'plots_grayscale'\n", "os.makedirs(PLOTS_GRAY, exist_ok=True)"]}, {"cell_type": "code", "execution_count": 12, "metadata": {"_kg_hide-input": true}, "outputs": [], "source": ["def normalize(c):\n", " return np.nan_to_num(c / c.max())\n", "\n", "def log_normalize(c):\n", " return normalize(np.log1p(c))\n", "\n", "# WARNING! This stretches each color channel independently\n", "# - it loses relative scale between them\n", "def log_normalize_chan(c):\n", " return np.dstack([log_normalize(a.T) for a in c.T])\n", "\n", "def write_img(fname, img):\n", " img *= 255\n", " if len(img.shape) == 3:\n", " img = cv2.cvtColor(img.astype('uint8'), cv2.COLOR_RGB2BGR)\n", " cv2.imwrite(fname, img)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Simple counting code.\n", "\n", "Using [`np.add.at`][1] is not strictly necessary here.\n", "It ensures that results are accumulated for elements that are indexed more than once, so this code is more generically useful.\n", "However, all buildings have at most one observation per hour, so simple indexing `cc[df.y, df.x] += df.meter_reading` would also work.\n", "\n", "\n", "[1]: https://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.at.html"]}, {"cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": ["for (m, b), df in train.groupby(['meter', 'building_id']):\n", " cc = np.zeros((HEIGHT, WIDTH), dtype=float)\n", " np.add.at(cc, (df.y, df.x), df.meter_reading)\n", " write_img(f'{PLOTS_GRAY}/{m}_{b}_log.png', log_normalize(cc))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Sneak preview - main display is later..."]}, {"cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": ["display(Image(f'{PLOTS_GRAY}/1_161_log.png', width=WIDTH*4))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Over 2000 plots :)"]}, {"cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": ["!ls -1 $PLOTS_GRAY | wc -l"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# RGB Plots\n", "\n", "One per building - three meters in red, green and blue channels. Each channel is log transformed and stretched separately - so some information about absolute energy usage per channel is lost."]}, {"cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": ["PLOTS_RGB = 'plots_rgb'\n", "os.makedirs(PLOTS_RGB, exist_ok=True)"]}, {"cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": ["METERS = ['electricity', 'chilledwater', 'steam', 'hotwater']\n", "M_DICT = {k:i for i, k in enumerate(METERS)}\n", "\n", "train.meter.value_counts()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Values are: {0: electricity, 1: chilledwater, 2: steam, 3: hotwater}\n", "\n", "Simplest decision is: ignore hotwater (least frequent), plot the others;\n", "\n", " - red = **electricity**\n", " - green = **chilledwater**\n", " - blue = **steam**\n", "\n", "The brightness of each color channel indicates the energy usage, whilst black could indicate zero energy usage or (much more likely) missing data."]}, {"cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": ["rgb_pngs = {}\n", "for b, df in train.query('meter<3').groupby('building_id'):\n", " cc = np.zeros((HEIGHT, WIDTH, 3), dtype=float)\n", " np.add.at(cc, (df.y, df.x, df.meter), df.meter_reading)\n", " png = f'{PLOTS_RGB}/{b}_log.png'\n", " rgb_pngs[b] = png\n", " write_img(png, log_normalize_chan(cc))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Over 1400 plots (well, one for each building... :)"]}, {"cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": ["!ls -1 $PLOTS_RGB | wc -l"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Matplotlib Display\n", "\n", "Show some plots selected for interesting features...\n", "\n", "Many will contain 5, 6 or 7 clearly visible daily peaks - however the exact start and end of those peaks can change between buildings, though is quite correlated for buildings in the same site. Similarly - the length of peak usage per day can vary within a building, from looking dormant, to just being shorter than a 'usual' day."]}, {"cell_type": "code", "execution_count": 20, "metadata": {"_kg_hide-input": true}, "outputs": [], "source": ["# Find month offsets\n", "month_y_min = train.groupby(train.timestamp.dt.month).y.min()\n", "ylabels = np.asarray(calendar.month_abbr)[month_y_min.index]\n", "yticks = month_y_min.values\n", "\n", "# Tick to mark days - each 24 hours\n", "days = list(calendar.day_abbr)\n", "xlabels = days\n", "xticks = np.arange(7) * 24\n", "WEATHER_COLS = ['air_temperature', 'precip_depth_1_hr', 'sea_level_pressure']\n", "\n", "def count_gray(df):\n", " cc = np.zeros((HEIGHT, WIDTH), dtype=float)\n", " np.add.at(cc, (df.y, df.x), df.meter_reading)\n", " return cc\n", "\n", "def count_rgb(df):\n", " cc = np.zeros((HEIGHT, WIDTH, 3), dtype=float)\n", " np.add.at(cc, (df.y, df.x, df.meter), df.meter_reading)\n", " return cc\n", "\n", "# use rank transform of selected columns (fork to try others!)\n", "def weather_rgb(df):\n", " cc = np.zeros((HEIGHT, WIDTH, 3), dtype=float)\n", " for i, c in enumerate(WEATHER_COLS):\n", " cc[df.y, df.x, i] = df[c].rank(pct=True)\n", " return cc\n", "\n", "def detail_list(series):\n", " return ''.join([f'
  • {k}: {v}'\n", " for k,v in series.dropna().items()])\n", "\n", "def display_plot(plotdata, title):\n", " fig, ax = plt.subplots(figsize=(14, 6))\n", " c = ax.imshow(plotdata)\n", " ax.set_xticks(xticks, False)\n", " ax.set_xticklabels(xlabels)\n", " ax.set_yticks(yticks)\n", " ax.set_yticklabels(ylabels)\n", " ax.set_title(title)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "def show_plot(src_df, building_id, comment):\n", " df = src_df.query(f'(building_id=={building_id}) and (meter<3)')\n", " display(HTML(f'

    Building {building_id}

    '))\n", " display(HTML(detail_list(buildings.loc[building_id])))\n", " display(df.groupby('meter').meter_reading.agg(['count', 'mean', 'max']))\n", " display(HTML(f'
    {comment}'))\n", " p = log_normalize_chan(count_rgb(df))\n", " display_plot(p, f'Building {building_id}')\n", "\n", "COMMENTS = {\n", " 647: \"Some have a lot of missing data\",\n", " 675: \"Nice intricate pattern\",\n", " 677: \"(Thursday) Here's Why I Did Not Go to Work Today (I didn't know this song - Google suggested it.)\",\n", " 182: \"Single hour spikes much higher than the mean (see stats above) - what's he building in there? (See comments! ↓)\",\n", " 822: \"Looks more like an on/off pattern (constant usage for one day, no hourly variation), and reminds me of a ten minute switching pattern in chinese mobile advert click patterns :) \",\n", " 751: \"Remember, blue is steam: Some activity at weekends, but weekends have different seasonal pattern\",\n", " 1017: \"Abrupt drop in electricity, switches to chilledwater (green)?\",\n", " 1063: \"Monthly pattern in weekends?\",\n", " 747: \"Who likes Fruit Salads?\",\n", " 1247: \"Interesting interaction between hour of day and position in year - phase shifts over time.\",\n", " 1355: \"Looks like four or more separate regimes - longer days in the later half (that start one hour later), and missing data around March/April that is common to most site 15 buildings.\",\n", "}"]}, {"cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": ["# Add this back if you want separate weather plots\n", "# for site, df in weather.groupby('site_id'):\n", "# display_plot(weather_rgb(df), f'Weather at site {site}')"]}, {"cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": ["for building_id, comment in COMMENTS.items():\n", " show_plot(train, building_id, comment)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# BY SITE\n", "\n", "Enough commenting on individual patterns: just see for yourself - here they **ALL** are :D\n", "\n", "**Note**: they are sorted by **primary_use**, then PNG file size - so more complicated plots (harder to compress → higher file size) come first.\n", "\n", "This clearly shows sites tend to have missing observations (black areas) at the same time, e.g. site 15 as shown above, and in site 0, nearly all buildings have missing observations for January/February at the top of each plot.\n", "\n", "Weather for each site is displayed first. Notice how much noisier plots from the natural world look compared to the robotic world of power consumption (and also different missing data patterns).\n", "\n", "Colors for weather:\n", "\n", " - red = **air_temperature**\n", " - green = **precip_depth_1_hr**\n", " - blue = **sea_level_pressure**\n", " \n", " \n", "## Tips\n", "\n", "Apple Mac users can run \"Digital Color Meter\" to see the RGB components of pixels in these plots to examine exact values. (Windows & Linux must have equivalents — if anyone knows please add a comment.)\n", "\n", "MacOS also has an extremely useful full screen zoom feature that is as simple as **<control> + scroll gesture**.\n", "\n", "If this does not work it can be enabled in the Accessibility settings:\n", "\n", " System Preferences\n", " >> Accessibility\n", " >> Zoom\n", " >> \"Use scroll gesture with modifier keys to zoom\"\n"]}, {"cell_type": "code", "execution_count": 23, "metadata": {"_kg_hide-input": true}, "outputs": [], "source": ["import base64\n", "\n", "buildings['png'] = pd.Series(rgb_pngs)\n", "buildings['png_size'] = pd.Series(rgb_pngs).map(os.path.getsize)\n", "\n", "IMG_PER_ROW = 4\n", "\n", "def img_src(fname):\n", " with open(fname, 'rb') as f:\n", " return base64.b64encode(f.read()).decode('utf-8')\n", "\n", "# Return HTML for a table of plots\n", "def make_table(df, per_row, inline=True):\n", " src = \"\"\n", " for i, (bid, row) in enumerate(df.iterrows()):\n", " if (i%per_row) == 0:\n", " if i:\n", " src += \"\"\n", " src += \"\"\n", " if inline:\n", " dat = f\"data:image/png;base64,{img_src(row.png)}\"\n", " else:\n", " dat = row.png\n", " src += (f''\n", " f'
    [{bid}] {row.primary_use}')\n", " src += \"\"\n", " return f\"{src}
    \"\n", "\n", "# Write single HTML table with all plots\n", "with open('pngs_sorted_by_size.html', 'w') as f:\n", " table_src = make_table(buildings.sort_values('png_size'), per_row=8, inline=False)\n", " print(f\"ASHRAE Plots\"\n", " f\"{table_src}\\n\", file=f)\n", "\n", "# Generate output section per site\n", "for site, df in buildings.sort_values(['primary_use', 'png_size'], ascending=[True, False]).groupby('site_id'):\n", " display(HTML(f'

    Site {site}

    '))\n", " \n", " display(HTML(f'

    Weather

    '))\n", " display_plot(weather_rgb(weather.loc[weather.site_id==site]), f'Weather site {site}')\n", " display(weather.loc[weather.site_id == site].agg({\n", " WEATHER_COLS[0]: ['min', 'mean', 'max'],\n", " WEATHER_COLS[1]: ['min', 'mean', 'max'],\n", " WEATHER_COLS[2]: ['min', 'mean', 'max']\n", " }).T.round(1)\n", " )\n", " \n", " display(HTML(f'

    Stats

    '))\n", " display(\n", " df.groupby('primary_use').agg({\n", " 'square_feet': ['count', 'mean'],\n", " 'year_built': ['min', 'median', 'max'],\n", " 'floor_count': ['min', 'median', 'max']\n", " }).T.dropna(how='all').T.fillna('').style.background_gradient(axis=0)\n", " )\n", " \n", " display(HTML(f'

    Buildings

    '))\n", " src = make_table(df, IMG_PER_ROW, inline=True)\n", " src += \"
    \"\n", " display(HTML(src))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Test Set Structure Mini EDA\n", "\n", "I was going to plot the missing/present patterns in the test set but it turns out there is only one pattern: the test set is a perfect grid of all days in 2017-2018, all buildings, with selected meters per building..."]}, {"cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": ["test = pd.read_csv(f'{INPUT}/test.csv', usecols=['building_id', 'meter'], dtype=DTYPE)\n", "test.shape"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We expect hourly observations, so two years is:"]}, {"cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": ["2*365*24"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Buildings are present in multiples of 17520:"]}, {"cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": ["test.groupby(['building_id']).size().value_counts().sort_index()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The test set is simply 2380 sets of `['building_id', 'meter']` :"]}, {"cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": ["test.groupby(['building_id', 'meter']).size().value_counts()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["However... there *will* be missing data in the test set - so some submission rows will be *public*, some will be *private* and some will be *ignored*, scores are only computed for rows marked public or private.\n", "\n", "See this old [Home Depot Product Search Relevance thread][1] for an example of a [Kaggle solution file][2], containing the ground truth (`Relevance`) and the `Usage` column :)\n", "\n", "(Side note: the Home Depot competition, like many others, had a different type of *ignored* test set data - rows that are added just to make leaderboard probing harder.)\n", "\n", "***UPDATE:***\n", "\n", "It has now been announced [the set of rows that are *Ignored* will change before the deadline][3]. Sites with (test set era) 2017-2018 data publicly available will not be included in the private LB. \n", "\n", "\n", " [1]: https://www.kaggle.com/c/home-depot-product-search-relevance/discussion/20587\n", " [2]: https://storage.googleapis.com/kaggle-forum-message-attachments/117897/4154/solution.csv\n", " [3]: https://www.kaggle.com/c/ashrae-energy-prediction/discussion/117357"]}, {"cell_type": "markdown", "metadata": {}, "source": ["# Clean Up\n", "\n", "Compress the generated pngs - Kaggle does not allow more than 500 output files.\n", "\n", " -bd : disable progress indicator\n", " -mmt[N] : set number of CPU threads\n", " -sdel : delete files after compression"]}, {"cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": ["!7z a -bd -mmt4 -sdel {PLOTS_GRAY}.zip {PLOTS_GRAY} >>compress.log"]}, {"cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": ["!7z a -bd -mmt4 -sdel {PLOTS_RGB}.zip {PLOTS_RGB} >>compress.log"]}], "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.7.1"}}, "nbformat": 4, "nbformat_minor": 1}