02-numpy-and-pandas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br><br><br><br>\n",
    "\n",
    "# Numpy and Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with a plausible problem: analyze a dataset of daily Newark temperatures since 1883.\n",
    "\n",
    "import glob\n",
    "\n",
    "for filename in glob.glob(\"data/newark-temperature-*.txt\"):\n",
    "    print(\"-------------------------------\")\n",
    "    print(filename)\n",
    "    with open(filename) as file:\n",
    "        print(file.read()[:1000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the averages into arrays.\n",
    "\n",
    "tavg = []\n",
    "with open(\"data/newark-temperature-tavg.txt\") as file:\n",
    "    for line in file.readlines():\n",
    "        tavg.append(float(line))\n",
    "\n",
    "tmin = []\n",
    "with open(\"data/newark-temperature-tmin.txt\") as file:\n",
    "    for line in file.readlines():\n",
    "        tmin.append(float(line))\n",
    "\n",
    "tmax = []\n",
    "with open(\"data/newark-temperature-tmax.txt\") as file:\n",
    "    for line in file.readlines():\n",
    "        tmax.append(float(line))\n",
    "\n",
    "print(\"how many?\", len(tavg), len(tmin), len(tmax))    # having the same length is essential!\n",
    "print(\"starts with\", tavg[:3], tmin[:3], tmax[:3])\n",
    "print(\"ends with  \", tavg[-3:], tmin[-3:], tmax[-3:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The minima and maxima are more complete than the averages.\n",
    "\n",
    "import math\n",
    "\n",
    "print(\"fraction good in avrg:\", sum(0 if math.isnan(x) else 1 for x in tavg) / len(tavg))\n",
    "print(\"fraction good in mini:\", sum(0 if math.isnan(x) else 1 for x in tmin) / len(tmin))\n",
    "print(\"fraction good in maxi:\", sum(0 if math.isnan(x) else 1 for x in tmax) / len(tmax))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "\n",
    "# So let's \"impute\" averages: the measured average is ideal, but take (mini + maxi)/2 if unavailable.\n",
    "\n",
    "imputed = []\n",
    "for average, minimum, maximum in zip(tavg, tmin, tmax):\n",
    "    if math.isnan(average):\n",
    "        imputed.append(0.5*(minimum + maximum))\n",
    "    else:\n",
    "        imputed.append(average)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same thing in Numpy: load the data and impute missing averages.\n",
    "\n",
    "import numpy\n",
    "\n",
    "np_tavg = numpy.array(tavg)\n",
    "np_tmin = numpy.array(tmin)\n",
    "np_tmax = numpy.array(tmax)\n",
    "\n",
    "print(\"how many?\", len(np_tavg), len(np_tmin), len(np_tmax))\n",
    "print(\"starts with\", np_tavg[:3], np_tmin[:3], np_tmax[:3])\n",
    "print(\"ends with  \", np_tavg[-3:], np_tmin[-3:], np_tmax[-3:])\n",
    "\n",
    "print()\n",
    "print(\"fraction good in avrg:\", numpy.isnan(np_tavg).sum() / len(np_tavg))\n",
    "print(\"fraction good in mini:\", numpy.isnan(np_tmin).sum() / len(np_tmin))\n",
    "print(\"fraction good in maxi:\", numpy.isnan(np_tmax).sum() / len(np_tmax))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "\n",
    "#                        condition             if true                  if false\n",
    "np_imputed = numpy.where(numpy.isnan(np_tavg), 0.5*(np_tmin + np_tmax), np_tavg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br><br><br><br>\n",
    "\n",
    "Your milage may vary, but <tt>6.73 ms</tt> versus <tt>79.1 µs</tt> is a factor of 85!\n",
    "\n",
    "Factor 100‒1000 speedups from pure Python → Numpy are common, and they make the difference between 5 minutes (bathroom break) and 8 hours (overnight).\n",
    "\n",
    "<br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fundamentally different code order\n",
    "\n",
    "Also notice that we had to change the code from\n",
    "\n",
    "```python\n",
    "imputed = []\n",
    "for average, minimum, maximum in zip(tavg, tmin, tmax):\n",
    "    if math.isnan(average):\n",
    "        imputed.append(0.5*(minimum + maximum))\n",
    "    else:\n",
    "        imputed.append(average)\n",
    "```\n",
    "\n",
    "to\n",
    "\n",
    "```python\n",
    "np_imputed = numpy.where(numpy.isnan(np_tavg), 0.5*(np_tmin + np_tmax), np_tavg)\n",
    "```\n",
    "\n",
    "Pure Python is step-by-step, which can be good or bad. Numpy is all-at-once, which can be good or bad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step-by-step:**\n",
    "\n",
    "   * is **good** because you can insert breakpoints and watch variables to debug the code; it's like a microscopic view with no abstraction;\n",
    "   * is **bad** because the bigger picture can be lost when spread out among so many lines. (This is why I use list comprehensions.)\n",
    "\n",
    "**All-at-once:**\n",
    "\n",
    "   * is **good** because the composition of functions often reads like an English description of the problem to be solved;\n",
    "   * is **bad** because many indexes need to align; it's hard to break the process apart to debug it. (I usually get a line of Numpy right on the fifth try. Error messages are your friend.)\n",
    "\n",
    "**Trade-offs:**\n",
    "\n",
    "Pure Python is generally easier to _write,_ making it good for prototyping. Numpy is often easier to _read_.\n",
    "\n",
    "And, of course, Numpy is faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pure Python code order: acts on one DATUM at a time.\n",
    "# Numpy code order: acts on one ATTRIBUTE at a time.\n",
    "\n",
    "a = numpy.random.uniform(5, 10, 10000)\n",
    "b = numpy.random.uniform(10, 20, 10000)\n",
    "c = numpy.random.uniform(-0.1, 0.1, 10000)\n",
    "\n",
    "# Computes one quadratic formula on ai, bi, ci before moving on to the next one.\n",
    "roots1 = numpy.empty(10000, dtype=a.dtype)\n",
    "for i in range(10000):\n",
    "    roots1[i] = (-b[i] + math.sqrt(b[i]**2 - 4*a[i]*c[i])) / (2*a[i])\n",
    "\n",
    "# Computes one step in the quadratic formula for all 10000 before moving on to the next step.\n",
    "roots2 = (-b + numpy.sqrt(b**2 - 4*a*c)) / (2*a)\n",
    "\n",
    "print(roots1[:10])\n",
    "print(roots2[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Numpy expression (-b + numpy.sqrt(b**2 - 4*a*c)) / (2*a) actually computes something like:\n",
    "\n",
    "tmp1 = numpy.negative(b)            # -b\n",
    "tmp2 = numpy.square(b)              # b**2\n",
    "tmp3 = numpy.multiply(4, a)         # 4*a\n",
    "tmp4 = numpy.multiply(tmp3, c)      # tmp3*c\n",
    "tmp5 = numpy.subtract(tmp2, tmp4)   # tmp2 - tmp4\n",
    "tmp6 = numpy.sqrt(tmp5)             # sqrt(tmp5)\n",
    "tmp7 = numpy.add(tmp1, tmp6)        # tmp1 + tmp6\n",
    "tmp8 = numpy.multiply(2, a)         # 2*a\n",
    "roots3 = numpy.divide(tmp7, tmp8)   # tmp7 / tmp8\n",
    "\n",
    "print(roots1[:10])\n",
    "print(roots2[:10])\n",
    "print(roots3[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Not only are operations like + and numpy.sqrt array-at-a-time, but equality checks are, too.\n",
    "\n",
    "print(roots1 == roots2)\n",
    "print(roots2 == roots3)\n",
    "\n",
    "# To get a single answer to the question, \"Are these arrays equal?\" we have to reduce it with all().\n",
    "\n",
    "print()\n",
    "print((roots1 == roots2).all())\n",
    "print((roots2 == roots3).all())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Not all roots1 == roots2? Why not? (Note: on your platform, they might be identical...)\n",
    "\n",
    "# Which ones fail? (Indexes where roots1 != roots2 is nonzero/true.)\n",
    "failures, = numpy.nonzero(roots1 != roots2)\n",
    "print(failures[:20])\n",
    "\n",
    "# Show me the elements at the first index that fails!\n",
    "print()\n",
    "print(\"from roots1:\", roots1[failures[0]])\n",
    "print(\"from roots2:\", roots2[failures[0]])\n",
    "\n",
    "# How big is the difference?\n",
    "print()\n",
    "print(\"       diff:\", roots1[failures[0]] - roots2[failures[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Oh! Are they all these tiny differences?\n",
    "\n",
    "print(roots1[failures] - roots2[failures])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br><br>\n",
    "\n",
    "The difference came from `math.sqrt` versus `numpy.sqrt`: the latter is implemented differently on arrays. (Even using `numpy.sqrt` in both places doesn't help, because it runs different code on a single scalar than it does on arrays.)\n",
    "\n",
    "<br>\n",
    "\n",
    "**Moral:** do not expect pure Python math to agree perfectly with Numpy math!\n",
    "\n",
    "<br>\n",
    "\n",
    "(It wasn't the issue in this case, but in general, [floating point math is not associative](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)—the order matters at a microscopic scale.)\n",
    "\n",
    "<br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br><br><br><br>\n",
    "\n",
    "### Pandas: like Numpy, but convenient for analysis\n",
    "\n",
    "<br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Our last example was based on pre-formatted data. The data I downloaded from NOAA looked like this:\n",
    "\n",
    "with open(\"data/newark-temperature.csv\") as file:\n",
    "    print(file.read()[:10000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To use this with Numpy, we'd have to parse the CSV (import csv), handle the header line, etc.\n",
    "\n",
    "# Instead, we use Pandas.\n",
    "import pandas\n",
    "\n",
    "df = pandas.read_csv(\"data/newark-temperature.csv\", parse_dates=[\"DATE\"])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computations on Pandas columns work like Numpy (and mostly defer to Numpy for the actual work).\n",
    "\n",
    "(df[\"TMIN\"] + df[\"TMAX\"])/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# But there are many built-in conveniences. The fillna method replaces values only if they are NaN.\n",
    "\n",
    "df[\"imputed\"] = df[\"TAVG\"].fillna((df[\"TMIN\"] + df[\"TMAX\"]) / 2)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# It also has built-in plotting.\n",
    "\n",
    "df[\"imputed\"].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These things would be convenient by themselves, but there's also an essential difference:\n",
    "# Pandas data are INDEXED, meaning that there's a special column indicating what each row means.\n",
    "\n",
    "print(df.index)\n",
    "print()\n",
    "\n",
    "df2 = df.set_index(df[\"DATE\"])\n",
    "df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The DataFrame is now indexed by date, and so is each column (called a Series because it's more\n",
    "# than an array).\n",
    "\n",
    "df2[\"imputed\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Not only do plots now have a labeled x axis, but Pandas now knows about gaps in the data.\n",
    "\n",
    "df2[\"imputed\"].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arrays are implicitly \"indexed\" by integers, starting with zero or counting backward from -1.\n",
    "print(tmin[0])\n",
    "print(tmin[1])\n",
    "print(\"...\")\n",
    "print(tmin[-2])\n",
    "print(tmin[-1])\n",
    "# DataFrames and Series are indexed by their index type: timestamps in this case.\n",
    "print()\n",
    "print(df2[\"TMIN\"].loc[pandas.Timestamp(\"1893-01-01\")])\n",
    "print(df2[\"TMIN\"].loc[pandas.Timestamp(\"1893-01-02\")])\n",
    "print(\"...\")\n",
    "print(df2[\"TMIN\"].loc[pandas.Timestamp(\"2018-11-03\")])\n",
    "print(df2[\"TMIN\"].loc[pandas.Timestamp(\"2018-11-04\")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arrays are implicitly \"indexed\" by integers, starting with zero or counting backward from -1.\n",
    "print(tmin[0])\n",
    "print(tmin[1])\n",
    "print(\"...\")\n",
    "print(tmin[-2])\n",
    "print(tmin[-1])\n",
    "# There's also an iloc if you want to select by position, ignoring the index.\n",
    "print()\n",
    "print(df2[\"TMIN\"].iloc[0])\n",
    "print(df2[\"TMIN\"].iloc[1])\n",
    "print(\"...\")\n",
    "print(df2[\"TMIN\"].iloc[-2])\n",
    "print(df2[\"TMIN\"].iloc[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This means that you can select your data without converting data coordinates to indexes.\n",
    "\n",
    "df2[\"imputed\"].loc[pandas.Timestamp(\"2017-11-04\"):].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This allows for some concise expressions, keeping the meaning of the analysis logic clear.\n",
    "\n",
    "df2[[\"imputed\", \"TAVG\"]].loc[pandas.Timestamp(\"2011-01-01\"):pandas.Timestamp(\"2015-01-01\")].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Local warming! (Remember, this is only Newark.)\n",
    "\n",
    "df2[[\"imputed\", \"TAVG\"]].rolling(3*365).mean().plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Doing the same analysis in Numpy, without Pandas, is possible but more difficult.\n",
    "# (Note the lack of legend and proper x axis values/labels.)\n",
    "\n",
    "imputed = df2[\"imputed\"].values   # values gives us a Numpy array: no Pandas, no index\n",
    "\n",
    "import matplotlib.pyplot\n",
    "matplotlib.pyplot.plot(numpy.convolve(imputed, numpy.ones(3*365) / (3*365), mode=\"valid\"))\n",
    "matplotlib.pyplot.plot(numpy.convolve(tavg,    numpy.ones(3*365) / (3*365), mode=\"valid\"));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}