lectures/lecture3/Type-stability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 2: Type Stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defeating type inference: Type instabilities\n",
    "\n",
    "To get good performance, there are some fairly simple rules that you need to follow in Julia code to avoid defeating the compiler's type inference.   See also the [performance tips section of the Julia manual](http://docs.julialang.org/en/stable/manual/performance-tips/).\n",
    "\n",
    "Three of the most important are:\n",
    "\n",
    "* Don't use (non-constant) global variables in critical code — put your critical code into a function (this is good advice anyway, from a software-engineering standpoint).  The compiler assumes that a **global variable can change type at any time**, so it is always stored in a \"box\", and \"taints\" anything that depends on it.\n",
    "\n",
    "* Local variables should be \"type-stable\": **don't change the type of a variable inside a function**.  Use a new variable instead.\n",
    "\n",
    "* Functions should be \"type-stable\": **a function's return type should only depend on the argument types, not on the argument values**.\n",
    "\n",
    "To diagnose all of these problems, the `@code_warntype` macro that we used above is your friend.  If it labels any variables (or the function's return value) as `Any` or `Union{...}`, it means that the compiler couldn't figure out a precise type.\n",
    "\n",
    "The third point, type-stability of functions, leads to lots of important but subtle choices in library design.  For example, consider the (built-in) `sqrt(x)` function, which computes $\\sqrt{x}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqrt(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might think that `sqrt(-1)` should return $i$ (or `im`, in Julia syntax).  (Matlab's `sqrt` function does this.)  Instead, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "DomainError:\nsqrt will only return a complex result if called with a complex argument. Try sqrt(complex(x)).",
     "output_type": "error",
     "traceback": [
      "DomainError:\nsqrt will only return a complex result if called with a complex argument. Try sqrt(complex(x)).",
      "",
      "Stacktrace:",
      " [1] \u001b[1msqrt\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Int64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./math.jl:434\u001b[22m\u001b[22m",
      " [2] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "sqrt(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0 + 1.0im"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqrt(-1 + 0im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why did Julia implement `sqrt` in this silly way, throwing an error for negative arguments unless you add a zero imaginary part?  Any reasonable person wants an imaginary result from `sqrt(-1)`, surely?\n",
    "\n",
    "The problem is that defining `sqrt` to return an imaginary result from `sqrt(-1)` would **not be type stable**: `sqrt(x)` would return a real result for non-negative real `x`, and a complex result for negative real `x`, so the **return type would depend on the value of `x`** and **not just its type.**\n",
    "\n",
    "That would defeat type inference, not just for the `sqrt` function, but for **anything the sqrt function touches**.  Unless the compiler can somehow figure out `x ≥ 0`, it will have to either store the result in a \"box\" or compile two branches of the result.  Let's see how that works by defining our own square-root function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mysqrt (generic function with 2 methods)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysqrt(x::Complex) = sqrt(x)\n",
    "mysqrt(x::Real) = x < 0 ? sqrt(complex(x)) : sqrt(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This definition is an example of Julia's [multiple dispatch style](http://docs.julialang.org/en/stable/manual/methods/), which in some sense is a generalization of object-oriented programming but focuses on \"verbs\" (functions) rather than nouns.  We will discuss this more in a later lecture.\n",
    "\n",
    "The `::Complex` and `::Real` are argument-type declarations.  Such declarations are **not related to performance**, but instead **act as a \"filter\"** to allow us to have one version of `mysqrt` for complex arguments and another for real arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4142135623730951"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysqrt(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0 + 1.4142135623730951im"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysqrt(-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0 + 1.4142135623730951im"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysqrt(-2+0im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks great, right?  But let's see what happens to type inference in a function that calls `mysqrt` instead of `sqrt`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variables:\n",
      "  #self#::#slowfun\n",
      "  x::Int64\n",
      "  #temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m\n",
      "  #temp#@_4::Core.MethodInstance\n",
      "  #temp#@_5\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m\n",
      "\n",
      "Body:\n",
      "  begin \n",
      "      $(Expr(:inbounds, false))\n",
      "      # meta: location In[6] mysqrt 2\n",
      "      unless (Base.slt_int)(x::Int64, 0)::Bool goto 6\n",
      "      #temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m = $(Expr(:invoke, MethodInstance for sqrt(::Complex{Float64}), :(Base.sqrt), :($(Expr(:new, Complex{Float64}, :((Base.sitofp)(Float64, x)::Float64), :((Base.sitofp)(Float64, 0)::Float64))))))\n",
      "      goto 8\n",
      "      6: \n",
      "      #temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m = (Base.Math.sqrt_llvm)((Base.sitofp)(Float64, x::Int64)::Float64)::Float64\n",
      "      8: \n",
      "      # meta: pop location\n",
      "      $(Expr(:inbounds, :pop))\n",
      "      unless (#temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m isa Complex{Float64})::Bool goto 14\n",
      "      #temp#@_4::Core.MethodInstance = MethodInstance for +(::Complex{Float64}, ::Int64)\n",
      "      goto 23\n",
      "      14: \n",
      "      unless (#temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m isa Float64)::Bool goto 18\n",
      "      #temp#@_4::Core.MethodInstance = MethodInstance for +(::Float64, ::Int64)\n",
      "      goto 23\n",
      "      18: \n",
      "      goto 20\n",
      "      20: \n",
      "      #temp#@_5\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m = (#temp#@_3\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m + 1)\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m\n",
      "      goto 25\n",
      "      23: \n",
      "      #temp#@_5\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m = $(Expr(:invoke, :(#temp#@_4), :(Main.+), :(#temp#@_3), 1))\n",
      "      25: \n",
      "      return #temp#@_5\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m\n",
      "  end\u001b[1m\u001b[91m::Union{Complex{Float64}, Float64}\u001b[39m\u001b[22m\n"
     ]
    }
   ],
   "source": [
    "slowfun(x) = mysqrt(x) + 1\n",
    "@code_warntype slowfun(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the compiler **doesn't know at compile-time that x is positive** (at compile-time it **uses only types, not values**, it doesn't know whether the result is real (`Float64`) or complex (`Complex{Float64}`) and has to store it in a \"box\".  This kills performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Int64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3665183052406099839"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i^100000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9223372036854775807"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i = typemax(Int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9223372036854775808"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining our own types\n",
    "\n",
    "Let's define our own type to represent a **\"point\" in two dimensions**.  Each point will have an $(x,y)$ location.  So that we can use the points with our `sum` functions above, we'll also define `+` and `zero` functions to do the obvious **vector addition**.\n",
    "\n",
    "One such definition in Julia is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(3, 4)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mutable struct Point1\n",
    "    x\n",
    "    y\n",
    "end\n",
    "Base.:+(p::Point1, q::Point1) = Point1(p.x + q.x, p.y + q.y)\n",
    "Base.zero(::Type{Point1}) = Point1(0,0)\n",
    "\n",
    "Point1(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(8, 10)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point1(3,4) + Point1(5,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our type is very generic, and can hold any type of `x` and `y` values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(3.7, 4 + 5im)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point1(3.7, 4+5im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps too generic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(\"x\", [3, 4, 5])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point1(\"x\", [3,4,5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since `x` and `y` can be *anything*, they must be **pointers to \"boxes\"**.  This is **bad news for performance**.\n",
    "\n",
    "A `mutable struct` is *mutable*, which means we can create a `Point1` object and then change `x` or `y`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(3, 4)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Point1(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(7, 4)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p.x = 7\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point1(7, 4)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that every reference to a `Point1` object must be a *pointer* to an object stored elsewhere in memory, because *how else would we \"know\" when an object changes?*  Furthermore, an **array of `Point1` objects must be an array of pointers** (which is **bad news for performance** again):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Point1,1}:\n",
       " Point1(7, 4)\n",
       " Point1(7, 4)\n",
       " Point1(7, 4)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = [p,p,p]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Point1,1}:\n",
       " Point1(7, 8)\n",
       " Point1(7, 8)\n",
       " Point1(7, 8)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Point1,1}:\n",
       " Point1(7, 8)\n",
       " Point1(7, 8)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q = [deepcopy(p), deepcopy(p)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Point1,1}:\n",
       " Point1(-3, 8)\n",
       " Point1(-3, 8)\n",
       " Point1(-3, 8)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p.x = -3\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Point1,1}:\n",
       " Point1(7, 8)\n",
       " Point1(7, 8)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test this out by creating an array of `Point1` objects and summing it.  Ideally, this would be about twice as slow as summing an equal-length array of numbers, since there are twice as many numbers to sum.  But because of all of the boxes and pointer-chasing, it should be far slower.\n",
    "\n",
    "To create the array, we'll call the `Point1(x,y)` constructor with our array `a`, using Julia's [\"dot-call\" syntax](http://docs.julialang.org/en/stable/manual/functions/#dot-syntax-for-vectorizing-functions) that applies a function \"element-wise\" to arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000000-element Array{Point1,1}:\n",
       " Point1(0.45576, 0.45576)    \n",
       " Point1(0.199227, 0.199227)  \n",
       " Point1(0.664066, 0.664066)  \n",
       " Point1(0.106669, 0.106669)  \n",
       " Point1(0.108969, 0.108969)  \n",
       " Point1(0.146101, 0.146101)  \n",
       " Point1(0.066834, 0.066834)  \n",
       " Point1(0.933081, 0.933081)  \n",
       " Point1(0.901608, 0.901608)  \n",
       " Point1(0.680544, 0.680544)  \n",
       " Point1(0.891394, 0.891394)  \n",
       " Point1(0.128344, 0.128344)  \n",
       " Point1(0.354404, 0.354404)  \n",
       " ⋮                           \n",
       " Point1(0.838992, 0.838992)  \n",
       " Point1(0.940758, 0.940758)  \n",
       " Point1(0.416696, 0.416696)  \n",
       " Point1(0.17713, 0.17713)    \n",
       " Point1(0.493455, 0.493455)  \n",
       " Point1(0.44981, 0.44981)    \n",
       " Point1(0.0298092, 0.0298092)\n",
       " Point1(0.894299, 0.894299)  \n",
       " Point1(0.270014, 0.270014)  \n",
       " Point1(0.32509, 0.32509)    \n",
       " Point1(0.191285, 0.191285)  \n",
       " Point1(0.643986, 0.643986)  "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = rand(10^7)\n",
    "a1 = Point1.(a, a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools, Compat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  517.472 ms (29999997 allocations: 610.35 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Point1(5.000915878780597e6, 5.000915878780597e6)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime sum($a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time is at least **50× slower** than we would like, but consistent with our other timing results on \"boxed\" values from last lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An imperfect solution: A concrete immutable type\n",
    "\n",
    "We can avoid these two problems by:\n",
    "\n",
    "* Declare the types of `x` and `y` to be *concrete* types, so that they don't need to be pointers to boxes.\n",
    "* Declare our Point to be an [immutable](https://en.wikipedia.org/wiki/Immutable_object) type (`x` and `y` cannot change), so that Julia is not forced to make every reference to a Point into a pointer: just `struct`, not `mutable struct`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point2(3.0, 4.0)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "struct Point2\n",
    "    x::Float64\n",
    "    y::Float64\n",
    "end\n",
    "Base.:+(p::Point2, q::Point2) = Point2(p.x + q.x, p.y + q.y)\n",
    "Base.zero(::Type{Point2}) = Point2(0.0,0.0)\n",
    "\n",
    "Point2(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point2(8.0, 10.0)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point2(3,4) + Point2(5,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Point2,1}:\n",
       " Point2(3.0, 4.0)\n",
       " Point2(3.0, 4.0)\n",
       " Point2(3.0, 4.0)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Point2(3,4)\n",
    "P = [p,p,p]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mtype Point2 is immutable\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mtype Point2 is immutable\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "p.x = 6 # gives an error since p is immutable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.0"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sizeof(P) / length(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this is working as we hope, then summation should be much faster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  11.792 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Point2(5.000915878780597e6, 5.000915878780597e6)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a2 = Point2.(a,a)\n",
    "@btime sum($a2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the time is **only about 10ms**, only slightly more than twice the cost of summing an array of individual numbers of the same length!\n",
    "\n",
    "Unfortunately, we paid a big price for this performance: our `Point2` type only works with *a single numeric type* (`Float64`), much like a C implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The best of both worlds: Parameterized immutable types\n",
    "\n",
    "How do we get a `Point` type that works for *any* type of `x` and `y`, but at the same time allows us to have an array of points that is concrete and homogeneous (every point in the array is forced to be the same type)?  At first glance, this seems like a contradiction in terms.\n",
    "\n",
    "The answer is not to define a *single* type, but rather to **define a whole family of types** that are *parameterized* by the type `T` of `x` and `y`.  In computer science, this is known as [parametric polymorphism](https://en.wikipedia.org/wiki/Parametric_polymorphism).  (An example of this can be found in [C++ templates](https://en.wikipedia.org/wiki/Template_%28C%2B%2B%29).)\n",
    "\n",
    "In Julia, we will define such a family of types as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point3{Int64}(3, 4)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "struct Point3{T<:Real}\n",
    "    x::T\n",
    "    y::T\n",
    "end\n",
    "Base.:+(p::Point3, q::Point3) = Point3(p.x + q.x, p.y + q.y)\n",
    "Base.zero{T}(::Type{Point3{T}}) = Point3(zero(T),zero(T))\n",
    "\n",
    "Point3(3,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `Point3` is actually a family of subtypes `Point{T}` for different types `T`.   The notation `<:` in Julia means \"is a subtype of\", and hence `T<:Real` means that we are constraining `T` to be a `Real` type (a built-in *abstract type* in Julia that includes e.g. integers or floating point)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point3{Float64}(8.6, 11.8)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Point3(3,4) + Point3(5.6, 7.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now, let's make an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000000-element Array{Point3{Float64},1}:\n",
       " Point3{Float64}(0.45576, 0.45576)    \n",
       " Point3{Float64}(0.199227, 0.199227)  \n",
       " Point3{Float64}(0.664066, 0.664066)  \n",
       " Point3{Float64}(0.106669, 0.106669)  \n",
       " Point3{Float64}(0.108969, 0.108969)  \n",
       " Point3{Float64}(0.146101, 0.146101)  \n",
       " Point3{Float64}(0.066834, 0.066834)  \n",
       " Point3{Float64}(0.933081, 0.933081)  \n",
       " Point3{Float64}(0.901608, 0.901608)  \n",
       " Point3{Float64}(0.680544, 0.680544)  \n",
       " Point3{Float64}(0.891394, 0.891394)  \n",
       " Point3{Float64}(0.128344, 0.128344)  \n",
       " Point3{Float64}(0.354404, 0.354404)  \n",
       " ⋮                                    \n",
       " Point3{Float64}(0.838992, 0.838992)  \n",
       " Point3{Float64}(0.940758, 0.940758)  \n",
       " Point3{Float64}(0.416696, 0.416696)  \n",
       " Point3{Float64}(0.17713, 0.17713)    \n",
       " Point3{Float64}(0.493455, 0.493455)  \n",
       " Point3{Float64}(0.44981, 0.44981)    \n",
       " Point3{Float64}(0.0298092, 0.0298092)\n",
       " Point3{Float64}(0.894299, 0.894299)  \n",
       " Point3{Float64}(0.270014, 0.270014)  \n",
       " Point3{Float64}(0.32509, 0.32509)    \n",
       " Point3{Float64}(0.191285, 0.191285)  \n",
       " Point3{Float64}(0.643986, 0.643986)  "
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a3 = Point3.(a,a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array{Point3{Float64},1}"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(a3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the type of this array is `Array{Point3{Float64},1}` (we could equivalently write this as `Vector{Point3{Float64}}`, since `Vector{T}` is a synonym for `Array{T,1}`).  You should learn a few things from this:\n",
    "\n",
    "* An `Array{T,N}` in Julia is itself a parameterized type, parameterized by the element type `T` and the dimensionality `N`.\n",
    "\n",
    "* Since the element type `T` is encoded in the `Array{T,N}` type, the element type does not need to be stored in each element.  That means that the `Array` is free to store an array of \"inlined\" elements, rather than an array of pointers to boxes.  (This is why `Array{Float64,1}` earlier could be stored in memory like a C `double*`.\n",
    "\n",
    "* It is still important that the element type be `immutable`, since an array of mutable elements would still need to be an array of pointers (so that it could \"notice\" if another reference to an element mutates it)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  11.805 ms (0 allocations: 0 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Point3{Float64}(5.000915878780597e6, 5.000915878780597e6)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime sum($a3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hooray! It is again **only about 10ms**, the same time as our completely concrete and inflexible `Point2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point3{Int8}(3, 5)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Point3{Int8}(3,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Point3{Integer}(3, 5)"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Point3{Integer}(3,Int8(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Int64"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(p.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Int8"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(p.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3 + 7im"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3 + 7im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Complex{Int64}"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "typeof(3 + 7im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "sqrt(z::<b>Complex</b>) at <a href=\"https://github.com/JuliaLang/julia/tree/903644385b91ed8d95e5e3a5716c089dd1f1b08a/base/complex.jl#L416\" target=\"_blank\">complex.jl:416</a>"
      ],
      "text/plain": [
       "sqrt(z::Complex) in Base at complex.jl:416"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@which sqrt(3 + 4im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Any"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "supertype(Point3{Int8})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another type of \"type instability\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mysum_slow (generic function with 1 method)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum_slow(a)\n",
    "    s = 0\n",
    "    for x in a\n",
    "        s += x\n",
    "    end\n",
    "    return s\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysum_slow([3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.3"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mysum_slow([3.3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variables:\n",
      "  #self#::#mysum_slow\n",
      "  a::Array{Float64,1}\n",
      "  x::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",
      "      x::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 + x::Float64)::Float64\n",
      "      goto 25\n",
      "      23: \n",
      "      #temp#@_7::Float64 = $(Expr(:invoke, :(#temp#@_6), :(Main.+), :(s), :(x)))\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(mysum_slow([3.3,4,5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mysum_fast (generic function with 1 method)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mysum_fast(a)\n",
    "    s = zero(eltype(a))\n",
    "    for x in a\n",
    "        s += x\n",
    "    end\n",
    "    return s\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variables:\n",
      "  #self#::#mysum_fast\n",
      "  a::Array{Float64,1}\n",
      "  x::Float64\n",
      "  #temp#::Int64\n",
      "  s::Float64\n",
      "\n",
      "Body:\n",
      "  begin \n",
      "      s::Float64 = (Base.sitofp)(Float64, 0)::Float64 # line 3:\n",
      "      #temp#::Int64 = 1\n",
      "      4: \n",
      "      unless (Base.not_int)((#temp#::Int64 === (Base.add_int)((Base.arraylen)(a::Array{Float64,1})::Int64, 1)::Int64)::Bool)::Bool goto 14\n",
      "      SSAValue(2) = (Base.arrayref)(a::Array{Float64,1}, #temp#::Int64)::Float64\n",
      "      SSAValue(3) = (Base.add_int)(#temp#::Int64, 1)::Int64\n",
      "      x::Float64 = SSAValue(2)\n",
      "      #temp#::Int64 = SSAValue(3) # line 4:\n",
      "      s::Float64 = (Base.add_float)(s::Float64, x::Float64)::Float64\n",
      "      12: \n",
      "      goto 4\n",
      "      14:  # line 6:\n",
      "      return s::Float64\n",
      "  end::Float64\n"
     ]
    }
   ],
   "source": [
    "@code_warntype mysum_fast([3.3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mMethodError: no method matching zero(::Type{Array{Float64,2}})\u001b[0m\nClosest candidates are:\n  zero(\u001b[91m::Type{Base.LibGit2.GitHash}\u001b[39m) at libgit2/oid.jl:106\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}\u001b[39m) at pkg/resolve/versionweight.jl:82\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}\u001b[39m) at pkg/resolve/versionweight.jl:124\n  ...\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mMethodError: no method matching zero(::Type{Array{Float64,2}})\u001b[0m\nClosest candidates are:\n  zero(\u001b[91m::Type{Base.LibGit2.GitHash}\u001b[39m) at libgit2/oid.jl:106\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}\u001b[39m) at pkg/resolve/versionweight.jl:82\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}\u001b[39m) at pkg/resolve/versionweight.jl:124\n  ...\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mmysum_fast\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Array{Float64,2},1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./In[75]:2\u001b[22m\u001b[22m",
      " [2] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "mysum_fast([rand(3,3) for i = 1:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 5.11819  4.99643  3.24856\n",
       " 5.22759  5.39936  5.74627\n",
       " 6.72741  4.78524  4.62801"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([rand(3,3) for i = 1:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mMethodError: no method matching zero(::Type{Array{Float64,2}})\u001b[0m\nClosest candidates are:\n  zero(\u001b[91m::Type{Base.LibGit2.GitHash}\u001b[39m) at libgit2/oid.jl:106\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}\u001b[39m) at pkg/resolve/versionweight.jl:82\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}\u001b[39m) at pkg/resolve/versionweight.jl:124\n  ...\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mMethodError: no method matching zero(::Type{Array{Float64,2}})\u001b[0m\nClosest candidates are:\n  zero(\u001b[91m::Type{Base.LibGit2.GitHash}\u001b[39m) at libgit2/oid.jl:106\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}\u001b[39m) at pkg/resolve/versionweight.jl:82\n  zero(\u001b[91m::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}\u001b[39m) at pkg/resolve/versionweight.jl:124\n  ...\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mmr_empty\u001b[22m\u001b[22m at \u001b[1m./reduce.jl:243\u001b[22m\u001b[22m [inlined]",
      " [2] \u001b[1m_mapreduce\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Base.#identity, ::Base.#+, ::IndexLinear, ::Array{Array{Float64,2},1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./reduce.jl:265\u001b[22m\u001b[22m",
      " [3] \u001b[1msum\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Array{Float64,2},1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./reduce.jl:359\u001b[22m\u001b[22m",
      " [4] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "sum(Matrix{Float64}[])"
   ]
  }
 ],
 "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
}