taxi-use-case.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Pilosa for transportation data\n",
    "This notebook demonstrates the use of Pilosa for analyzing taxi ride data, in the tradition of these guys:\n",
    "- http://tech.marksblogg.com/benchmarks.html\n",
    "- http://toddwschneider.com/posts/analyzing-1-1-billion-nyc-taxi-and-uber-trips-with-a-vengeance/\n",
    "\n",
    "To use this, you'll need an instance of Pilosa with an index populated according to the schema defined here: https://github.com/pilosa/pdk/blob/master/usecase/taxi/main.go\n",
    "\n",
    "* Install Pilosa: https://www.pilosa.com/docs/installation/\n",
    "* Getting started: https://www.pilosa.com/docs/getting-started/\n",
    "* To import data yourself, follow the instructions for the taxi usecase at https://github.com/pilosa/pdk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Dataset and Pilosa schema\n",
    "We're working with the csv files listed here: http://www.nyc.gov/html/tlc/html/about/trip_record_data.shtml. These have about 20 columns; about half of them are relevant for the benchmark queries we're looking at: cab type, distance, fare, number of passengers, dropoff and pickup time and location. We imported these fields, creating one or more Pilosa frames from each of them: \n",
    "\n",
    "| frame  | mapping |\n",
    "| ------ | ------- |\n",
    "|cab_type| direct map of enum int -> row ID |\n",
    "|dist_miles| round(dist) -> row ID |\n",
    "|total_amount_dollars| round(dist) -> row ID | \n",
    "|passenger_count| direct map of integer value -> row ID|\n",
    "|drop_grid_id| (lat, lon) -> 100x100 rectangular grid -> cell ID |\n",
    "|drop_year| year(timestamp) -> row ID |\n",
    "|drop_month| month(timestamp) -> row ID |\n",
    "|drop_day| day(timestamp) -> row ID |\n",
    "|drop_time| time of day mapped to one of 48 half-hour buckets |\n",
    "|pickup_grid_id| (lat, lon) -> 100x100 rectangular grid -> cell ID |\n",
    "|pickup_year| year(timestamp) -> row ID |\n",
    "|pickup_month| month(timestamp) -> row ID |\n",
    "|pickup_day| day(timestamp) -> row ID |\n",
    "|pickup_time| time of day mapped to one of 48 half-hour buckets |\n",
    "\n",
    "We also created two extra frames that represent the duration and average speed of each ride:\n",
    "\n",
    "| frame | mapping |\n",
    "| ----- | ------- |\n",
    "|duration_minutes| round(drop_timestamp - pickup_timestamp) -> row ID |\n",
    "|speed_mph| round(dist_miles / (drop_timestamp - pickup_timestamp)) -> row ID|\n",
    "\n",
    "This is just one example of a schema definition we can use to answer some interesting queries with Pilosa. In other cases, we might prefer different mappings to Pilosa bitmaps. Other examples include:\n",
    "\n",
    "* Mapping numeric values to custom bins instead of rounding (which produces regularly spaced bins).\n",
    "* Mapping timestamps to bins spanning the full range of times represented in the data, instead of using multiple frames for the timestamp components.\n",
    "* Using Pilosa's built-in support for timestamps.\n",
    "* Mapping (latitude, longitude) to real-world regions, such as NYC boroughs or neighborhoods.\n",
    "* Mapping string fields by matching or containment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import json\n",
    "import requests\n",
    "from itertools import product\n",
    "import pandas as pd\n",
    "import pilosa\n",
    "\n",
    "index = 'taxi'\n",
    "host = 'localhost:10101'\n",
    "\n",
    "qurl = 'http://%s/index/%s/query' % (host, index)\n",
    "client = pilosa.Client(host)\n",
    "pindex = pilosa.Index(index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "One common action is getting the total number of columns present in the index. We'll use this to demonstrate a few different ways of talking to the server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "224809\n",
      "224809\n"
     ]
    }
   ],
   "source": [
    "def get_ride_count():\n",
    "    # briefly demonstrate querying pilosa with \"raw http\"\n",
    "    resp = requests.post(qurl, data='TopN(frame=cab_type)')\n",
    "    data = json.loads(resp.content)['results'][0]\n",
    "    counts = [r['count'] for r in data]\n",
    "    return sum(counts)\n",
    "\n",
    "def get_ride_count_client():\n",
    "    # briefly demonstrate querying pilosa with the python client\n",
    "    pframe = pindex.frame('cab_type')\n",
    "    resp = client.query(pframe.topn(0))\n",
    "    data = resp.result.count_items\n",
    "    counts = [r.count for r in data]\n",
    "    return sum(counts)\n",
    "\n",
    "# try these out\n",
    "print(get_ride_count())\n",
    "print(get_ride_count_client())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now  let's try to query the index, using the four benchmarks described by Mark (http://tech.marksblogg.com/billion-nyc-taxi-rides-bigquery.html) for inspiration.\n",
    "\n",
    "---\n",
    "\n",
    "## Query 1: count per cab type\n",
    "\n",
    "This is a simple query, just a count of number of rides, by cab type. Pilosa's TopN query does exactly this, so it's fast and easy. `TopN(frame=foo)`, the simplest form of the TopN query, returns a list of row IDs in the frame \"foo\", in descending order of the *count of set bits*. \n",
    "\n",
    "This should produce a result like the following:\n",
    "```\n",
    "cab_type    count\n",
    "       0  1204143\n",
    "0.006066 sec\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cab_type   count\n",
      "       0  224809\n",
      "query: 0.007023 sec\n"
     ]
    }
   ],
   "source": [
    "# Build and post query\n",
    "t0 = time.time()\n",
    "q = 'TopN(frame=cab_type)'\n",
    "resp = requests.post(qurl, data=q)\n",
    "t1 = time.time()\n",
    "\n",
    "# Display results\n",
    "res = resp.json()['results'][0]\n",
    "df = pd.DataFrame({\n",
    "    'count': [d['count'] for d in res],\n",
    "    'cab_type': [d['id'] for d in res],\n",
    "})\n",
    "print(df.to_string(index=False))\n",
    "print('query: %f sec' % (t1-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "It took 11 lines of python to generate this result table, but the query is accomplished with only two or three lines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Query 2: average(total_amount) per passenger_count\n",
    "\n",
    "This query requires some postprocessing. We won't compute an average directly, rather we'll take advantage of the schema to get a few bitmap counts, then use those to compute an average on the client side. We'll use another form of the TopN query: `TopN(bar, frame=foo)`. This returns a list of rows from frame \"foo\", in descending order of the *count of set bits in the intersection with the bitmap \"bar\"*. \n",
    "\n",
    "We send one of these TopN queries for each row in the passenger_count frame, intersecting with the total_amount_dollars frame. We then compute one average(total_amount) from each of these queries. The denominator is the sum of `count` values. The numerator is the sum of `count` values, weighted by the corresponding total_amount value, which is approximately equal to the row ID (`key` in the result). This works because of our schema which maps fare values to row IDs by rounding to the nearest integer. \n",
    "\n",
    "Note: this calculation is approximate. However, it should produce a reasonable estimate, as long as the distribution of fare values is somewhat well behaved. Yellow taxi fares, for example, are not quite uniform, but they aren't skewed in a way that should badly skew this estimate.\n",
    "\n",
    "Pilosa also supports storage of scalar attributes per row or column; exact fare value can be stored this way. One feature on the Pilosa roadmap is to support aggregation over these attributes. When this is available, the averages in this query can be calculated internally.\n",
    "\n",
    "This should produce a result like the following:\n",
    "\n",
    "```average_amount  passenger_count\n",
    "      6.266355                0\n",
    "     14.318400                1\n",
    "     14.826454                2\n",
    "     14.830131                3\n",
    "     15.197937                4\n",
    "     14.529424                5\n",
    "     15.251607                6\n",
    "     17.000000                7\n",
    "      9.888889                8\n",
    "     19.470588                9\n",
    "query: 1.936771 sec\n",
    "postprocess: 0.002440 sec```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average_amount  passenger_count\n",
      "      6.666667                0\n",
      "     14.735619                1\n",
      "     15.044188                2\n",
      "     15.244204                3\n",
      "     16.557385                4\n",
      "     14.899677                5\n",
      "     15.535072                6\n",
      "      7.250000                7\n",
      "      9.000000                8\n",
      "     37.000000                9\n",
      "query: 0.364508 sec\n",
      "postprocess: 0.001544 sec\n"
     ]
    }
   ],
   "source": [
    "# Build and post query\n",
    "t0 = time.time()\n",
    "qs = ''\n",
    "pcounts = range(10)\n",
    "for i in pcounts:\n",
    "    qs += \"TopN(Bitmap(rowID=%d, frame='passenger_count'), frame=total_amount_dollars)\" % i\n",
    "resp = requests.post(qurl, data=qs)\n",
    "t1 = time.time()\n",
    "\n",
    "# Compute averages\n",
    "average_amounts = []\n",
    "for pcount, topn in zip(pcounts, resp.json()['results']):\n",
    "    # each `topn` corresponds to one passenger_count value, and contains \n",
    "    # a list of dicts, each containing a count and a key.\n",
    "    # The key is the bitmap id for the total_amount, and the count is the\n",
    "    # total number of rides matching that (passenger_count, total_amount)\n",
    "    # combination. Thanks to our schema, which maps these decimal fare \n",
    "    # values to integer-valued bitmap IDs by simple rounding, we can use \n",
    "    # the key itself as the value or weight in the average calculation. \n",
    "    wsum = sum([r['count'] * r['id'] for r in topn])\n",
    "    count = sum([r['count'] for r in topn])\n",
    "    if count == 0:\n",
    "        pcounts.remove(pcount)\n",
    "    else:\n",
    "        average_amounts.append(float(wsum)/count)\n",
    "t2 = time.time()\n",
    "\n",
    "# Display results\n",
    "df = pd.DataFrame({\n",
    "    'passenger_count': pcounts,\n",
    "    'average_amount': average_amounts,\n",
    "})\n",
    "df = df.reset_index(drop=True)\n",
    "print(df.to_string(index=False))\n",
    "print('query: %f sec' % (t1-t0))\n",
    "print('postprocess: %f sec' % (t2-t1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Query 3: count per (passenger_count, year)\n",
    "\n",
    "This one is relatively simple. The output should be a table with one row per (passenger_count, year) pair, and we need to send one `Count` query to produce the count value for each of those rows. Each (passenger_count, year) pair corresponds to an intersection of two rows, one from each frame.\n",
    "\n",
    "This should produce a result like the following:\n",
    "\n",
    "```Count  passenger_count  year\n",
    "953470                1  2013\n",
    " 92598                2  2013\n",
    " 27533                3  2013\n",
    " 11246                4  2013\n",
    "104491                5  2013\n",
    " 14312                6  2013\n",
    "    30                7  2013\n",
    "    18                8  2013\n",
    "    17                9  2013\n",
    "query: 0.153679 sec```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Count  passenger_count  year\n",
      "166309                1  2013\n",
      " 18557                2  2013\n",
      "  5823                3  2013\n",
      "  1943                4  2013\n",
      " 29395                5  2013\n",
      "  2723                6  2013\n",
      "     4                7  2013\n",
      "     1                8  2013\n",
      "     3                9  2013\n",
      "query: 0.045291 sec\n"
     ]
    }
   ],
   "source": [
    "# build and execute query\n",
    "t0 = time.time()\n",
    "qs = ''\n",
    "years = range(2009, 2016)\n",
    "pcounts = range(1, 10)\n",
    "\n",
    "for year, pcount in product(years, pcounts):\n",
    "    bmps = [\n",
    "        \"Bitmap(rowID=%d, frame='pickup_year')\" % year,\n",
    "        \"Bitmap(rowID=%d, frame='passenger_count')\" % pcount,\n",
    "    ]\n",
    "    qs += \"Count(Intersect(%s))\" % ', '.join(bmps)\n",
    "\n",
    "resp = requests.post(qurl, data=qs)\n",
    "t1 = time.time()\n",
    "\n",
    "# display results\n",
    "counts = resp.json()['results']\n",
    "df = pd.DataFrame({\n",
    "    'year': [x[0] for x in product(years, pcounts)],\n",
    "    'passenger_count': [x[1] for x in product(years, pcounts)],\n",
    "    'Count': counts\n",
    "})\n",
    "df = df[df.Count > 0]\n",
    "df = df.reset_index(drop=True)\n",
    "print(df.to_string(index=False))\n",
    "print('query: %f sec' % (t1-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Query 4: count per (passenger_count, year, round(trip_distance)) ordered by (year, count)\n",
    "\n",
    "To produce the full result table, we could iterate over the full product of all bitmaps in the three frames (passenger_count, year, distance), and send a `Count(Intersect())` query for each one - but this would be slow, and would only get worse if these frames were growing. A slightly smarter alternative is to make a quick guess at the size of these intersections, then query for the counts in that order. We then have the option to quit early, after a set number of result rows have been produced, or a set percentage of index rows have been accounted for (this latter approach is used here, controlled with the `pct_thresh` parameter).\n",
    "\n",
    "We can get the \"quick guess\" by sending TopN queries for each of the three frames, then combining the results, estimating the maximum possible size of the intersection as the minimum count of the three rows. We then sort this list on the size guess, and check the actual size in descending order. \n",
    "\n",
    "Generating the full result in this way is somewhat slow, but getting the top few rows is quick. The implementation here also sends the second group of requests serially, for simplicity - batching or parallelizing these would improve efficiency. That's exactly what we've done in our taxi demo (https://github.com/pilosa/demo-taxi/), which is a slightly more presentable version of these queries.\n",
    "\n",
    "This should produce a result like the following:\n",
    "\n",
    "```1204143\n",
    "Count  distance  passenger_count  year\n",
    "303011         1                1  2013\n",
    "202342         2                1  2013\n",
    "119368         3                1  2013\n",
    " 77763         4                1  2013\n",
    " 62200         0                1  2013\n",
    "   ...\n",
    "  1337         4                6  2013\n",
    "  1201         6                3  2013\n",
    "   944         7                3  2013\n",
    "   640         8                3  2013\n",
    "   427         9                3  2013\n",
    "fast query: 0.005938 sec\n",
    "slow query (43): 1.549860 sec\n",
    "postprocess: 0.004469 sec```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Count  distance  passenger_count  year\n",
      "47971         1                1  2013\n",
      "34382         2                1  2013\n",
      "20818         3                1  2013\n",
      "14176         0                1  2013\n",
      "13775         4                1  2013\n",
      " 9411         5                1  2013\n",
      " 8976         1                5  2013\n",
      " 6922         6                1  2013\n",
      " 6230         2                5  2013\n",
      " 5470         7                1  2013\n",
      " 4900         1                2  2013\n",
      " 3966         8                1  2013\n",
      " 3752         3                5  2013\n",
      " 3675         2                2  2013\n",
      " 2650         9                1  2013\n",
      " 2363         4                5  2013\n",
      " 2307         3                2  2013\n",
      " 2014         0                2  2013\n",
      " 1727         0                5  2013\n",
      " 1667         5                5  2013\n",
      " 1581         4                2  2013\n",
      " 1539         1                3  2013\n",
      " 1290         6                5  2013\n",
      " 1278         2                3  2013\n",
      " 1034         5                2  2013\n",
      " 1006         7                5  2013\n",
      "  841         6                2  2013\n",
      "  796         1                6  2013\n",
      "  751         3                3  2013\n",
      "  716         8                5  2013\n",
      "  569         7                2  2013\n",
      "  525         0                3  2013\n",
      "  520         2                6  2013\n",
      "  490         9                5  2013\n",
      "  455         4                3  2013\n",
      "  402         8                2  2013\n",
      "  376         3                6  2013\n",
      "  330         5                3  2013\n",
      "  312         9                2  2013\n",
      "  276         6                3  2013\n",
      "  245         4                6  2013\n",
      "  197         7                3  2013\n",
      "  187         5                6  2013\n",
      "  183         0                6  2013\n",
      "  149         8                3  2013\n",
      "  110         6                6  2013\n",
      "   84         9                3  2013\n",
      "   75         7                6  2013\n",
      "   69         8                6  2013\n",
      "   50         9                6  2013\n",
      "fast query: 0.005385 sec\n",
      "slow query (49): 0.478308 sec\n",
      "postprocess: 0.005201 sec\n"
     ]
    }
   ],
   "source": [
    "verbose = False  # set to True to print slow query updates\n",
    "pct_thresh = 95.0\n",
    "\n",
    "# build and post preliminary queries\n",
    "t0 = time.time()\n",
    "years = range(2009, 2016)\n",
    "pcounts = range(1, 9)\n",
    "dists = range(50)\n",
    "topns = [\n",
    "    \"TopN(frame='pickup_year')\"\n",
    "    \"TopN(frame='passenger_count')\"\n",
    "    \"TopN(frame='dist_miles')\"\n",
    "]\n",
    "qs = ', '.join(topns)\n",
    "resp = requests.post(qurl, data=qs)\n",
    "res = resp.json()['results']\n",
    "t1 = time.time()\n",
    "\n",
    "# assemble TopN results into candidates\n",
    "year_data = [(x['id'], x['count']) for x in res[0]]\n",
    "pcount_data = [(x['id'], x['count']) for x in res[1]]\n",
    "dist_data = [(x['id'], x['count']) for x in res[2]]\n",
    "\n",
    "cands = []\n",
    "for (k_year, c_year), (k_pcount, c_pcount), (k_dist, c_dist) in product(year_data, pcount_data, dist_data):\n",
    "    # this ugly loop creates a list of candidate tuples of the form\n",
    "    # (year, passenger_count, round(distance), max_intersection_count)\n",
    "    cands.append([k_year, k_pcount, k_dist, min([c_year, c_pcount, c_dist])])\n",
    "\n",
    "cands.sort(key=lambda x: -x[3])  # sort candidates by max_intersection_count\n",
    "\n",
    "# iterate over candidates in order of estimated largest\n",
    "n = 0\n",
    "total = 0\n",
    "num_rides = get_ride_count()\n",
    "years, pcounts, dists, counts = [], [], [], []\n",
    "for year, pcount, dist, maxcount in cands:\n",
    "    bmps = [\n",
    "        \"Bitmap(rowID=%d, frame='pickup_year')\" % year,\n",
    "        \"Bitmap(rowID=%d, frame='passenger_count')\" % pcount,\n",
    "        \"Bitmap(rowID=%d, frame='dist_miles')\" % dist,\n",
    "    ]\n",
    "    q = \"Count(Intersect(%s))\" % ', '.join(bmps)\n",
    "    resp = requests.post(qurl, data=q)\n",
    "    count = json.loads(resp.content)['results'][0]\n",
    "    total += count\n",
    "\n",
    "    years.append(year)\n",
    "    pcounts.append(pcount)\n",
    "    dists.append(dist)\n",
    "    counts.append(count)\n",
    "\n",
    "    pct = (100.*total)/num_rides\n",
    "    if verbose:\n",
    "        print('%d. %.2f%% (max %.2f%%)' % (n, pct, pct_thresh))\n",
    "    if pct >= pct_thresh:\n",
    "        break\n",
    "    n += 1\n",
    "t2 = time.time()\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    'year': years,\n",
    "    'passenger_count': pcounts,\n",
    "    'distance': dists,\n",
    "    'Count': counts,\n",
    "})\n",
    "df = df[df.Count > 0]\n",
    "df = df.sort_values(by=['year', 'Count'], ascending=[0, 0])\n",
    "t3 = time.time()\n",
    "\n",
    "print(df.to_string(index=False))\n",
    "print('fast query: %f sec' % (t1-t0))\n",
    "print('slow query (%d): %f sec' % (n, t2-t1))\n",
    "print('postprocess: %f sec' % (t3-t2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}