modern_4_performance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance\n",
    "\n",
    "[Wes McKinney](https://twitter.com/wesmckinn), the creator of pandas, is kind of obsessed with performance. From micro-optimizations for element access, to [embedding](https://github.com/pydata/pandas/tree/master/pandas/src/klib) a fast hash table inside pandas, we all benefit from his and others' hard work.\n",
    "This post will focus mainly on making efficient use of pandas and NumPy.\n",
    "\n",
    "One thing I'll explicitly not touch on is storage formats.\n",
    "Performance is just one of many factors that go into choosing a storage format.\n",
    "Just know that pandas can talk to [many formats](http://pandas.pydata.org/pandas-docs/version/0.18.0/io.html), and the format that strikes the right balance between performance, portability, data-types, metadata handling, etc., is an [ongoing](http://blog.cloudera.com/blog/2016/03/feather-a-fast-on-disk-format-for-data-frames-for-r-and-python-powered-by-apache-arrow/) topic of discussion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "if int(os.environ.get(\"MODERN_PANDAS_EPUB\", 0)):\n",
    "    import prep # noqa\n",
    "\n",
    "sns.set_style('ticks')\n",
    "sns.set_context('talk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructors\n",
    "\n",
    "It's pretty common to have many similar sources (say a bunch of CSVs) that need to be combined into a single DataFrame. There are two routes to the same end:\n",
    "\n",
    "1. Initialize one DataFrame and append to that\n",
    "2. Make many smaller DataFrames and concatenate at the end\n",
    "\n",
    "For pandas, the second option is faster.\n",
    "DataFrame appends are expensive relative to a list append.\n",
    "Depending on the values, pandas might have to recast the data to a different type.\n",
    "And indexes are immutable, so each time you append pandas has to create an entirely new one.\n",
    "\n",
    "In the last section we downloaded a bunch of weather files, one per state, writing each to a separate CSV.\n",
    "One could imagine coming back later to read them in, using the following code.\n",
    "\n",
    "The idiomatic python way\n",
    "\n",
    "```python\n",
    "files = glob.glob('weather/*.csv')\n",
    "columns = ['station', 'date', 'tmpf', 'relh', 'sped', 'mslp',\n",
    "           'p01i', 'vsby', 'gust_mph', 'skyc1', 'skyc2', 'skyc3']\n",
    "\n",
    "# init empty DataFrame, like you might for a list\n",
    "weather = pd.DataFrame(columns=columns)\n",
    "\n",
    "for fp in files:\n",
    "    city = pd.read_csv(fp, columns=columns)\n",
    "    weather.append(city)\n",
    "```\n",
    "\n",
    "This is pretty standard code, quite similar to building up a list of tuples, say.\n",
    "The only nitpick is that you'd probably use a list-comprehension if you were just making a list.\n",
    "But we don't have special syntax for DataFrame-comprehensions (if only), so you'd fall back to the \"initialize empty container, append to said container\" pattern.\n",
    "\n",
    "But there's a better, pandorable, way\n",
    "\n",
    "```python\n",
    "files = glob.glob('weather/*.csv')\n",
    "weather_dfs = [pd.read_csv(fp, names=columns) for fp in files]\n",
    "weather = pd.concat(weather_dfs)\n",
    "```\n",
    "\n",
    "Subjectively this is cleaner and more beautiful.\n",
    "There's fewer lines of code.\n",
    "You don't have this extraneous detail of building an empty DataFrame.\n",
    "And objectively the pandorable way is faster, as we'll test next."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll define two functions for building an identical DataFrame. The first `append_df`, creates an empty DataFrame and appends to it. The second, `concat_df`,  creates many DataFrames, and concatenates them at the end. We also write a short decorator that runs the functions a handful of times and records the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "size_per = 5000\n",
    "N = 100\n",
    "cols = list('abcd')\n",
    "\n",
    "def timed(n=30):\n",
    "    '''\n",
    "    Running a microbenchmark. Never use this.\n",
    "    '''\n",
    "    def deco(func):\n",
    "        def wrapper(*args, **kwargs):\n",
    "            timings = []\n",
    "            for i in range(n):\n",
    "                t0 = time.time()\n",
    "                func(*args, **kwargs)\n",
    "                t1 = time.time()\n",
    "                timings.append(t1 - t0)\n",
    "            return timings\n",
    "        return wrapper\n",
    "    return deco\n",
    "    \n",
    "@timed(60)\n",
    "def append_df():\n",
    "    '''\n",
    "    The pythonic (bad) way\n",
    "    '''\n",
    "    df = pd.DataFrame(columns=cols)\n",
    "    for _ in range(N):\n",
    "        df.append(pd.DataFrame(np.random.randn(size_per, 4), columns=cols))\n",
    "    return df\n",
    "\n",
    "@timed(60)\n",
    "def concat_df():\n",
    "    '''\n",
    "    The pandorabe (good) way\n",
    "    '''\n",
    "    dfs = [pd.DataFrame(np.random.randn(size_per, 4), columns=cols)\n",
    "           for _ in range(N)]\n",
    "    return pd.concat(dfs, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>level_0</th>\n",
       "      <th>Method</th>\n",
       "      <th>Time (s)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>Append</td>\n",
       "      <td>0.202072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>Concat</td>\n",
       "      <td>0.116064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>Append</td>\n",
       "      <td>0.153113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>Concat</td>\n",
       "      <td>0.107477</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>Append</td>\n",
       "      <td>0.161898</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   level_0  Method  Time (s)\n",
       "0        0  Append  0.202072\n",
       "1        0  Concat  0.116064\n",
       "2        1  Append  0.153113\n",
       "3        1  Concat  0.107477\n",
       "4        2  Append  0.161898"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_append = append_df()\n",
    "t_concat = concat_df()\n",
    "\n",
    "timings = (pd.DataFrame({\"Append\": t_append, \"Concat\": t_concat})\n",
    "             .stack()\n",
    "             .reset_index()\n",
    "             .rename(columns={0: 'Time (s)',\n",
    "                              'level_1': 'Method'}))\n",
    "timings.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAGcCAYAAAAs6p7QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xtc1GW+B/DPMAw66CaIoBaGgC+xAhUzvMCmhFqaooSI\nx9Bjq7ZQHC9djoYXyCtoq6aUhOtZMmC9lbRtq3XEwjUjQ/cQSpQ21OIFkAAzGWQGnvNHy6wToA8I\nvxmHz/v14vWCZ5555juDfHx+t+enEkIIEBFJsLN0AUR092BgEJE0BgYRSWNgEJE0BgYRSWNgEJE0\nxQOjsLAQERER8Pf3R1hYGPLz85vtt2/fPjz++OMYPnw4IiIikJeXZ3osLy8PM2bMwPDhwzFhwgTs\n3btXqfKJOjehoBs3bohHH31U7NmzRxiNRnHgwAExatQoUVNTY9YvNzdXjBw5UhQVFQkhhDh48KAY\nPny4qK6uFlevXhUBAQHiww8/FEIIcfbsWREQECBOnDih5Fsh6pQUnWHk5uZCrVYjMjISarUa4eHh\ncHFxQU5Ojlm/0tJSzJ8/Hz4+PgCAadOmwc7ODufOncOlS5cwduxYTJo0CQDw4IMPYsSIEfjHP/6h\n5Fsh6pTslXwxnU4Hb29vszZPT0/odDqztqlTp5r9fOrUKdTU1GDAgAFwcnJCUlKS6bGrV68iLy8P\nYWFhHVc4EQFQeB+GXq+HVqs1a9NqtaitrW3xOefPn8eiRYuwaNEiODk5mT127do1REdHw8/PD8HB\nwR1SMxH9m6IzjObCQa/Xw9HRsdn+x48fxwsvvIB58+Zh/vz5Zo+VlJQgJiYGHh4e2LJlS6vqqKqq\nQnV1tVlbfX09bty4AR8fH9jbK/qxEN01FP3L8PLyQkZGhllbcXExQkNDm/R99913sWHDBqxevdq0\nv6LR2bNnsWDBAkydOhVLly5tdR3p6elITk5u9rHs7Gy4u7u3ekyiTkHJPayNR0nS09OFwWAQ+/fv\nF6NHjxZ6vd6s34kTJ8TgwYNFXl5ekzGuXLkiRo0aJXbu3NnmOiorK4VOpzP7OnHihBg4cKAoKSlp\n87hEtk4lhLKXt3/77bdYtWoVzp07Bw8PDyQkJGDw4MGIj4+HSqVCQkIC5s2bh9zcXHTp0qUx1KBS\nqbBt2zacPXsWW7duhVarRWPpKpUKc+bMweLFi9tc14ULFxASEsIZBtEtKB4Y1oqBQXR7PDWciKQx\nMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhI\nGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgOD\niKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQxMIhIGgODiKQx\nMIhIGgODiKQxMIhIGgODiKQxMIhIGgPDRhQUFKCgoMDSZZCNY2DYiMzMTGRmZlq6DLJxDAwbUFBQ\ngDNnzuDMmTOcZVCHYmDYgJtnFpxlUEdSPDAKCwsREREBf39/hIWFIT8/v9l++/btw+OPP47hw4cj\nIiICeXl5rR6DiNqZUNCNGzfEo48+Kvbs2SOMRqM4cOCAGDVqlKipqTHrl5ubK0aOHCmKioqEEEIc\nPHhQDB8+XFRXV0uP0VolJSVi4MCBoqSk5I7GsYSvvvpKTJ48WUyePFl89dVXli6HbJiiM4zc3Fyo\n1WpERkZCrVYjPDwcLi4uyMnJMetXWlqK+fPnw8fHBwAwbdo02NnZ4dy5c9JjdCZ+fn7w9fWFr68v\n/Pz8LF0O2TB7JV9Mp9PB29vbrM3T0xM6nc6sberUqWY/nzp1CjU1NRgwYACysrKkxuhsZs2aZekS\nqBNQNDD0ej20Wq1Zm1arRW1tbYvPOX/+PBYtWoRFixbBycmpTWN0BpxZkBIUDYzm/rD1ej0cHR2b\n7X/8+HG88MILmDdvHubPn9+mMZpTVVWF6upqs7bS0lLp5xN1VooGhpeXFzIyMszaiouLERoa2qTv\nu+++iw0bNmD16tWYNGlSm8ZoSXp6OpKTk1tZPRFZ5ChJenq6MBgMYv/+/WL06NFCr9eb9Ttx4oQY\nPHiwyMvLa/MYt1JZWSl0Op3Z14kTJ+7aoyRESlEJIYSSAfXtt99i1apVOHfuHDw8PJCQkIDBgwcj\nPj4eKpUKCQkJmDdvHnJzc9GlS5fGUINKpcK2bdsQFBTU4hh34sKFCwgJCUF2djbc3d3b460S2RzF\nA8NaMTCIbo+nhhORNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCR\nNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYBCRNAYGEUljYNiIgoIC\nFBQUWLoMsnEMDBuRmZmJzMxMS5dBNo6BYQMKCgpw5swZnDlzhrMM6lAMDBtw88yCswzqSAwMIpLG\nwLABs2bNavZ7ovZmb+kC6M75+fnB19fX9D1RR2Fg2AjOLEgJDAwbwZkFKYH7MGwET9wiJTAwbARP\n3CIlMDBsAE/cIqUwMGwAT9wipTAwiEgaA8MG8MQtUgoPq9oAnrhFSmFg2AjOLEgJ3CQhImkMDBvB\n8zBICQwMG8DzMEgpDAwbwPMwSCkMDCKSxsCwATwPg5TCwCAiaQwMG8B9GKQUxQOjsLAQERER8Pf3\nR1hYGPLz82/ZPy0tDQsXLjRrO3r0KKZMmYJhw4Zh4sSJ+Otf/9qRJRPRvygaGHV1dYiJicH06dOR\nl5eHqKgoxMTEQK/XN+mr1+uxceNGJCUlQaVSmdpra2uxePFiLFy4EKdPn8aaNWuwbNkyXLp0Scm3\nYlW4D4OUomhg5ObmQq1WIzIyEmq1GuHh4XBxcUFOTk6TvrGxsSgpKcHMmTPN2lUqFbp16waDwWD6\nWaPRQK1WK/IerFHjtSS+vr68loQ6lKLXkuh0Onh7e5u1eXp6QqfTNembmJgIV1dXJCcno7Ky0tTe\npUsXJCYmYuHChXj55ZchhMC6devQu3fvDq/fmnFmQUpQdIah1+uh1WrN2rRaLWpra5v0dXV1bXaM\nixcv4sUXX8S6deuQn5+PHTt2YN26dfjmm286pGYi+jdFA6O5cNDr9XB0dJQe48iRI3jggQcwefJk\n2NvbY8yYMRg7diyysrKkx6iqqkJxcbHZV0lJifTzrRGvJSElKLpJ4uXlhYyMDLO24uJihIaGSo/R\npUsX1NXVmbXZ29vD3l7+raSnpyM5OVm6v7VrvJak8Xvux6COougMY+TIkairq0NGRgaMRiMOHDiA\nyspKBAUFSY8xduxY6HQ6HDx4EABw8uRJHDlyBBMnTpQeIyoqCocPHzb7SktLa+3bsRo8D4OUougM\nw8HBATt37sSqVauwefNmeHh4YMeOHejatSvi4+OhUqmQkJBwyzH69OmDlJQUJCYmYv369ejTpw+S\nkpLw4IMPStfh7OwMZ2dnszaNRtOWt0TUqaiEEMLSRViDCxcuICQkBNnZ2XB3d7d0Oa1SUFCAuLg4\nAMD69eu5SUIdhkv02QCu6UlKYWDYCJ6HQUpgYFiQwWBARUVFu4zl5OQEALh8+XK7jAcAvXr14r4d\nMsPAsBCDwYDo6GiUl5dbupQWubm5ISUlhaFBJry8nYikcYZhIRqNBikpKe2ySVJeXo4VK1YAANau\nXQs3N7c7HhPgJgk1xcCwII1Gg759+7brmG5ubu0+JlEjbpIQkTQGBhFJY2AQkTQGBhFJY2AQkTQG\nBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJ\nY2AQkTQGBhFJY2AQkTQGBhFJu+1tBvR6PbKysnDs2DEUFBSguroaKpUKLi4u8PX1xZgxY/Dkk0/C\n0dFRiXqJyIJaDAyj0Yi33noLb7/9Nu6//34EBQXhySefhJOTE+rr61FVVYVvvvkG+/btw2uvvYbZ\ns2fj2WefhYODg5L1E5GCWgyMGTNmICgoCB988AF69+59y0G+//577NmzBzNmzEBWVla7F0lE1qHF\nwHjrrbfg6uoqNUj//v2xbNkyq76xMBHduRZ3et4uLC5fvoz6+nqztva6pycRWSepoyRlZWX4r//6\nL5w9exY3btzArFmzEBwcjODgYBQVFXV0jURkJaQC49VXX0V1dTWcnZ1x8OBBnDt3Dnv37sW4ceOw\nbt26jq6R6K5RUFCAgoICS5fRYaTu3p6bm4sDBw7g3nvvxZEjRxAcHIwhQ4agZ8+emDJlSkfXaBWM\nRiOuXLli6TKadfO+I2vej+Tq6gp7e6l/cnetzMxMAMCGDRssXEnHkPrtaTQa1NfX4/r16zh58qTp\nwygvL+80519cuXIFzz77rKXLuK0VK1ZYuoQWpaamom/fvpYuo8MUFBTgzJkzpu/9/PwsXFH7k9ok\nGT16NOLi4hATEwMHBweMHTsWx44dQ1xcHMaMGdPRNRLdFRpnF7/+3pZIzTDWrFmD119/HRcvXsSb\nb76Jbt26obCwEAEBAVi2bFlH12h1tPcHw05jXTMrIRoAACqVdZ3t32Cogf6fn1i6DGonLQaGEAIq\nlQoA0L17dyxfvtzs8ejo6Fs+x5bZaRxh5/AbS5dBVmbWrFmIi4szfW+LWvzvaPr06Thy5IjUIA0N\nDfjb3/6G8PDwdiuM6G7j5+cHX19f+Pr62uT+C+AWM4ytW7dizZo1WLduHUJCQhAYGIgBAwbA2dkZ\nQghUVVWhqKgIX375JQ4dOoRBgwZh69atStZOZHVsdWbRqMXA6NevH1JTU/HVV18hIyMDcXFxqKqq\nMtvkcHFxwaOPPoo333wTgwcPVqRgImtmqzOLRrfd6Tl48GBTGFy8eBE//vgj7Ozs4OrqetuL0ojI\ntrTqLJr77rsP9913X0fVQkRWzrqOwRGRVbPt83Q7SIPhuqVLuGvws7ItigdGYWEh4uPjcf78efTv\n3x8JCQkYMmRIi/3T0tJw+vRpbNu2zdRWVlaG+Ph4fPnll/jNb36DefPmYfbs2R1at9FoNH2v/+en\nHfpaturmz5DuTq3aJCkrK0Nubi5qa2tRUVHR6herq6tDTEwMpk+fjry8PERFRSEmJgZ6vb5JX71e\nj40bNyIpKanJyWDPPfccBgwYgC+//BK7du1CcnIy/u///q/V9RC1N16tCqCmpgbLly/HoUOHYGdn\nh48++ggbNmxAdXU1kpOT0bNnT6kXy83NhVqtRmRkJAAgPDwcaWlpyMnJwRNPPGHWNzY2Fo6Ojpg5\ncyYqKytN7fn5+bhy5QpefPFFqFQqeHt7Y+/evXB2dpZ9z21y81WW2vvHwk7TrUNfz1Y0GK6bZmS2\nfqUqwKtVAQCbNm1CaWkpDh06hKeeegoA8NJLL2Hp0qVYv349XnvtNakX0+l08Pb2Nmvz9PSETqdr\n0jcxMRGurq5ITk42C4yzZ89iwIAB2LhxIz744AN0794d0dHRmDZtmlQN7cFO042nhlMTvFr1X7Kz\ns/HKK6/A09PT1Obl5YWEhAT8/e9/l34xvV4PrVZr1qbValFbW9ukb0tLBF69ehVffPEFevbsiU8/\n/RQbNmzAmjVrcOrUKek6iDoCr1b9l59//hndu3dv0q5SqVq1I6u5cNDr9a1aU8PBwQFOTk5YsGAB\nAMDf3x8TJkxAdnY2Hn74YakxqqqqUF1dbdZWWloqXQNRZyU1wwgKCkJKSorZor9VVVXYtGkTAgMD\npV/My8sLxcXFZm3FxcUYMGCA9Bienp4wGo0QQpjaGhoazH6+nfT0dDzxxBNmX3PnzpV+PlFzbr6O\nxFavKZEKjBUrVuD777/HqFGjUFtbi/nz5yM4OBhXr15tctn7rYwcORJ1dXXIyMiA0WjEgQMHUFlZ\niaCgIOkxAgMDodVqkZycjPr6epw+fRpHjhzBxIkTpceIiorC4cOHzb7S0tKkn0/UnE59terN3Nzc\nsG/fPnz++efQ6XQwGo3w9vZGYGBgq9a/cHBwwM6dO7Fq1Sps3rwZHh4e2LFjB7p27Yr4+HioVCok\nJCTccowuXbrgnXfewauvvorRo0eje/fuWLlyZasufnN2dm5yVEWj0Ug/n6gltjqzaKQSrZjLCyFg\nMBiatNvC7REvXLiAkJAQZGdnw93dvcnjly9fNq3p2c37SR4lkdRQdw3Xv/sQgO2v6dkZSM0wTp48\niVdffRXff/89GhoaTO2NK2x9/fXXHVagNWow1Fi6hCaseYk+sh1SgbFy5UoMGDAAS5cuRdeuXTu6\nJqvHNSqps5IKjPLycqSkpJidh0FEnY9UYIwfPx45OTmdOjBcXV2Rmppq6TKaVV5ebrofydq1a632\nHreyN/cm6yUVGC+88AJCQ0Px17/+Ff369YOdnfl28h/+8IcOKc6a2Nvb3xU77Nzc3O6KOunuJBUY\ny5cvh0qlgru7O/dhEN3C+++/DwCYOnWqhSvpGFKBkZeXh/T0dJs9GYWovTReQ2KrgSF1DM7DwwN1\ndXUdXQvRXe39999HTU0NampqTDMNWyM1w4iJicGyZcswe/Zs3H///U3WNWjNqd1EturXV6va4ixD\nKjCWLFkCAFi/fn2TxzrjiVtEnZVUYBQVFXV0HUR3vVmzZuGPf/yj6Xtb1OI+jJv3WdTV1d3yi4h+\n2dHp6OgIR0dHm9wcAW4xwxgyZAiOHz8OFxcXDB48uNmrUjvrtSRELbHVmUWjFgPj7bffRo8ePQAA\nu3fvVqwgoruZrc4sGrUYGCdPnoSfnx/s7e0REBCgZE1EZKVa3IfxxhtvoKaGlyYT0b+1GBitWSOT\niDqHWx5WNRgMUkdBbGHFLaL20HjXM1u9jOKWgREcHCw1CI+SEP2iU9/5bNu2baYjJUR0a53hzmct\nBoZKpcKwYcPg4uKiZD1Ed61fX0tii7OMFgODOz07nsFgQEVFxR2PU15e3uz3d6pXr168/QKZaTEw\nwsLC0KVLFyVr6VQMBgOio6Pb9Q8cgGmpvvbg5uaGlJQUhoakWbNmIS4uzvS9LWoxMGxxOkXUkRrv\nfNb4vS2SulqV2p9Go0FKSkq7bJIAMN0U+9drldwJbpK0nq3OLBoxMCxIo9FwwV4bY6szi0bWdZss\nIrJqDAwiksbAICJpDAyidlRQUGC6nsQWMTCI2lFmZqbZGZ+2hoFB1E4aryU5c+aMzc4yGBhE7eTX\n15LYIgYGEUljYBC1k5vP8rTVMz55pidRO+G1JETUKrY6s2jEwCBqR7Y6s2jEfRhEJI2BQUTSGBhE\nJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSFA+MwsJCREREwN/fH2FhYcjPz79l/7S0\nNCxcuLDZxyoqKjB69Gjk5OR0RKlE9CuKBkZdXR1iYmIwffp05OXlISoqCjExMdDr9U366vV6bNy4\nEUlJSVCpVM2Ot3z5cly9erWjyyaif1E0MHJzc6FWqxEZGQm1Wo3w8HC4uLg0O0OIjY1FSUkJZs6c\n2exYe/bsQbdu3dCnT5+OLpuI/kXRwNDpdPD29jZr8/T0hE6na9I3MTER27dvh4uLS5PHiouL8ac/\n/QkJCQm8yzyRghS9vF2v10Or1Zq1abVa1NbWNunr6ura7Bj19fVYunQpVq5ciXvuuadNdVRVVaG6\nutqsrbS0tE1jEXUmigZGc+Gg1+vh6OgoPcYbb7yBBx54AEFBQW2uIz09HcnJyW1+PlFnpWhgeHl5\nISMjw6ytuLgYoaGh0mMcOnQIFRUVOHToEADg2rVrWLJkCWJiYrBgwQKpMaKiojB58mSzttLSUsyd\nO1e6DqLOSNHAGDlyJOrq6pCRkYHIyEhkZWWhsrKyVbOFxqBo9NhjjyE+Ph5jxoyRHsPZ2RnOzs5m\nbRqNRvr5RJ2Vojs9HRwcsHPnTnzwwQcYMWIEMjMzsWPHDnTt2hXx8fFISEho9ZgtHXIlovanEjzM\nAAC4cOECQkJCkJ2dDXd3d0uXQ2SVeGo4UTvizZiJSBpvxkxEUngzZiKSxpsxExHdhIFB1E54M2Yi\nksabMRNRq9jqzKIRA4OoHdnqzKIR92EQkTQGBhFJY2AQkTQGBhFJY2AQkTQeJaFOzWAwoKKiot3G\nMxqNAAB7+/b50+rVq5dVLe7EwKBOy2AwIDo6GuXl5ZYupUVubm5ISUmxmtDgJgkRSeMMgzotjUaD\nlJSUdtskKS8vx4oVKwAAa9euhZub2x2PyU0SIiui0WjQt2/fdh/Xzc2tQ8a1NG6SEJE0BgYRSWNg\nEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0\nBgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYR\nSVM8MAoLCxEREQF/f3+EhYUhPz//lv3T0tKwcOFCs7a8vDzMmDEDw4cPx4QJE7B3796OLJmI/kXR\nwKirq0NMTAymT5+OvLw8REVFISYmBnq9vklfvV6PjRs3IikpCSqVytT+008/4fnnn8fcuXORl5eH\nrVu3YvPmzfj888+VfCtEnZKigZGbmwu1Wo3IyEio1WqEh4fDxcUFOTk5TfrGxsaipKQEM2fONGu/\ndOkSxo4di0mTJgEAHnzwQYwYMQL/+Mc/FHkPRJ2ZooGh0+ng7e1t1ubp6QmdTtekb2JiIrZv3w4X\nFxez9kGDBiEpKcn089WrV5GXl4cHHnigY4omIhNFA0Ov10Or1Zq1abVa1NbWNunr6up62/GuXbuG\n6Oho+Pn5ITg4uN3qJKLm2Sv5Ys2Fg16vh6OjY6vHKikpQUxMDDw8PLBly5ZWPbeqqgrV1dVmbaWl\npa2ugaizUTQwvLy8kJGRYdZWXFyM0NDQVo1z9uxZLFiwAFOnTsXSpUtbXUd6ejqSk5Nb/Tyizk7R\nwBg5ciTq6uqQkZGByMhIZGVlobKyEkFBQdJjVFRUYMGCBfjd736H+fPnt6mOqKgoTJ482ayttLQU\nc+fObdN4RJ2FooHh4OCAnTt3YtWqVdi8eTM8PDywY8cOdO3aFfHx8VCpVEhISLjlGO+++y6qqqrw\n5ptv4o033gAAqFQqzJkzB4sXL5aqw9nZGc7OzmZtGo2mTe+JqDNRCSGEpYuwBhcuXEBISAiys7Ph\n7u5u6XKoBUajEVeuXLF0Gc0qLy/HihUrAABr166Fm5ubhStqnqurK+zt2zZXUHSGQXSnrly5gmef\nfdbSZdxWY3BYo9TUVPTt27dNz+W1JEQkjTMMumvdE9QXdo7W9U9YNPyyha+yU92mp7Iaaoz46fjl\nOx7Huj5tolawc7SHujt3ViuJmyREJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhE\nJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhEJI2BQUTSGBhEJI1retJdq77GaOkS\n7hrt9VkxMOiuYjT++x/+tXZYBbszuvkzbC1ukhCRNM4w6K5y8y3+fhPUF2oruy+JtaqvMZpmZG29\nTSLAwKC7mJr3JVEcN0mISBpnGHTXarDCoyTWfKvE9sDAoLtWe9wrlFqHmyREJI0zDLqruLq6IjU1\n1dJlNKu8vBwrVqwAAKxduxZubm4Wrqh5rq6ubX4uA4PuKvb29ujbt6+ly7gtNze3u6LO1uImCRFJ\nY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTQGBhFJY2AQkTRefEadmsFg\nQEVFRbuMVV5e3uz3d6JXr17QaKxnGUIGBnVaBoMB0dHR7fbHfbPGy9zvlJubG1JSUqwmNBTfJCks\nLERERAT8/f0RFhaG/Pz8W/ZPS0vDwoUL72gMImofis4w6urqEBMTg+eeew7Tp09HVlYWYmJikJ2d\nDa1Wa9ZXr9dj+/bt+NOf/oQJEya0aQyiW9FoNEhJSWm3TRLg3zcJupOl/G/WqTdJcnNzoVarERkZ\nCQAIDw9HWloacnJy8MQTT5j1jY2NhaOjI2bOnInKyso2jUF0OxqNxiYXuukoim6S6HQ6eHt7m7V5\nenpCp9M16ZuYmIjt27fDxcWlzWMQUftSNDD0en2TzQatVova2tomfVtad7A1YxBR+1J0k6S5P2y9\nXg9HR0dFx6iqqkJ1dbVZ28WLFwEApaWl0uMQ2bI+ffo02RejaGB4eXkhIyPDrK24uBihoaGKjpGe\nno7k5ORmH3v66aelxyGyZdnZ2XB3dzdrUzQwRo4cibq6OmRkZCAyMhJZWVmorKxEUFCQomNERUVh\n8uTJZm11dXW4dOkSvLy8oFarpceyFiUlJZg7dy7S0tLQr18/S5fTadnS76FPnz5N2hQNDAcHB+zc\nuROrVq3C5s2b4eHhgR07dqBr166Ij4+HSqVCQkJCm8eQ5ezsDGdn5ybtPj4+rX1LVsNgMAD45Zf8\n6/8VSDm2/ntQCSGEpYugO1dcXIwnnngChw8fhqenp6XL6bRs/ffAi8+ISBoDg4ikMTBshJOTE2Jj\nY+Hk5GTpUjo1W/89cB8GEUnjDIOIpDEwiEgaA4OIpDEwiEgaA4OIpDEwiEgaA4PISl24cMHSJTTB\nwLACL7/8Mnx9fXHlyhVLl2KmoqICgwYNsnQZijl27Bjmzp2LESNGYMSIEZg/fz7OnDljkVqOHj2K\nJUuWWOS1b4WBYWE//fQTjh07hokTJ+LPf/6zpcsxI4SASqWydBmK2LdvH+Li4vDMM8/gxIkT+Pvf\n/47AwED853/+J7777jvF66muroY1nlPJ+5JYWFZWFh555BE8/fTTiI2NxXPPPQd7e3skJyejtLQU\nOp0OX3/9NXx8fPDqq6/Cx8cHJ0+exJYtW9C/f38cPnwYbm5uWLx4MSZOnAgAuHz5MlavXo3Tp0/D\nyckJv//97/HUU08BAGbPno2HH34Yn376Kf75z3/ioYceQlJSEu69914IIbB161bs27cPdnZ2+I//\n+A9LfjSpdh8rAAAKH0lEQVSKqa2tRVJSEjZv3owxY8YAANRqNZ555hlUVVXhu+++g5OTE9avX4/P\nPvsMWq0WkyZNwuLFi6HRaPDKK6+gW7du+Prrr1FYWAhvb2+sWbMGDzzwAAAgMzMT//M//4Pq6mr4\n+flh9erV6NevHwoLC7Fp0yacO3cO169fx7Bhw7Bp0yZcvHgRCQkJMBqNCAoKwvHjxy358ZgTZFGT\nJ08Wn3zyiRBCiClTpogPPvhACCHE9u3bxUMPPSQ+/fRTYTAYxLZt20RwcLAwGAziiy++ED4+PmLb\ntm3CYDCInJwc4efnJ86fPy/q6+tFaGio2Lx5szAajaKoqEj89re/FV988YUQQoioqCgxfvx4cfHi\nRXHt2jXx9NNPi1WrVgkhhMjMzBQTJkwQJSUl4ueffxYxMTFi0KBBFvlclPTZZ5+JoUOHivr6+hb7\nREZGipdeeknU1NSIsrIyMX36dPGHP/xBCCHEsmXLREBAgPjmm2/EjRs3xJIlS8S8efOEEELk5OSI\ngIAAUVBQIBoaGsSmTZvEzJkzhRBCjB8/XqSnpwshhKiurhYRERHi9ddfF0II8d5774nw8PCOfNtt\nwk0SCzp9+jSuXbtm+l9t5syZSE9PNz0+btw4jBkzBvb29nj++edx/fp1nD59GsAv96uIjY2Fvb09\nHn30UQQFBeHQoUMoKChAaWkplixZArVaDR8fH8yYMQP79u0zjRsaGop7770X3bt3x7hx4/DDDz8A\nAP72t78hKioK7u7u6NatG1566SUFPw3Lqaqqwj333AM7u+b/HEpKSpCfn48VK1ZAq9XCzc0NixYt\nwnvvvWfq89hjj2HgwIFwcHDApEmTzD7TsLAw+Pr6QqVSITY21nRXtF27dmHWrFnQ6/W4fPkynJ2d\nUVZW1vFv+A5wk8SC9u3bh6qqKvz2t78F8MtNcK5evYqzZ88CADw8PEx97ezs4ObmhoqKCvTq1Qvu\n7u5m+xf69u2LK1eu4NKlS7h27RoCAgIA/LIfoqGhAQ899JCp782rjWk0GjQ0NAD45X6gNy/LZosr\nRjWnV69euHr1Kurr65ssz3jt2jWUlZVBq9WiR48epvZ7770XP/74I+rr6wGYf6b29vamz/TXO467\ndu1q+l3k5+djwYIFqKmpwcCBA/HTTz+hZ8+eHfY+2wMDw0J+/vlnHD58GG+//bbZ2o/r1q1Deno6\n7rvvPrMVzOvr601/0Eajscn9QC9evIihQ4fC1dUVffr0wdGjR02P3XwjqFvp3bu3afV0AFb/v117\n8ff3h0ajwbFjxxAcHGz2WFxcHK5fvw69Xo+rV6+aQqOkpAQ9evS47fqvvXv3Nvscr1+/jjfeeAOz\nZ8/GsmXL8Oc//xl+fn6m1xJWuKPzZtwksZCsrCz0798fQ4cOhYuLi+krPDwcH374IaqqqvDxxx/j\n1KlTMBgM2L59O3r27Al/f38AwKVLl7B7927U19fj008/xRdffIEnn3wSQ4cORdeuXbFr1y4YjUaU\nlZVh7ty5Zps6LXnqqafwzjvvQKfToaamBps3b+7oj8EqODg4YMmSJVi5ciVycnJQX1+P69evIzk5\nGZ9//jni4uIwatQorF+/HjU1NSgrK8P27dulVqqfMmUKsrKyUFRUBKPRiJSUFOTn50Ov1wOAaS3a\nnJwcHD582HSrRQcHB1y/fr3j3nQbcYZhIfv378eUKVOatI8ePRo9e/bE/v378cgjj2Dr1q0oLCzE\nkCFDkJqaatoM6dOnD77++muMGjUKvXv3RnJysmmm8tZbb2Ht2rVITU2FRqPBk08+ieeffx4AbnmY\ndNq0aaioqMCcOXNgNBoxd+5cHD58uAPevfWZNWsWevTogeTkZLz88suws7PDkCFDkJ6ejgEDBuC1\n117D2rVrERISApVKhalTp0qdJzFy5Ei8/PLLWLRoESorKzFs2DBs3rwZvXv3xnPPPYc5c+ZAo9Eg\nICAAL774It59910AwCOPPAIhBB555BF89tlncHBw6OiPQAoX0LFSycnJKCsrw5o1a5o8dvLkSaxc\nuRIfffSRBSqjzoybJEQkjYFBRNK4SUJE0jjDICJpDAwiksbAICJpDAwiksbAoFYbNGgQBg0ahPPn\nzzd5rKCgAIMGDcKcOXOkx/v4449Np7ofPHgQQUFB7VYrADzzzDNITk5u1zE7KwYGtYlGo8GRI0ea\ntH/88cctXvXZnEuXLmHhwoVWeRo0NcXAoDYJCAhoMTCGDh0qPU5DQ0OnWdXLFjAwqE3Gjx+PwsJC\nsysxv/32W9PKUTfT6XSYN28ehg4dinHjxuH11183XRY+btw4AMCkSZOQlZVlek5qaiqCgoLg7++P\nV155BXV1dabHjh8/joiICAwdOhTjx4/Hnj17zF7vvffew7hx4+Dv74/Vq1ebXovuHAOD2sTd3R0+\nPj5ms4yPP/4Y48ePN5sx1NXVYf78+fDx8cFf/vIXrF+/Hh999BG2bNkC4JeL8IQQyMzMxKRJkwD8\nsobEV199hd27dyM5ORmHDh3C/v37AQBffvkloqOjMXHiRLz//vuIjo5GYmKi6SK5EydOYNWqVfj9\n73+P9957Dw0NDTh58qRSH4vNY2BQm4WEhCA7O9v08//+7/9iwoQJZn3+8pe/oFu3bvjv//5v3H//\n/QgICMDy5cuxe/duCCFMC8Y4OTmZrshUq9VITEyEl5cXAgMDERgYiMLCQgDAO++8g8ceewy/+93v\n4OHhgfDwcERFRSE1NRUAsHfvXjz++OOIiIiAp6cnVq5cib59+yrxcXQKDAxqs/Hjx+PkyZP4+eef\n8cMPP6CsrAwjRoww6/Pdd99Bp9PB39/f9BUbGwuj0djifTd69OiB7t27m36+5557UFtbaxpv8ODB\nZv2HDRuG4uJi0+M3ry6mVqtNi/HSneN6GNRmPj4+uO+++/DJJ5+gvLwcISEhTY6Q1NfX4+GHH8a6\ndeuaPL9v377Nrup1q6MsXbp0adLW0NBgWnhGpVI1WbXK3p7/zNsLZxh0R0JCQnD06NFmN0cAwNvb\nGz/88AP69OmDfv36oV+/frh8+TJee+21Nh0h8fT0RH5+vlnbqVOn0L9/fwDAwIEDzR4XQqCoqKj1\nb4yaxcCgOzJ+/Hjk5OSguLgYgYGBTR5vXMZu2bJlOH/+PPLy8rBixQpoNBo4ODjA0dERAFBUVISa\nmprbvt68efPwySef4I9//CN++OEHHDhwAHv27DGdKDZ79mwcPXoUu3fvxvfff48NGzaYrVNKd4aB\nQa1286xg6NCh6Natm+l2CL+m1Wqxa9cuVFdXIyIiAosWLUJgYCDWrl0L4JednWFhYVi6dKnZrRBa\n8uCDD2Lr1q14//33MWXKFOzatQtxcXGIiIgw1bNlyxZkZmZi2rRp+PHHH02HbunOcT0MIpLGGQYR\nSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSWNgEJE0BgYRSft/VlNAxNVVclwAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1103e48d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4, 6))\n",
    "sns.boxplot(x='Method', y='Time (s)', data=timings)\n",
    "sns.despine()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datatypes\n",
    "\n",
    "The pandas type system essentially [NumPy's](http://docs.scipy.org/doc/numpy/reference/arrays.scalars.html) with a few extensions (`categorical`, `datetime64` with timezone, `timedelta64`).\n",
    "An advantage of the DataFrame over a 2-dimensional NumPy array is that the DataFrame can have columns of various types within a single table.\n",
    "That said, each column should have a specific dtype; you don't want to be mixing bools with ints with strings within a single column.\n",
    "For one thing, this is slow.\n",
    "It forces the column to be have an `object` dtype (the fallback python-object container type), which means you don't get any of the type-specific optimizations in pandas or NumPy.\n",
    "For another, it means you're probably violating the maxims of tidy data, which we'll discuss next time.\n",
    "\n",
    "When should you have `object` columns?\n",
    "There are a few places where the NumPy / pandas type system isn't as rich as you might like.\n",
    "There's no integer NA (at the moment anyway), so if you have any missing values, represented by `NaN`, your otherwise integer column will be floats.\n",
    "There's also no `date` dtype (distinct from `datetime`).\n",
    "Consider the needs of your application: can you treat an integer `1` as `1.0`?\n",
    "Can you treat `date(2016, 1, 1)` as `datetime(2016, 1, 1, 0, 0)`?\n",
    "In my experience, this is rarely a problem other than when writing to something with a stricter schema like a database.\n",
    "But at that point it's fine to cast to one of the less performant types, since you're just not doing numeric operations anymore.\n",
    "\n",
    "The last case of `object` dtype data is text data.\n",
    "Pandas doesn't have any fixed-width string dtypes, so you're stuck with python objects.\n",
    "There is an important exception here, and that's low-cardinality text data, for which you'll want to use the `category` dtype (see below).\n",
    "\n",
    "If you have object data (either strings or python objects) that needs to be converted, checkout the [`to_numeric`](http://pandas.pydata.org/pandas-docs/version/0.18.0/generated/pandas.to_numeric.html), [`to_datetime`](http://pandas.pydata.org/pandas-docs/version/0.18.0/generated/pandas.to_datetime.html) and [`to_timedelta`](http://pandas.pydata.org/pandas-docs/version/0.18.0/generated/pandas.to_timedelta.html) methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iteration, Apply, And Vectorization\n",
    "\n",
    "We know that [\"Python is slow\"](https://jakevdp.github.io/blog/2014/05/09/why-python-is-slow/) (scare  quotes since that statement is too broad to be meaningful).\n",
    "There are various steps that can be taken to improve your code's performance from relatively simple changes, to rewriting your code in a lower-level language, to trying to parallelize it.\n",
    "And while you might have many options, there's typically an order you would proceed in.\n",
    "\n",
    "First (and I know it's cliché to say so, but still) benchmark your code.\n",
    "Make sure you actually need to spend time optimizing it.\n",
    "There are [many](https://github.com/nvdv/vprof) [options](https://jiffyclub.github.io/snakeviz/) [for](https://github.com/rkern/line_profiler) [benchmarking](https://docs.python.org/3.5/library/timeit.html) and visualizing where things are slow.\n",
    "\n",
    "Second, consider your algorithm.\n",
    "Make sure you aren't doing more work than you need to.\n",
    "A common one I see is doing a full sort on an array, just to select the `N` largest or smallest items.\n",
    "Pandas has methods for that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"data/627361791_T_ONTIME.csv\")\n",
    "delays = df['DEP_DELAY']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62914     1461.0\n",
       "455195    1482.0\n",
       "215520    1496.0\n",
       "454520    1500.0\n",
       "271107    1560.0\n",
       "Name: DEP_DELAY, dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Select the 5 largest delays\n",
    "delays.nlargest(5).sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "307517   -112.0\n",
       "40118     -85.0\n",
       "36065     -46.0\n",
       "86280     -44.0\n",
       "27749     -42.0\n",
       "Name: DEP_DELAY, dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delays.nsmallest(5).sort_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We follow up the `nlargest` or `nsmallest` with a sort (the result of `nlargest/smallest` is unordered), but it's much easier to sort 5 items that 500,000. The timings bear this out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 66.4 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit delays.sort_values().tail(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 12.8 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit delays.nlargest(5).sort_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Use the right algorithm\" is easy to say, but harder to apply in practice since you have to actually figure out the best algorithm to use.\n",
    "That one comes down to experience.\n",
    "\n",
    "Assuming you're at a spot that needs optimizing, and you've got the correct algorithm, *and* there isn't a readily available optimized version of what you need in pandas/numpy/scipy/scikit-learn/statsmodels/..., then what?\n",
    "\n",
    "The first place to turn is probably a vectorized NumPy implementation.\n",
    "Vectorization here means operating directly on arrays, rather than looping over lists scalars.\n",
    "This is generally much less work than rewriting it in something like Cython, and you can get pretty good results just by making *effective* use of NumPy and pandas.\n",
    "While not every operation can be vectorized, many can.\n",
    "\n",
    "Let's work through an example calculating the [Great-circle distance](https://en.wikipedia.org/wiki/Great-circle_distance) between airports.\n",
    "Grab the table of airport latitudes and longitudes from the [BTS website](http://www.transtats.bts.gov/DL_SelectFields.asp?Table_ID=288&DB_Short_Name=Aviation%20Support%20Tables) and extract it to a CSV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import requests\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'os' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-4eae145ee563>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     16\u001b[0m }\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexists\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/airports.csv.zip'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     19\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'url_4.txt'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m         \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'os' is not defined"
     ]
    }
   ],
   "source": [
    "headers = {\n",
    "    'Pragma': 'no-cache',\n",
    "    'Origin': 'http://www.transtats.bts.gov',\n",
    "    'Accept-Encoding': 'gzip, deflate',\n",
    "    'Accept-Language': 'en-US,en;q=0.8',\n",
    "    'Upgrade-Insecure-Requests': '1',\n",
    "    'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_11_4) AppleWebKit/537.36'\\\n",
    "                  '(KHTML, like Gecko) Chrome/51.0.2704.103 Safari/537.36',\n",
    "    'Content-Type': 'application/x-www-form-urlencoded',\n",
    "    'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8',\n",
    "    'Cache-Control': 'no-cache',\n",
    "    'Referer': 'http://www.transtats.bts.gov/DL_SelectFields.asp?Table_ID=288&DB_Short_'\n",
    "               'Name=Aviation%20Support%20Tables',\n",
    "    'Connection': 'keep-alive',\n",
    "    'DNT': '1',\n",
    "}\n",
    "\n",
    "if not os.path.exists('data/airports.csv.zip'):\n",
    "    with open('url_4.txt') as f:\n",
    "        data = f.read().strip()\n",
    "\n",
    "    r = requests.post('http://www.transtats.bts.gov/DownLoad_Table.asp?Table_ID=288&Has'\n",
    "                      '_Group=0&Is_Zipped=0', data=data, headers=headers)\n",
    "\n",
    "    with open('data/airports.csv.zip', 'wb') as f:\n",
    "        f.write(r.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zf = zipfile.ZipFile('data/airports.csv.zip')\n",
    "fp = zf.extract(zf.filelist[0], path='data')\n",
    "airports = pd.read_csv(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATITUDE</th>\n",
       "      <th>LONGITUDE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AIRPORT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8F3</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>58.457500</td>\n",
       "      <td>-154.023333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>60.482222</td>\n",
       "      <td>-146.582222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>63.541667</td>\n",
       "      <td>-150.993889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>59.331667</td>\n",
       "      <td>-135.896667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          LATITUDE   LONGITUDE\n",
       "AIRPORT                       \n",
       "8F3      33.623889 -101.240833\n",
       "A03      58.457500 -154.023333\n",
       "A09      60.482222 -146.582222\n",
       "A18      63.541667 -150.993889\n",
       "A24      59.331667 -135.896667"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coord = (pd.read_csv(fp, index_col=['AIRPORT'],\n",
    "                     usecols=['AIRPORT', 'LATITUDE', 'LONGITUDE'])\n",
    "           .groupby(level=0).first()\n",
    "           .dropna()\n",
    "           .sample(n=500, random_state=42)\n",
    "           .sort_index())\n",
    "\n",
    "coord.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For whatever reason, suppose we're interested in all the pairwise distances (I've limited it to just a sample of 500 airports to make this manageable.\n",
    "In the real world you *probably* don't need *all* the pairwise distances and would be better off with a [tree](http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KDTree.html). Remember: think about what you actually need, and find the right algorithm for that).\n",
    "\n",
    "MultiIndexes have an alternative `from_product` constructor for getting the [Cartesian product](https://en.wikipedia.org/wiki/Cartesian_product) of the arrays you pass in.\n",
    "We'll give it `coords.index` twice (to get its Cartesian product with itself).\n",
    "That gives a MultiIndex of all the combination.\n",
    "With some minor reshaping of `coords` we'll have a DataFrame with all the latitude/longitude pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>LATITUDE_1</th>\n",
       "      <th>LONGITUDE_1</th>\n",
       "      <th>LATITUDE_2</th>\n",
       "      <th>LONGITUDE_2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>origin</th>\n",
       "      <th>dest</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">8F3</th>\n",
       "      <th>8F3</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>58.457500</td>\n",
       "      <td>-154.023333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>60.482222</td>\n",
       "      <td>-146.582222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>63.541667</td>\n",
       "      <td>-150.993889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>59.331667</td>\n",
       "      <td>-135.896667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             LATITUDE_1  LONGITUDE_1  LATITUDE_2  LONGITUDE_2\n",
       "origin dest                                                  \n",
       "8F3    8F3    33.623889  -101.240833   33.623889  -101.240833\n",
       "       A03    33.623889  -101.240833   58.457500  -154.023333\n",
       "       A09    33.623889  -101.240833   60.482222  -146.582222\n",
       "       A18    33.623889  -101.240833   63.541667  -150.993889\n",
       "       A24    33.623889  -101.240833   59.331667  -135.896667"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx = pd.MultiIndex.from_product([coord.index, coord.index],\n",
    "                                 names=['origin', 'dest'])\n",
    "\n",
    "pairs = pd.concat([coord.add_suffix('_1').reindex(idx, level='origin'),\n",
    "                   coord.add_suffix('_2').reindex(idx, level='dest')],\n",
    "                  axis=1)\n",
    "pairs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "idx = idx[idx.get_level_values(0) <= idx.get_level_values(1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "125250"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(idx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll break that down a bit, but don't lose sight of the real target: our great-circle distance calculation.\n",
    "\n",
    "The `add_suffix` (and `add_prefix`) method is handy for quickly renaming the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>LATITUDE_1</th>\n",
       "      <th>LONGITUDE_1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AIRPORT</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8F3</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>58.457500</td>\n",
       "      <td>-154.023333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>60.482222</td>\n",
       "      <td>-146.582222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>63.541667</td>\n",
       "      <td>-150.993889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>59.331667</td>\n",
       "      <td>-135.896667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         LATITUDE_1  LONGITUDE_1\n",
       "AIRPORT                         \n",
       "8F3       33.623889  -101.240833\n",
       "A03       58.457500  -154.023333\n",
       "A09       60.482222  -146.582222\n",
       "A18       63.541667  -150.993889\n",
       "A24       59.331667  -135.896667"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coord.add_suffix('_1').head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively you could use the more general `.rename` like `coord.rename(columns=lambda x: x + '_1')`.\n",
    "\n",
    "Next, we have the `reindex`.\n",
    "Like I mentioned in the prior chapter, indexes are crucial to pandas.\n",
    "`.reindex` is all about aligning a Series or DataFrame to a given index.\n",
    "In this case we use `.reindex` to align our original DataFrame to the new\n",
    "MultiIndex of combinations.\n",
    "By default, the output will have the original value if that index label was already present, and `NaN` otherwise.\n",
    "If we just called `coord.reindex(idx)`, with no additional arguments, we'd get a DataFrame of all `NaN`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>LATITUDE</th>\n",
       "      <th>LONGITUDE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>origin</th>\n",
       "      <th>dest</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">8F3</th>\n",
       "      <th>8F3</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             LATITUDE  LONGITUDE\n",
       "origin dest                     \n",
       "8F3    8F3        NaN        NaN\n",
       "       A03        NaN        NaN\n",
       "       A09        NaN        NaN\n",
       "       A18        NaN        NaN\n",
       "       A24        NaN        NaN"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coord.reindex(idx).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's because there weren't any values of `idx` that were in `coord.index`,\n",
    "which makes sense since `coord.index` is just a regular one-level Index, while `idx` is a MultiIndex.\n",
    "We use the `level` keyword to handle the transition from the original single-level Index, to the two-leveled `idx`.\n",
    "\n",
    "> `level` : int or name\n",
    ">\n",
    "  Broadcast across a level, matching Index values on the\n",
    "    passed MultiIndex level\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>LATITUDE</th>\n",
       "      <th>LONGITUDE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>origin</th>\n",
       "      <th>dest</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">8F3</th>\n",
       "      <th>8F3</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>58.457500</td>\n",
       "      <td>-154.023333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>60.482222</td>\n",
       "      <td>-146.582222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>63.541667</td>\n",
       "      <td>-150.993889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>59.331667</td>\n",
       "      <td>-135.896667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              LATITUDE   LONGITUDE\n",
       "origin dest                       \n",
       "8F3    8F3   33.623889 -101.240833\n",
       "       A03   58.457500 -154.023333\n",
       "       A09   60.482222 -146.582222\n",
       "       A18   63.541667 -150.993889\n",
       "       A24   59.331667 -135.896667"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coord.reindex(idx, level='dest').head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you ever need to do an operation that mixes regular single-level indexes with Multilevel Indexes, look for a level keyword argument.\n",
    "For example, all the arithmatic methods (`.mul`, `.add`, etc.) have them.\n",
    "\n",
    "This is a bit wasteful since the distance from airport `A` to `B` is the same as `B` to `A`.\n",
    "We could easily fix this with a `idx = idx[idx.get_level_values(0) <= idx.get_level_values(1)]`, but we'll ignore that for now.\n",
    "\n",
    "\n",
    "Quick tangent, I got some... let's say skepticism, on my last piece about the value of indexes.\n",
    "Here's an alternative version for the skeptics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from itertools import product, chain\n",
    "coord2 = coord.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"0\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>LATITUDE_1</th>\n",
       "      <th>LONGITUDE_1</th>\n",
       "      <th>LATITUDE_1</th>\n",
       "      <th>LONGITUDE_2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>origin</th>\n",
       "      <th>dest</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">8F3</th>\n",
       "      <th>8F3</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A03</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>58.457500</td>\n",
       "      <td>-154.023333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A09</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>60.482222</td>\n",
       "      <td>-146.582222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A18</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>63.541667</td>\n",
       "      <td>-150.993889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A24</th>\n",
       "      <td>33.623889</td>\n",
       "      <td>-101.240833</td>\n",
       "      <td>59.331667</td>\n",
       "      <td>-135.896667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             LATITUDE_1  LONGITUDE_1  LATITUDE_1  LONGITUDE_2\n",
       "origin dest                                                  \n",
       "8F3    8F3    33.623889  -101.240833   33.623889  -101.240833\n",
       "       A03    33.623889  -101.240833   58.457500  -154.023333\n",
       "       A09    33.623889  -101.240833   60.482222  -146.582222\n",
       "       A18    33.623889  -101.240833   63.541667  -150.993889\n",
       "       A24    33.623889  -101.240833   59.331667  -135.896667"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = product(coord2.add_suffix('_1').itertuples(index=False),\n",
    "            coord2.add_suffix('_2').itertuples(index=False))\n",
    "y = [list(chain.from_iterable(z)) for z in x]\n",
    "\n",
    "df2 = (pd.DataFrame(y, columns=['origin', 'LATITUDE_1', 'LONGITUDE_1',\n",
    "                                'dest', 'LATITUDE_1', 'LONGITUDE_2'])\n",
    "       .set_index(['origin', 'dest']))\n",
    "df2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's also readable (it's Python after all), though a bit slower.\n",
    "To me the `.reindex` method seems more natural.\n",
    "My thought process was, \"I need all the combinations of origin & destination (`MultiIndex.from_product`).\n",
    "Now I need to align this original DataFrame to this new MultiIndex (`coords.reindex`).\"\n",
    "\n",
    "With that diversion out of the way, let's turn back to our great-circle distance calculation.\n",
    "Our first implementation is pure python.\n",
    "The algorithm itself isn't too important, all that matters is that we're doing math operations on scalars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def gcd_py(lat1, lng1, lat2, lng2):\n",
    "    '''\n",
    "    Calculate great circle distance between two points.\n",
    "    http://www.johndcook.com/blog/python_longitude_latitude/\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lat1, lng1, lat2, lng2: float\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    distance:\n",
    "      distance from ``(lat1, lng1)`` to ``(lat2, lng2)`` in kilometers.\n",
    "    '''\n",
    "    # python2 users will have to use ascii identifiers (or upgrade)\n",
    "    degrees_to_radians = math.pi / 180.0\n",
    "    ϕ1 = (90 - lat1) * degrees_to_radians\n",
    "    ϕ2 = (90 - lat2) * degrees_to_radians\n",
    "    \n",
    "    θ1 = lng1 * degrees_to_radians\n",
    "    θ2 = lng2 * degrees_to_radians\n",
    "    \n",
    "    cos = (math.sin(ϕ1) * math.sin(ϕ2) * math.cos(θ1 - θ2) +\n",
    "           math.cos(ϕ1) * math.cos(ϕ2))\n",
    "    # round to avoid precision issues on identical points causing ValueErrors\n",
    "    cos = round(cos, 8)\n",
    "    arc = math.acos(cos)\n",
    "    return arc * 6373  # radius of earth, in kilometers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second implementation uses NumPy.\n",
    "Aside from numpy having a builtin `deg2rad` convenience function (which is probably a bit slower than multiplying by a constant $\\frac{\\pi}{180}$), basically all we've done is swap the `math` prefix for `np`.\n",
    "Thanks to NumPy's broadcasting, we can write code that works on scalars or arrays of conformable shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gcd_vec(lat1, lng1, lat2, lng2):\n",
    "    '''\n",
    "    Calculate great circle distance.\n",
    "    http://www.johndcook.com/blog/python_longitude_latitude/\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lat1, lng1, lat2, lng2: float or array of float\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    distance:\n",
    "      distance from ``(lat1, lng1)`` to ``(lat2, lng2)`` in kilometers.\n",
    "    '''\n",
    "    # python2 users will have to use ascii identifiers\n",
    "    ϕ1 = np.deg2rad(90 - lat1)\n",
    "    ϕ2 = np.deg2rad(90 - lat2)\n",
    "    \n",
    "    θ1 = np.deg2rad(lng1)\n",
    "    θ2 = np.deg2rad(lng2)\n",
    "    \n",
    "    cos = (np.sin(ϕ1) * np.sin(ϕ2) * np.cos(θ1 - θ2) +\n",
    "           np.cos(ϕ1) * np.cos(ϕ2))\n",
    "    arc = np.arccos(cos)\n",
    "    return arc * 6373"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the python version on our DataFrame, we can either iterate..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 916 ms, sys: 14 ms, total: 930 ms\n",
      "Wall time: 931 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "origin  dest\n",
       "8F3     8F3         0.000000\n",
       "        A03      4744.967448\n",
       "        A09      4407.533212\n",
       "        A18      4744.593127\n",
       "        A24      3820.092688\n",
       "                    ...     \n",
       "ZXV     ZAZ      6748.190727\n",
       "        ZBF      1736.084217\n",
       "        ZBX       832.642824\n",
       "        ZKB     12843.096516\n",
       "        ZXV         0.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "pd.Series([gcd_py(*x) for x in pairs.itertuples(index=False)],\n",
    "          index=pairs.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or use `DataFrame.apply`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.2 s, sys: 105 ms, total: 19.3 s\n",
      "Wall time: 19.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "r = pairs.apply(lambda x: gcd_py(x['LATITUDE_1'], x['LONGITUDE_1'],\n",
    "                                 x['LATITUDE_2'], x['LONGITUDE_2']), axis=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But as you can see, you don't want to use apply, especially with `axis=1` (calling the function on each row).  It's doing a lot more work handling dtypes in the background, and trying to infer the correct output shape that are pure overhead in this case. On top of that, it has to essentially use a for loop internally.\n",
    "\n",
    "You *rarely* want to use `DataFrame.apply` and almost never should use it with `axis=1`. Better to write functions that take arrays, and pass those in directly. Like we did with the vectorized version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 49.4 ms, sys: 10.5 ms, total: 60 ms\n",
      "Wall time: 49.3 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "r = gcd_vec(pairs['LATITUDE_1'], pairs['LONGITUDE_1'],\n",
    "            pairs['LATITUDE_2'], pairs['LONGITUDE_2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "origin  dest\n",
       "8F3     8F3        0.000000\n",
       "        A03     4744.967484\n",
       "        A09     4407.533240\n",
       "        A18     4744.593111\n",
       "        A24     3820.092639\n",
       "dtype: float64"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I try not to use the word \"easy\" when teaching, but that optimization was easy right?\n",
    "Why then, do I come across uses of `apply`, in my code and others', even when the vectorized version is available?\n",
    "The difficulty lies in knowing about broadcasting, and seeing where to apply it.\n",
    "\n",
    "For example, the README for [lifetimes](https://github.com/CamDavidsonPilon/lifetimes) (by Cam Davidson Pilon, also author of [Bayesian Methods for Hackers](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers), [lifelines](https://github.com/CamDavidsonPilon/lifelines), and [Data Origami](https://dataorigami.net)) used to have an example of passing [this method](https://github.com/CamDavidsonPilon/lifetimes/blob/5b4f7de0720413b6951ac0a4b0082bd50255a231/lifetimes/estimation.py#L249) into a `DataFrame.apply`.\n",
    "\n",
    "```python\n",
    "data.apply(lambda r: bgf.conditional_expected_number_of_purchases_up_to_time(\n",
    "    t, r['frequency'], r['recency'], r['T']), axis=1\n",
    ")\n",
    "```\n",
    "\n",
    "If you look at the function [I linked to](https://github.com/CamDavidsonPilon/lifetimes/blob/5b4f7de0720413b6951ac0a4b0082bd50255a231/lifetimes/estimation.py#L249), it's doing a fairly complicated computation involving a negative log likelihood and the Gamma function from `scipy.special`.\n",
    "But crucially, it was already vectorized.\n",
    "We were able to change the example to just pass the arrays (Series in this case) into the function, rather than applying the function to each row.\n",
    "\n",
    "```python\n",
    "bgf.conditional_expected_number_of_purchases_up_to_time(\n",
    "    t, data['frequency'], data['recency'], data['T']\n",
    ")\n",
    "```\n",
    "\n",
    "This got us another 30x speedup on the example dataset.\n",
    "I bring this up because it's very natural to have to translate an equation to code and think, \"Ok now I need to apply this function to each row\", so you reach for `DataFrame.apply`.\n",
    "See if you can just pass in the NumPy array or Series itself instead.\n",
    "\n",
    "Not all operations this easy to vectorize.\n",
    "Some operations are iterative by nature, and rely on the results of surrounding computations to proceed. In cases like this you can hope that one of the scientific python libraries has implemented it efficiently for you, or write your own solution using Numba / C / Cython / Fortran.\n",
    "\n",
    "Other examples take a bit more thought or knowledge to vectorize.\n",
    "Let's look at [this](http://nbviewer.jupyter.org/github/jreback/pydata2015-london/blob/master/notebooks/idioms.ipynb)\n",
    "example, taken from Jeff Reback's PyData London talk, that groupwise normalizes a dataset by subtracting the mean and dividing by the standard deviation for each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def create_frame(n, n_groups):\n",
    "    # just setup code, not benchmarking this\n",
    "    stamps = pd.date_range('20010101', periods=n, freq='ms')\n",
    "    random.shuffle(stamps.values)    \n",
    "    return pd.DataFrame({'name': np.random.randint(0,n_groups,size=n),\n",
    "                         'stamp': stamps,\n",
    "                         'value': np.random.randint(0,n,size=n),\n",
    "                         'value2': np.random.randn(n)})\n",
    "\n",
    "\n",
    "df = create_frame(1000000,10000)\n",
    "\n",
    "def f_apply(df):\n",
    "    # Typical transform\n",
    "    return df.groupby('name').value2.apply(lambda x: (x-x.mean())/x.std())\n",
    "\n",
    "def f_unwrap(df):\n",
    "    # \"unwrapped\"\n",
    "    g = df.groupby('name').value2\n",
    "    v = df.value2\n",
    "    return (v-g.transform(np.mean))/g.transform(np.std)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Timing it we see that the \"unwrapped\" version, get's quite a bit better performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 4.24 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit f_apply(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 75.4 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit f_unwrap(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas GroupBy objects intercept calls for common functions like mean, sum, etc. and substitutes them with optimized Cython versions.\n",
    "So the unwrapped `.transform(np.mean)` and `.transform(np.std)` are fast, while the `x.mean` and `x.std` in the `.apply(lambda x: x - x.mean()/x.std())` aren't.\n",
    "\n",
    "`Groupby.apply` is always going to be around, beacuse it offers maximum flexibility. If you need to [fit a model on each group and create additional columns in the process](http://stackoverflow.com/q/35924126/1889400), it can handle that. It just might not be the fastest (which may be OK sometimes).\n",
    "\n",
    "This last example is admittedly niche.\n",
    "I'd like to think that there aren't too many places in pandas where the natural thing to do `.transform((x - x.mean()) / x.std())` is slower than the less obvious alternative.\n",
    "Ideally the user wouldn't have to know about GroupBy having special fast implementations of common methods.\n",
    "But that's where we are now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categoricals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Thanks to some great work by [Jan Schulz](https://twitter.com/janschulz), [Jeff Reback](https://twitter.com/janschulz), and others, pandas 0.15 gained a new [Categorical](http://pandas.pydata.org/pandas-docs/version/0.18.0/categorical.html) data type. Categoricals are nice for many reasons beyond just efficiency, but we'll focus on that here.\n",
    "\n",
    "Categoricals are an efficient way of representing data (typically strings) that have a low *cardinality*, i.e. relatively few distinct values relative to the size of the array. Internally, a Categorical stores the categories once, and an array of `codes`, which are just integers that indicate which category belongs there. Since it's cheaper to store a `code` than a `category`, we save on memory (shown next).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "800.00 KB\n"
     ]
    }
   ],
   "source": [
    "import string\n",
    "\n",
    "s = pd.Series(np.random.choice(list(string.ascii_letters), 100000))\n",
    "print('{:0.2f} KB'.format(s.memory_usage(index=False) / 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100.42 KB\n"
     ]
    }
   ],
   "source": [
    "c = s.astype('category')\n",
    "print('{:0.2f} KB'.format(c.memory_usage(index=False) / 1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beyond saving memory, having codes and a fixed set of categories offers up a bunch of algorithmic optimizations that pandas and others can take advantage of.\n",
    "\n",
    "[Matthew Rocklin](https://twitter.com/mrocklin) has a very nice [post](http://matthewrocklin.com/blog/work/2015/06/18/Categoricals) on using categoricals, and optimizing code in general."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Going Further\n",
    "\n",
    "The pandas documentation has a section on [enhancing performance](http://pandas.pydata.org/pandas-docs/version/0.18.0/enhancingperf.html), focusing on using Cython or `numba` to speed up a computation. I've focused more on the lower-hanging fruit of picking the right algorithm, vectorizing your code, and using pandas or numpy more effetively. There are further optimizations availble if these aren't enough."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This post was more about how to make effective use of numpy and pandas, than writing your own highly-optimized code.\n",
    "In my day-to-day work of data analysis it's not worth the time to write and compile a cython extension.\n",
    "I'd rather rely on pandas to be fast at what matters (label lookup on large arrays, factorizations for groupbys and merges, numerics).\n",
    "If you want to learn more about what pandas does to make things fast, checkout Jeff Tratner' talk from PyData Seattle [talk](http://www.jeffreytratner.com/slides/pandas-under-the-hood-pydata-seattle-2015.pdf) on pandas' internals.\n",
    "\n",
    "Next time we'll look at a differnt kind of optimization: using the Tidy Data principles to facilitate efficient data analysis.\n"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}