lectures/lecture1/Boxes-and-registers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 1: Boxes and Registers\n",
    "\n",
    "This is an [IJulia notebook](https://github.com/JuliaLang/IJulia.jl) (the [Jupyter](http://jupyter.org/)-based front-end for Julia) for [18.S096 at MIT in IAP 2017](https://math.mit.edu/classes/18.S096/iap17/), designed to accompany the first lecture.\n",
    "\n",
    "The basic goal of this lecture is to understand why some code (in some languages and/or styles of coding) is slow while other code is fast, based on whether it can be compiled to efficiently use the CPU registers and low-level arithmetic instructions, or whether it relies on \"boxed\" types that force \"dynamic\" computations.\n",
    "\n",
    "To illustrate this, we will implement a **sum** function `sum(a)`, which computes\n",
    "\n",
    "$$\n",
    "\\mathrm{sum}(a) = \\sum_{i=1}^n a_i\n",
    "$$\n",
    "\n",
    "for an array `a` with `n` elements.   We will use the built-in `sum` functions in Julia and Python along with hand-coded implementations in C, Python, and Julia.\n",
    "\n",
    "We will use some tricks so that we can write and benchmark C, Python, and Julia code all in the same notebook.  In the case of Python, this will rely on the [PyCall](https://github.com/JuliaPy/PyCall.jl) package to call Python from Julia.  We will use the [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl) Julia package to collect benchmarking statistics for us."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Low-level C code\n",
    "\n",
    "To start with, we will write a baseline implementation in the low-level C programming language.  Our C function `c_sum` will only work for a single data type: an array `X` of double-precision floating-point values (`double` in C, or `Float64` in Julia).\n",
    "\n",
    "(In contrast, our Julia code, and some of our Python code, will work for any numeric type; we'll see whether we pay a price for this.)\n",
    "\n",
    "Julia can easily call C functions in shared libraries via its `ccall` syntax.  So, we'll take our C routine (in a string) and pipe it through the C compiler `gcc` to produce a shared library file that we can load and call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "c_sum (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_code = \"\"\"\n",
    "#include <stddef.h>\n",
    "double c_sum(size_t n, double *X) {\n",
    "    double s = 0.0;\n",
    "    for (size_t i = 0; i < n; ++i) {\n",
    "        s += X[i];\n",
    "    }\n",
    "    return s;\n",
    "}\n",
    "\"\"\"\n",
    "# compile to a shared library by piping C_code to gcc:\n",
    "# (only works if you have gcc installed)\n",
    "const Clib = tempname()\n",
    "open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * \".\" * Libdl.dlext) -`, \"w\") do f\n",
    "    print(f, C_code)\n",
    "end\n",
    "c_sum(X::Array{Float64}) = ccall((\"c_sum\", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we should first check whether our function is correct, by comparing it to Julia's built-in `sum` function on an array of $10^7$ random numbers.  Different floating-point algorithms for the `sum` function give slightly different results (Julia's `sum` algorithm is actually *much more accurate* than the one here, but that's a story for another day), so we'll compute their \"fractional difference\" or \"relative error\" and make sure that this is small.  \n",
    "\n",
    "(Double-precision floating-point arithmetic keeps about 15 decimal digits, so any relative error close to $10^{-15}$ is a reasonable amount of roundoff error.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4676984614603977e-13"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define a function to compute the relative (fractional) error |x-y| / mean(|x|,|y|)\n",
    "relerr(x,y) = abs(x - y) * 2 / (abs(x) + abs(y))\n",
    "\n",
    "a = rand(10^7) # array of random numbers in [0,1)\n",
    "relerr(c_sum(a), sum(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collecting accurate benchmarking statistics can be a tricky business, so we'll use the Julia `BenchmarkTools` package to do most of the work.   If you don't have it installed, you may need to type `Pkg.add(\"BenchmarkTools\")` to tell Julia to download and install it.\n",
    "\n",
    "It defines a *macro* `@btime` that takes some Julia code and *transforms* it into a benchmark measuring the speed of that code.   We pass the argument `a` of the `c_sum` function to be benchmarked by the special syntax `$a` for technical reasons, basically to make sure that Julia's analysis of the variable `a` happens *before* the benchmark starts.  Macro syntax and **metaprogramming** will be a topic for another lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  9.860 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using BenchmarkTools, Compat\n",
    "\n",
    "c_bench = @btime c_sum($a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the  time is around **10 ms** for summing $10^7$ numbers, or about **1 billion additions per second**.  That sounds like a lot, but on my **2.5GHz laptop** it is well below the peak rate at which the computer can perform arithmetic.   It doesn't reach the peak arithmetic rate because for each floating-point addition, the processor needs to perform several additional calculations to load the next element of the array from memory, not to mention the time for the memory access itself.  Of course, you may get a slightly different number if you run this benchmark on a different computer.\n",
    "\n",
    "This **10 ms** number for type-specific compiled C code is a baseline against which we will compare our other implementations of summation, below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python sum functions\n",
    "\n",
    "Now, we'll call the Python functions and benchmark them.   The PyCall package allows to load Python as a library and to call it directly from Julia, sharing memory with Python and passing data and functions back and forth.  There is very little overhead to this, and in any case we will be summing $10^7$ numbers so the overhead of the Julia/Python interface is negligible compared to the cost of the summation itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "v\"2.7.14\""
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyCall\n",
    "PyCall.pyversion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### built-in `sum` of a Python `list`\n",
    "\n",
    "To start with, I will convert our array `a` into a Python `list` (the built-in Python array-like data structure), and sum it with the built-in Python `sum` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PyObject <built-in function sum>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# call a low-level PyCall function to get a Python list, because\n",
    "# by default PyCall will convert to a NumPy array instead (we benchmark NumPy below):\n",
    "apy_list = PyCall.array2py(a, 1, 1)\n",
    "# get the Python built-in \"sum\" function:\n",
    "pysum = pybuiltin(\"sum\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that we can call it, and that it computes the same answer as the Julia `sum`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4676984614603977e-13"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relerr(pysum(apy_list), sum(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll benchmark it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  101.506 ms (17 allocations: 512 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py_list_bench = @btime $pysum($apy_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It takes around **100 ms**, or almost **10x slower** than the C routine above.  This is true **even though the Python sum** function is [written in almost 200 lines of C code](https://github.com/python/cpython/blob/b5be30c92ff8cc0ac83f48014fe78bc048141021/Python/bltinmodule.c#L2191-L2365).   The problem is that the Python code *pays a price for being generic*: it handles arbitrary iterable data structures of arbitrary numeric \"boxes\" (`PyObject` pointers), and has to perform lots of computations both to fetch each item and also to perform each addition.\n",
    "\n",
    "### NumPy `sum` of a NumPy `array`\n",
    "\n",
    "You can do *much better* if you can take advantage of the fact that *all of the elements are the same type*.  Then, you can store the array as the actual floating-point data stored consecutively in memory (not an array of pointers to boxes), and your inner loop can be fast because the type checks can occur *outside* the loop. In Python, this kind of **homogeneous array** is exactly what is provided by [NumPy](http://www.numpy.org/). Internally, a NumPy array is essentially just a wrapper around a C-like `double*` array.   NumPy\n",
    "also provides `numpy.sum` function that can sum a NumPy array quickly.\n",
    "\n",
    "There is a catch, though: NumPy itself is written mostly in C, not Python.  And because C code is not type-generic, in order to handle a wide variety of NumPy array types (integers, double precision, single precision, etcetera), NumPy uses rather tricky **auto-generated C code**.  And even then it can only handle a small set of commonly used types; you can't define your own types and sum them quickly.\n",
    "\n",
    "Anyway, we can easily convert a Julia array to a NumPy array with PyCall (in fact, PyCall does this by\n",
    "by default), and benchmark `numpy.sum`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.341 ms (22 allocations: 720 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117e6"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy_sum = pyimport(\"numpy\")[\"sum\"]\n",
    "apy_numpy = PyObject(a) # converts to a numpy array by default\n",
    "py_numpy_bench = @btime $numpy_sum($apy_numpy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "WOW, it is actually **roughly twice as fast** as our C function!\n",
    "\n",
    "The reason for this extra boost is that the NumPy functions exploit [SIMD instructions](https://en.wikipedia.org/wiki/Streaming_SIMD_Extensions): special CPU instructions that can perform multiple additions at once, which we didn't use in our C code above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hand-written Python `sum` function\n",
    "\n",
    "To complete the story, let's write our own `mysum` function in Python that sums an arbitrary Python list (or array, or any iterable Python container).\n",
    "\n",
    "Of course, you would never do this for summation — you would always use one of the built-in\n",
    "functions in practice.  But someday, you will inevitably run into a problem where the\n",
    "performance-critical code has not already been written for you, and you will need to write\n",
    "your own.  So it is a good exercise to see how easy it is to get performance that\n",
    "is comparable to the library routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PyObject <function mysum at 0x14a7eccf8>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py\"\"\"\n",
    "def mysum(a):\n",
    "    s = 0.0\n",
    "    for x in a:\n",
    "        s = s + x\n",
    "    return s\n",
    "\"\"\"\n",
    "mysum_py = py\"mysum\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, let's check that it works first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4676984614603977e-13"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relerr(mysum_py(apy_list), sum(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's time it on our Python list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.450 s (17 allocations: 512 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime $mysum_py($apy_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yikes, about **1 s**.  That's more than **10× slower than the Python `sum`** and about **100× slower than our C code**.\n",
    "\n",
    "Using our `mysum` function with the NumPy array is no better, and in fact is a bit worse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.290 s (22 allocations: 720 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime $mysum_py($apy_numpy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can't take advantage of the NumPy array format in Python itself — you still have to write the performance-critical code in C (or hope someone else has written it for you)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Built-in Julia `sum` function\n",
    "\n",
    "Now, let's try the same thing in Julia, starting with the built-in Julia `sum` function operating on our array `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.442 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.0002245561170075e6"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j_bench = @btime sum($a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hooray, about **2× the speed of the C code**, comparable to `numpy.sum`.  Again, you can guess that it must be using SIMD to beat the C code.  And, again, it must be taking advantage of the fact that the array is homogeneous.\n",
    "\n",
    "The type of `a` is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array{Float64,1}"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the Julia type for a 1-dimensional array of `Float64` values (64-bit \"double\" precision floating-point numbers, equivalent to C `double`).  Because the type of the elements is \"attached\" to the type of the array (a \"parameterized\" type, more on this later), Julia is able to store it as a \"flat\" array of consecutive `Float64` values in memory.\n",
    "\n",
    "In contrast, the Julia equivalent of a Python `list` is a `Vector{Any}` (a synonym for `Array{Any,1}`): internally, this is an array of pointers to \"boxes\" that can hold any type (`Any`).  This makes things *much* slower: each `+` computation on an `Any` value must dynamically look up the type of object, figure out what `+` function to call, and allocate a new \"box\" to store the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  334.033 ms (10032767 allocations: 153.09 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.0002245561170075e6"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_any = Vector{Any}(a)\n",
    "j_bench_any = @btime sum($a_any)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is about a **50× slowdown**.  It is more than **2× slower than the Python `sum(list)` code**, in fact.  The Python `sum` function is better optimized than the Julia `sum` function for dealing with untyped (`Any`) values, in part because in Julia it is expected that you will use \"concretely\" typed arrays in all performance-critical cases.\n",
    "\n",
    "Unlike NumPy, however, Julia allows you to make efficient homogeneous arrays for any data type, even data types you define yourself, and you can operate on them efficiently with code written in Julia itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hand-written Julia `sum` functions\n",
    "\n",
    "Let's try to write our own `sum` function in Julia, just as we wrote our own Python function.  We'll implement four different versions and see how they compare.  We'll start simple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4676984614603977e-13"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum1(A)\n",
    "    s = 0\n",
    "    for a in A\n",
    "        s += a\n",
    "    end\n",
    "    return s\n",
    "end\n",
    "relerr(mysum1(a), sum(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  302.440 ms (29999999 allocations: 457.76 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j1_bench = @btime mysum1($a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yikes!  The performance is **terrible**, more than **200ms**.\n",
    "\n",
    "However, it turns out we made a simple mistake: we had a **type instability** in our code, because the summation variable `s` **changes type** in our function.  It starts out as an *integer* `0` (a Julia `Int`), but once we do `s += a` (shorthand for `s = s + a`), the result *changes* to a floating-point value.  Julia's compiler is quite simple-minded about such things: **once it sees that a variable changes type, it assumes it can be any type** and stores it in a \"box\".\n",
    "\n",
    "One symptom of this is the \"allocs estimate: 30000000\" line above: it is doing *zillions of allocations* in order to allocate all the \"boxes\" for the untyped `s` values.  Another way of seeing this would be to use the `@code_warntype` macro provided with Julia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variables:\n",
      "  #self#::#mysum1\n",
      "  A::Array{Float64,1}\n",
      "  a::Float64\n",
      "  #temp#@_4::Int64\n",
      "  s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n",
      "  #temp#@_6::Core.MethodInstance\n",
      "  #temp#@_7::Float64\n",
      "\n",
      "Body:\n",
      "  begin \n",
      "      s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m = 0 # line 3:\n",
      "      #temp#@_4::Int64 = 1\n",
      "      4: \n",
      "      unless (Base.not_int)((#temp#@_4::Int64 === (Base.add_int)((Base.arraylen)(A::Array{Float64,1})::Int64, 1)::Int64)::Bool)::Bool goto 29\n",
      "      SSAValue(2) = (Base.arrayref)(A::Array{Float64,1}, #temp#@_4::Int64)::Float64\n",
      "      SSAValue(3) = (Base.add_int)(#temp#@_4::Int64, 1)::Int64\n",
      "      a::Float64 = SSAValue(2)\n",
      "      #temp#@_4::Int64 = SSAValue(3) # line 4:\n",
      "      unless (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m isa Int64)::Bool goto 14\n",
      "      #temp#@_6::Core.MethodInstance = MethodInstance for +(::Int64, ::Float64)\n",
      "      goto 23\n",
      "      14: \n",
      "      unless (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m isa Float64)::Bool goto 18\n",
      "      #temp#@_6::Core.MethodInstance = MethodInstance for +(::Float64, ::Float64)\n",
      "      goto 23\n",
      "      18: \n",
      "      goto 20\n",
      "      20: \n",
      "      #temp#@_7::Float64 = (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m + a::Float64)::Float64\n",
      "      goto 25\n",
      "      23: \n",
      "      #temp#@_7::Float64 = $(Expr(:invoke, :(#temp#@_6), :(Main.+), :(s), :(a)))\n",
      "      25: \n",
      "      s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m = #temp#@_7::Float64\n",
      "      27: \n",
      "      goto 4\n",
      "      29:  # line 6:\n",
      "      return s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n",
      "  end\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n"
     ]
    }
   ],
   "source": [
    "@code_warntype mysum1(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tip-off here is the `s::Any` lines, telling you that `s` is stored in a \"box\" that holds type `Any`.  In computer-science lingo, we would say that **compiler type inference has failed** to determine the type of `s`.   More on this below.\n",
    "\n",
    "One solution would be to initialize `s = 0.0`, i.e. start `s` out as a floating-point value.  This would work for our floating-point array `a`, but then wouldn't work for other types of arrays.  Instead, Julia provides an `eltype(A)` function to fetch the *type of the elements* of `A`, and a `zero` function to *initialize `s` to the correct type of zero* for `A`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4676984614603977e-13"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum2(A)\n",
    "    s = zero(eltype(A))\n",
    "    for a in A\n",
    "        s += a\n",
    "    end\n",
    "    return s\n",
    "end\n",
    "relerr(mysum2(a), sum(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  10.139 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000224556117741e6"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j2_bench = @btime mysum2($a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that's more like it!  It runs in about **10 ms**, essentially the **same speed as the hand-written C code**.\n",
    "\n",
    "*Unlike* the C code, however, it works for *any* type of array, and in fact just about any type of \"iterable container\" (as long as it provides an `eltype` method), and can **sum any type of value** (as long as `zero` and `+` are defined), including user-defined types.  (We'll give an example below).\n",
    "\n",
    "The performance does not quite match the Julia built-in `sum` function or the `numpy.sum` function, however.  Our guess above was that they were exploiting SIMD optimizations.  However, we can do that too, in our own Julia code, by using the `@simd` decorator to tell Julia's compiler to turn on SIMD optimizations for that loop.\n",
    "\n",
    "(SIMD optimizations are not turned on by default, because they only speed up very particular kinds of code, and turning them on everywhere would slow down the compiler too much.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.4701475983568594e-15"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum3(A)\n",
    "    s = zero(eltype(A))\n",
    "    @simd for a in A\n",
    "        s += a\n",
    "    end\n",
    "    return s\n",
    "end\n",
    "relerr(mysum3(a), sum(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  6.108 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.00022455611703e6"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j3_bench = @btime mysum3($a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Hooray!  Basically the same speed as Julia's built-in `sum` function and `numpy.sum`!   And it only required **7 lines of code**, some care with types, and a very minor bit of wizardry with the `@simd` tag to get the last factor of two.\n",
    "\n",
    "Moreover, the code is still **type generic**: it can sum any container of any type that works with addition.   For example, it works for complex numbers, which are about two times slower as you might expect (since each complex addition requires two real additions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  14.389 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5.000159463459768e6 + 4.999881407607564e6im"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = rand(Complex{Float64}, length(a));\n",
    "@btime mysum3($z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we **didn't have to declare any types** of any arguments or variables; the compiler figured everything out.  How?"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}