notebooks/02.03-Computation-on-arrays-ufuncs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PDSH-cover-small.png\">\n",
    "\n",
    "*This notebook contains an excerpt from the [Python Data Science Handbook](http://shop.oreilly.com/product/0636920034919.do) by Jake VanderPlas; the content is available [on GitHub](https://github.com/jakevdp/PythonDataScienceHandbook).*\n",
    "\n",
    "*The text is released under the [CC-BY-NC-ND license](https://creativecommons.org/licenses/by-nc-nd/3.0/us/legalcode), and code is released under the [MIT license](https://opensource.org/licenses/MIT). If you find this content useful, please consider supporting the work by [buying the book](http://shop.oreilly.com/product/0636920034919.do)!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [The Basics of NumPy Arrays](02.02-The-Basics-Of-NumPy-Arrays.ipynb) | [Contents](Index.ipynb) | [Aggregations: Min, Max, and Everything In Between](02.04-Computation-on-arrays-aggregates.ipynb) >\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/02.03-Computation-on-arrays-ufuncs.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computation on NumPy Arrays: Universal Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up until now, we have been discussing some of the basic nuts and bolts of NumPy; in the next few sections, we will dive into the reasons that NumPy is so important in the Python data science world.\n",
    "Namely, it provides an easy and flexible interface to optimized computation with arrays of data.\n",
    "\n",
    "Computation on NumPy arrays can be very fast, or it can be very slow.\n",
    "The key to making it fast is to use *vectorized* operations, generally implemented through NumPy's *universal functions* (ufuncs).\n",
    "This section motivates the need for NumPy's ufuncs, which can be used to make repeated calculations on array elements much more efficient.\n",
    "It then introduces many of the most common and useful arithmetic ufuncs available in the NumPy package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Slowness of Loops\n",
    "\n",
    "Python's default implementation (known as CPython) does some operations very slowly.\n",
    "This is in part due to the dynamic, interpreted nature of the language: the fact that types are flexible, so that sequences of operations cannot be compiled down to efficient machine code as in languages like C and Fortran.\n",
    "Recently there have been various attempts to address this weakness: well-known examples are the [PyPy](http://pypy.org/) project, a just-in-time compiled implementation of Python; the [Cython](http://cython.org) project, which converts Python code to compilable C code; and the [Numba](http://numba.pydata.org/) project, which converts snippets of Python code to fast LLVM bytecode.\n",
    "Each of these has its strengths and weaknesses, but it is safe to say that none of the three approaches has yet surpassed the reach and popularity of the standard CPython engine.\n",
    "\n",
    "The relative sluggishness of Python generally manifests itself in situations where many small operations are being repeated – for instance looping over arrays to operate on each element.\n",
    "For example, imagine we have an array of values and we'd like to compute the reciprocal of each.\n",
    "A straightforward approach might look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.16666667,  1.        ,  0.25      ,  0.25      ,  0.125     ])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "\n",
    "def compute_reciprocals(values):\n",
    "    output = np.empty(len(values))\n",
    "    for i in range(len(values)):\n",
    "        output[i] = 1.0 / values[i]\n",
    "    return output\n",
    "        \n",
    "values = np.random.randint(1, 10, size=5)\n",
    "compute_reciprocals(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This implementation probably feels fairly natural to someone from, say, a C or Java background.\n",
    "But if we measure the execution time of this code for a large input, we see that this operation is very slow, perhaps surprisingly so!\n",
    "We'll benchmark this with IPython's ``%timeit`` magic (discussed in [Profiling and Timing Code](01.07-Timing-and-Profiling.ipynb)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 2.91 s per loop\n"
     ]
    }
   ],
   "source": [
    "big_array = np.random.randint(1, 100, size=1000000)\n",
    "%timeit compute_reciprocals(big_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It takes several seconds to compute these million operations and to store the result!\n",
    "When even cell phones have processing speeds measured in Giga-FLOPS (i.e., billions of numerical operations per second), this seems almost absurdly slow.\n",
    "It turns out that the bottleneck here is not the operations themselves, but the type-checking and function dispatches that CPython must do at each cycle of the loop.\n",
    "Each time the reciprocal is computed, Python first examines the object's type and does a dynamic lookup of the correct function to use for that type.\n",
    "If we were working in compiled code instead, this type specification would be known before the code executes and the result could be computed much more efficiently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducing UFuncs\n",
    "\n",
    "For many types of operations, NumPy provides a convenient interface into just this kind of statically typed, compiled routine. This is known as a *vectorized* operation.\n",
    "This can be accomplished by simply performing an operation on the array, which will then be applied to each element.\n",
    "This vectorized approach is designed to push the loop into the compiled layer that underlies NumPy, leading to much faster execution.\n",
    "\n",
    "Compare the results of the following two:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.16666667  1.          0.25        0.25        0.125     ]\n",
      "[ 0.16666667  1.          0.25        0.25        0.125     ]\n"
     ]
    }
   ],
   "source": [
    "print(compute_reciprocals(values))\n",
    "print(1.0 / values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the execution time for our big array, we see that it completes orders of magnitude faster than the Python loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 4.6 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit (1.0 / big_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vectorized operations in NumPy are implemented via *ufuncs*, whose main purpose is to quickly execute repeated operations on values in NumPy arrays.\n",
    "Ufuncs are extremely flexible – before we saw an operation between a scalar and an array, but we can also operate between two arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.        ,  0.5       ,  0.66666667,  0.75      ,  0.8       ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(5) / np.arange(1, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And ufunc operations are not limited to one-dimensional arrays–they can also act on multi-dimensional arrays as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1,   2,   4],\n",
       "       [  8,  16,  32],\n",
       "       [ 64, 128, 256]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(9).reshape((3, 3))\n",
    "2 ** x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computations using vectorization through ufuncs are nearly always more efficient than their counterpart implemented using Python loops, especially as the arrays grow in size.\n",
    "Any time you see such a loop in a Python script, you should consider whether it can be replaced with a vectorized expression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring NumPy's UFuncs\n",
    "\n",
    "Ufuncs exist in two flavors: *unary ufuncs*, which operate on a single input, and *binary ufuncs*, which operate on two inputs.\n",
    "We'll see examples of both these types of functions here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array arithmetic\n",
    "\n",
    "NumPy's ufuncs feel very natural to use because they make use of Python's native arithmetic operators.\n",
    "The standard addition, subtraction, multiplication, and division can all be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x     = [0 1 2 3]\n",
      "x + 5 = [5 6 7 8]\n",
      "x - 5 = [-5 -4 -3 -2]\n",
      "x * 2 = [0 2 4 6]\n",
      "x / 2 = [ 0.   0.5  1.   1.5]\n",
      "x // 2 = [0 0 1 1]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(4)\n",
    "print(\"x     =\", x)\n",
    "print(\"x + 5 =\", x + 5)\n",
    "print(\"x - 5 =\", x - 5)\n",
    "print(\"x * 2 =\", x * 2)\n",
    "print(\"x / 2 =\", x / 2)\n",
    "print(\"x // 2 =\", x // 2)  # floor division"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also a unary ufunc for negation, and a ``**`` operator for exponentiation, and a ``%`` operator for modulus:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-x     =  [ 0 -1 -2 -3]\n",
      "x ** 2 =  [0 1 4 9]\n",
      "x % 2  =  [0 1 0 1]\n"
     ]
    }
   ],
   "source": [
    "print(\"-x     = \", -x)\n",
    "print(\"x ** 2 = \", x ** 2)\n",
    "print(\"x % 2  = \", x % 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, these can be strung together however you wish, and the standard order of operations is respected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.  , -2.25, -4.  , -6.25])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-(0.5*x + 1) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of these arithmetic operations are simply convenient wrappers around specific functions built into NumPy; for example, the ``+`` operator is a wrapper for the ``add`` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 3, 4, 5])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.add(x, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following table lists the arithmetic operators implemented in NumPy:\n",
    "\n",
    "| Operator\t    | Equivalent ufunc    | Description                           |\n",
    "|---------------|---------------------|---------------------------------------|\n",
    "|``+``          |``np.add``           |Addition (e.g., ``1 + 1 = 2``)         |\n",
    "|``-``          |``np.subtract``      |Subtraction (e.g., ``3 - 2 = 1``)      |\n",
    "|``-``          |``np.negative``      |Unary negation (e.g., ``-2``)          |\n",
    "|``*``          |``np.multiply``      |Multiplication (e.g., ``2 * 3 = 6``)   |\n",
    "|``/``          |``np.divide``        |Division (e.g., ``3 / 2 = 1.5``)       |\n",
    "|``//``         |``np.floor_divide``  |Floor division (e.g., ``3 // 2 = 1``)  |\n",
    "|``**``         |``np.power``         |Exponentiation (e.g., ``2 ** 3 = 8``)  |\n",
    "|``%``          |``np.mod``           |Modulus/remainder (e.g., ``9 % 4 = 1``)|\n",
    "\n",
    "Additionally there are Boolean/bitwise operators; we will explore these in [Comparisons, Masks, and Boolean Logic](02.06-Boolean-Arrays-and-Masks.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Absolute value\n",
    "\n",
    "Just as NumPy understands Python's built-in arithmetic operators, it also understands Python's built-in absolute value function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 1, 0, 1, 2])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([-2, -1, 0, 1, 2])\n",
    "abs(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding NumPy ufunc is ``np.absolute``, which is also available under the alias ``np.abs``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 1, 0, 1, 2])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.absolute(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 1, 0, 1, 2])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ufunc can also handle complex data, in which the absolute value returns the magnitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.,  5.,  2.,  1.])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([3 - 4j, 4 - 3j, 2 + 0j, 0 + 1j])\n",
    "np.abs(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trigonometric functions\n",
    "\n",
    "NumPy provides a large number of useful ufuncs, and some of the most useful for the data scientist are the trigonometric functions.\n",
    "We'll start by defining an array of angles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "theta = np.linspace(0, np.pi, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute some trigonometric functions on these values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta      =  [ 0.          1.57079633  3.14159265]\n",
      "sin(theta) =  [  0.00000000e+00   1.00000000e+00   1.22464680e-16]\n",
      "cos(theta) =  [  1.00000000e+00   6.12323400e-17  -1.00000000e+00]\n",
      "tan(theta) =  [  0.00000000e+00   1.63312394e+16  -1.22464680e-16]\n"
     ]
    }
   ],
   "source": [
    "print(\"theta      = \", theta)\n",
    "print(\"sin(theta) = \", np.sin(theta))\n",
    "print(\"cos(theta) = \", np.cos(theta))\n",
    "print(\"tan(theta) = \", np.tan(theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values are computed to within machine precision, which is why values that should be zero do not always hit exactly zero.\n",
    "Inverse trigonometric functions are also available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x         =  [-1, 0, 1]\n",
      "arcsin(x) =  [-1.57079633  0.          1.57079633]\n",
      "arccos(x) =  [ 3.14159265  1.57079633  0.        ]\n",
      "arctan(x) =  [-0.78539816  0.          0.78539816]\n"
     ]
    }
   ],
   "source": [
    "x = [-1, 0, 1]\n",
    "print(\"x         = \", x)\n",
    "print(\"arcsin(x) = \", np.arcsin(x))\n",
    "print(\"arccos(x) = \", np.arccos(x))\n",
    "print(\"arctan(x) = \", np.arctan(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponents and logarithms\n",
    "\n",
    "Another common type of operation available in a NumPy ufunc are the exponentials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x     = [1, 2, 3]\n",
      "e^x   = [  2.71828183   7.3890561   20.08553692]\n",
      "2^x   = [ 2.  4.  8.]\n",
      "3^x   = [ 3  9 27]\n"
     ]
    }
   ],
   "source": [
    "x = [1, 2, 3]\n",
    "print(\"x     =\", x)\n",
    "print(\"e^x   =\", np.exp(x))\n",
    "print(\"2^x   =\", np.exp2(x))\n",
    "print(\"3^x   =\", np.power(3, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse of the exponentials, the logarithms, are also available.\n",
    "The basic ``np.log`` gives the natural logarithm; if you prefer to compute the base-2 logarithm or the base-10 logarithm, these are available as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x        = [1, 2, 4, 10]\n",
      "ln(x)    = [ 0.          0.69314718  1.38629436  2.30258509]\n",
      "log2(x)  = [ 0.          1.          2.          3.32192809]\n",
      "log10(x) = [ 0.          0.30103     0.60205999  1.        ]\n"
     ]
    }
   ],
   "source": [
    "x = [1, 2, 4, 10]\n",
    "print(\"x        =\", x)\n",
    "print(\"ln(x)    =\", np.log(x))\n",
    "print(\"log2(x)  =\", np.log2(x))\n",
    "print(\"log10(x) =\", np.log10(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also some specialized versions that are useful for maintaining precision with very small input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exp(x) - 1 = [ 0.          0.0010005   0.01005017  0.10517092]\n",
      "log(1 + x) = [ 0.          0.0009995   0.00995033  0.09531018]\n"
     ]
    }
   ],
   "source": [
    "x = [0, 0.001, 0.01, 0.1]\n",
    "print(\"exp(x) - 1 =\", np.expm1(x))\n",
    "print(\"log(1 + x) =\", np.log1p(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When ``x`` is very small, these functions give more precise values than if the raw ``np.log`` or ``np.exp`` were to be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specialized ufuncs\n",
    "\n",
    "NumPy has many more ufuncs available, including hyperbolic trig functions, bitwise arithmetic, comparison operators, conversions from radians to degrees, rounding and remainders, and much more.\n",
    "A look through the NumPy documentation reveals a lot of interesting functionality.\n",
    "\n",
    "Another excellent source for more specialized and obscure ufuncs is the submodule ``scipy.special``.\n",
    "If you want to compute some obscure mathematical function on your data, chances are it is implemented in ``scipy.special``.\n",
    "There are far too many functions to list them all, but the following snippet shows a couple that might come up in a statistics context:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import special"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma(x)     = [  1.00000000e+00   2.40000000e+01   3.62880000e+05]\n",
      "ln|gamma(x)| = [  0.           3.17805383  12.80182748]\n",
      "beta(x, 2)   = [ 0.5         0.03333333  0.00909091]\n"
     ]
    }
   ],
   "source": [
    "# Gamma functions (generalized factorials) and related functions\n",
    "x = [1, 5, 10]\n",
    "print(\"gamma(x)     =\", special.gamma(x))\n",
    "print(\"ln|gamma(x)| =\", special.gammaln(x))\n",
    "print(\"beta(x, 2)   =\", special.beta(x, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "erf(x)  = [ 0.          0.32862676  0.67780119  0.84270079]\n",
      "erfc(x) = [ 1.          0.67137324  0.32219881  0.15729921]\n",
      "erfinv(x) = [ 0.          0.27246271  0.73286908         inf]\n"
     ]
    }
   ],
   "source": [
    "# Error function (integral of Gaussian)\n",
    "# its complement, and its inverse\n",
    "x = np.array([0, 0.3, 0.7, 1.0])\n",
    "print(\"erf(x)  =\", special.erf(x))\n",
    "print(\"erfc(x) =\", special.erfc(x))\n",
    "print(\"erfinv(x) =\", special.erfinv(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many, many more ufuncs available in both NumPy and ``scipy.special``.\n",
    "Because the documentation of these packages is available online, a web search along the lines of \"gamma function python\" will generally find the relevant information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced Ufunc Features\n",
    "\n",
    "Many NumPy users make use of ufuncs without ever learning their full set of features.\n",
    "We'll outline a few specialized features of ufuncs here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying output\n",
    "\n",
    "For large calculations, it is sometimes useful to be able to specify the array where the result of the calculation will be stored.\n",
    "Rather than creating a temporary array, this can be used to write computation results directly to the memory location where you'd like them to be.\n",
    "For all ufuncs, this can be done using the ``out`` argument of the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.  10.  20.  30.  40.]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(5)\n",
    "y = np.empty(5)\n",
    "np.multiply(x, 10, out=y)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can even be used with array views. For example, we can write the results of a computation to every other element of a specified array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  1.   0.   2.   0.   4.   0.   8.   0.  16.   0.]\n"
     ]
    }
   ],
   "source": [
    "y = np.zeros(10)\n",
    "np.power(2, x, out=y[::2])\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we had instead written ``y[::2] = 2 ** x``, this would have resulted in the creation of a temporary array to hold the results of ``2 ** x``, followed by a second operation copying those values into the ``y`` array.\n",
    "This doesn't make much of a difference for such a small computation, but for very large arrays the memory savings from careful use of the ``out`` argument can be significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregates\n",
    "\n",
    "For binary ufuncs, there are some interesting aggregates that can be computed directly from the object.\n",
    "For example, if we'd like to *reduce* an array with a particular operation, we can use the ``reduce`` method of any ufunc.\n",
    "A reduce repeatedly applies a given operation to the elements of an array until only a single result remains.\n",
    "\n",
    "For example, calling ``reduce`` on the ``add`` ufunc returns the sum of all elements in the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(1, 6)\n",
    "np.add.reduce(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, calling ``reduce`` on the ``multiply`` ufunc results in the product of all array elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.multiply.reduce(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we'd like to store all the intermediate results of the computation, we can instead use ``accumulate``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  3,  6, 10, 15])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.add.accumulate(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1,   2,   6,  24, 120])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.multiply.accumulate(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for these particular cases, there are dedicated NumPy functions to compute the results (``np.sum``, ``np.prod``, ``np.cumsum``, ``np.cumprod``), which we'll explore in [Aggregations: Min, Max, and Everything In Between](02.04-Computation-on-arrays-aggregates.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outer products\n",
    "\n",
    "Finally, any ufunc can compute the output of all pairs of two different inputs using the ``outer`` method.\n",
    "This allows you, in one line, to do things like create a multiplication table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2,  3,  4,  5],\n",
       "       [ 2,  4,  6,  8, 10],\n",
       "       [ 3,  6,  9, 12, 15],\n",
       "       [ 4,  8, 12, 16, 20],\n",
       "       [ 5, 10, 15, 20, 25]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.arange(1, 6)\n",
    "np.multiply.outer(x, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``ufunc.at`` and ``ufunc.reduceat`` methods, which we'll explore in [Fancy Indexing](02.07-Fancy-Indexing.ipynb), are very helpful as well.\n",
    "\n",
    "Another extremely useful feature of ufuncs is the ability to operate between arrays of different sizes and shapes, a set of operations known as *broadcasting*.\n",
    "This subject is important enough that we will devote a whole section to it (see [Computation on Arrays: Broadcasting](02.05-Computation-on-arrays-broadcasting.ipynb))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ufuncs: Learning More"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More information on universal functions (including the full list of available functions) can be found on the [NumPy](http://www.numpy.org) and [SciPy](http://www.scipy.org) documentation websites.\n",
    "\n",
    "Recall that you can also access information directly from within IPython by importing the packages and using IPython's tab-completion and help (``?``) functionality, as described in [Help and Documentation in IPython](01.01-Help-And-Documentation.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [The Basics of NumPy Arrays](02.02-The-Basics-Of-NumPy-Arrays.ipynb) | [Contents](Index.ipynb) | [Aggregations: Min, Max, and Everything In Between](02.04-Computation-on-arrays-aggregates.ipynb) >\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/jakevdp/PythonDataScienceHandbook/blob/master/notebooks/02.03-Computation-on-arrays-ufuncs.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}