Chapter3_Dispatch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TUTORIAL:\n",
    "# Julia for technical computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 3. Multiple dispatch\n",
    "\n",
    "What is multiple dispatch? **Dispatch** means that Julia will call a particular version of a function, based on the types of the arguments. **Multiple dispatch** means that Julia takes the types of all of the arguments into account, not just one. **Dynamic multiple dispatch** means that the same code will always result in the same function call, regardless of whether the Julia compiler was able to infer the types of the arguments at compile-time or had to do a check (*dynamically*) at run-time.\n",
    "\n",
    "It is through advanced **type inference** that Julia can often skip the run-time check. If functions unambiguously map the types of their input arguments to an output type, if they are *type stable*, then Julia can follow the flow of types through your code without actually executing the code.\n",
    "\n",
    "Of course we aim to avoid runtime type checks, because they slow the execution of our code. But we are not fanatic about it either. After all, **an interactive and dynamic programming language is mostly about convenience**. We do want to optimize the time-critical path in our code and Julia allows that, without resorting to another language for the parts that have to be fast. Outside of that path, by all means, use untyped fields and violate type-stability at will! Live your life the happy I-don't-know-about-types-and-I-do-not-want-to-know way, and enjoy it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we set out to illustrate the possibilities enabled by multiple dispatch in Julia. You will find little in this chapter that compares to Matlab or Python. We do elaborate on the difference between multiple dispatch and OOP concepts like method overloading and polymorphism. But that is my bias. Note that Julia is neither the first nor the only language to offer multiple dispatch, or *multi-methods*. Other examples include Common Lisp and Dylan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Dispatch and algorithm selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There are functions and methods. What's the difference?\n",
    "\n",
    "I did not tell you this before, but perhaps you've noticed: there are [functions](http://julia.readthedocs.org/en/release-0.4/manual/functions/) and there are [methods](http://julia.readthedocs.org/en/release-0.4/manual/methods/) in Julia. One function can have several methods. Each method corresponds to a combination of types of arguments.\n",
    "\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(x::Integer) = \"Acting on Integer\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(x::AbstractFloat) = \"Acting on AbstractFloat\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia tells me that `action` is a generic function with 2 methods. Here you see dispatch at work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can list all methods of a function, and you can also simply ask Julia which method it intends to invoke for a given set of arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(action)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which action(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which action(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Types match with their supertypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both methods we provided the argument had an abstract type. When Julia encounters a set of argument types and is looking for a function to dispatch to, an abstract type matches with all of its subtypes. However, if you have a specific implementation for a specific subtype, then you can provide a more specific method. Julia will dispatch to the most specific method available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(x::Int32) = \"For Int32 I have a better algorithm\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(action)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "action(Int32(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methods are algorithms are methods: Julia selects the best available algorithm for the problem at hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function overloading is possible in many typed languagues, but in Julia it is really pervasive. You add methods to existing functions all the time. Popular ones are the basic operators, like `+`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(+)   # careful, long output! Scroll down afterwards"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In my workspace at the time of writing, Julia reports that `+` has 155 methods. Each of them is optimized for a specific combination of types. You may have spotted a low-level method for the addition of two boolean values (`Bool`), as well as a high-level method for the addition of a lower-triangular and an upper-triangular matrix.\n",
    "Whenever you use `+`, or whenever you invoke another function that uses `+`, a specific optimized method is used. Even better, that specific method is compiled specifically for the types you are working with, enabling all sorts of optimizations. This all happens automatically and you don't have to think about it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will attempt to illustrate the power of this property throughout this chapter. The selection of the best algorithm does not happen just for `+`. It happens for all functions, all the time, for low-level as well as high-level operations. A good design of Julia code makes sure that different problems map to different types, or that different properties of a problem match to different types. You write the algorithms in terms of these types, and **Julia makes sure that the right algorithm is called at the right time**. As you see, this is not entirely for free. It is not really the case that Julia selects the right algorithm. Julia selects the method that best matches the type signature. It is up to you, the programmer, to make sure that this corresponds to the right algorithm.\n",
    "\n",
    "This does mean that **multiple dispatch affects the design of your code**. You design data structures with dispatch in mind, and you write algorithms that apply to types. This is not difficult to achieve once you're used to it, but it is quite pervasive and it is so pretty much from the start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An example: arithmetic modulo N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add our own method to the generic function `+`. We create a type of integers modulo N and define their addition.\n",
    "\n",
    "In the code below, we use a so-called *inner constructor*. The inner constructor is a function with the name of the type that is included in the definition of the type. The type parameter `N` is available inside the type definition. We use the constructor to make sure that field `n` is always betweeen `0` and `N-1`. Read more about constructors in the [manual](http://docs.julialang.org/en/release-0.4/manual/constructors/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "immutable ModInt{N}\n",
    "    n :: Int64\n",
    "    \n",
    "    ModInt(n) = new(mod(n,N))  # Inner constructor: make sure that 0 <= n < N\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a `ModInt` by explicitly calling the constructor with type parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = ModInt{7}(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, if you like numbers modulo 7 very much, you can define an alias just for this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typealias ModInt7 ModInt{7}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ModInt7(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enough syntax worries, that's not our focus: let's define addition. We include the parameter `N` in the definition of `+` for `ModInt`'s, so our method of `+` will be a parametric one. Note that the type parameter `N` is the same for the types of arguments `a` and `b`. The method will only apply to `ModInt`'s with the same parameter `N`. This is called **diagonal dispatch**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: +\n",
    "\n",
    "+{N}(a::ModInt{N}, b::ModInt{N}) = ModInt{N}(a.n+b.n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = ModInt{7}(4)\n",
    "b = ModInt{7}(5)\n",
    "a+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which a+b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There you are. Whenever you add two `ModInt`'s together anywhere in your code, this method will be called. Add a few other arithmetic operations, and chances are you can build matrices of `ModInt`'s and perform linear algebra operations.\n",
    "\n",
    "Julia does not know what to do when adding `ModInt`'s with different `N`, because we did not define it. (I would not know what to do mathematically in this case either.) The type signature does not match with any of the methods of `+`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ModInt{7}(2) + ModInt{5}(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's add another operation and see more of multiple dispatch at work. We can define multiplication by an integer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: *\n",
    "\n",
    "*{N}(a::Integer, b::ModInt{N}) = ModInt{N}(a*b.n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "4*ModInt{31}(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I have typed `a` to be an `Integer`. This means that in one line I have defined multiplication (only on the left!) with `Int64`, `Int32` and so on. If, for some reason, you have a better way of interacting with `Int32`, you can make a more specific method. Even with multiple arguments, Julia will select the most specific method, favouring subtypes over supertypes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function *{N}(a::Int32, b::ModInt{N})\n",
    "    println(\"Hi. I see you're multiplying with an Int32. Good, because I have a better algorithm for this case!\")\n",
    "    ModInt{N}(a*b.n)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "2 * ModInt{7}(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Int32(2) * ModInt{7}(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This was a silly example, but in case you *do* have a better example for a specific case, it is trivial to add it. Have a new structured matrix with a fast matrix-vector product? Define a type for the matrix, and define multiplication with a vector. That's it: your fast matrix-vector product will be used whenever `A*x` comes up. Hey, that's exactly what we will do in the next section! But first a final comment on dispatch."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Ambiguity warnings: what if different signatures match your argument types?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case there is ambiguity, the compiler will emit a warning. This often happens when duck-typing, i.e., if you don't specify the type of an argument. An untyped argument will match any type. For the purposes of dispatch, an untyped argument has type `Any`.\n",
    "\n",
    "Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "combine(a, b::ModInt{7}) = \"method 1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "combine(a::ModInt{7}, b) = \"method 2\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a warning (the first time you evaluate the code) because both type signatures match when `a` and `b` are both of type `ModInt`. There is no rule to prefer one method over the other, because neither is more specific than the other. If this case never happens at runtime then it will never be a problem. The compiler is simply saying it *could* happen. If it *does* happen, Julia will just pick one method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "combine(ModInt{7}(2), ModInt{7}(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution in this case is to do as the compiler suggested in the warning message: provide a definition for the more specific case. Since this definition is more specific, it takes precedence over the other two and the ambiguity is resolved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "combine(a::ModInt{7}, b::ModInt{7}) = \"method 3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "combine(ModInt{7}(2), ModInt{7}(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll encounter many ambiguity warnings soon enough. It is not always as simple to get rid of them, especially if one of your methods clashes with a method in somebody else's code.\n",
    "\n",
    "If you want to see what I mean, try the same definition of `ModInt` above, but make it a subtype of `Integer`. Now define some operations and try to get rid of all warnings. Good luck with that! (Hint: I've mentioned it before, but look into the [conversion-promotion](http://julia.readthedocs.org/en/release-0.4/manual/conversion-and-promotion/) system. Also, ask yourself: of course a `ModInt` is an integer, but is that sufficient reason to actually inherit from `Integer`?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2. Single dispatch versus multiple dispatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polymorphic member functions in OOP are like single dispatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Say in an object-oriented language you have an abstract or virtual `Matrix` class, and two subclasses called `DenseMatrix` and `LowRankMatrix`. Perhaps the `Matrix` class has a member function called `rank`, which returns the rank of the matrix. You would supply a definition of `rank` for `DenseMatrix` based on the singular value decomposition, for example, while the `LowRankMatrix` probably stores its rank as part of its representation and so `rank` would have a more efficient implementation than the general one.\n",
    "\n",
    "In any case, for an object called `mat` that is an instance of `Matrix`, you would probably write\n",
    "```\n",
    "mat.rank()\n",
    "```\n",
    "or\n",
    "```\n",
    "mat->rank()\n",
    "```\n",
    "and the version of `rank` that is invoked depends on the actual type of `mat` at run-time. This is **polymorphism**.\n",
    "\n",
    "In Julia you would typically write\n",
    "```\n",
    "rank(mat)\n",
    "```\n",
    "where the version of `rank` that is invoked depends on the actual type of `mat` at run-time. **This, too, is polymorphism**. It means exactly the same thing. It is just different syntax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An example: low-rank matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our examples are becoming more advanced as we go along. Just for fun, here is a simple definition of a low rank matrix with a fast matrix-vector product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: *, size, rank\n",
    "\n",
    "# A LowRankMatrix is an MxN matrix of the form A*B^t, where A is MxK and B is NxK.\n",
    "# K is the rank of the matrix.\n",
    "# It is good practice to include the element type as a type parameter: this makes\n",
    "# the types of a and b concrete (and therefore efficient), but still general.\n",
    "immutable LowRankMatrix{T}\n",
    "    a :: Array{T,2}\n",
    "    b :: Array{T,2}\n",
    "end\n",
    "\n",
    "size(A::LowRankMatrix) = ( size(A.a,1), size(A.b,1) )\n",
    "\n",
    "*{T}(A::LowRankMatrix{T}, x::Array{T,1}) = A.a * A.b' * x\n",
    "\n",
    "rank(A::LowRankMatrix) = size(A.a, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did you notice how short this code is? We leave it as an exercise to add proper dimensionality checks, for example in an inner constructor and using [assertions](http://julia.readthedocs.org/en/release-0.4/manual/metaprogramming/#building-an-advanced-macro). The code is also not optimal, because the matrix-vector product allocates memory. Still, the matrix-vector product has the right computational complexity: $O((N+M)K)$ instead of $O(NM)$ for multiplication with a dense $M \\times N$ matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = LowRankMatrix(rand(10,2), rand(10,2))\n",
    "size(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rank(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A*rand(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Arithmetics with low-rank matrices\n",
    "\n",
    "It may make sense to you to have `rank` be a member function of `mat`. Clearly, the rank of a matrix is an inherent property of that matrix. But things start to look differently if you mix several objects. What about addition?\n",
    "\n",
    "In the OOP statement\n",
    "```\n",
    "mat.add(mat2)\n",
    "```\n",
    "`mat` is more special than `mat2`, whereas in the Julia formulation\n",
    "```\n",
    "add(mat, mat2)\n",
    "```\n",
    "their roles are entirely equivalent. It is quite clear that the best algorithm to add two matrices is determined by the types of both matrices, not just the first one. Their order should not matter.\n",
    "\n",
    "In several typed languages, using function overloading you can define addition for two matrices outside of a class definition. That is how it works in Julia too, using multiple dispatch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: +\n",
    "\n",
    "# We simply concatenate to create a matrix of rank k1+k2, where k1=rank(mat1) and k2=rank(mat2)\n",
    "+(mat1::LowRankMatrix, mat2::LowRankMatrix) = LowRankMatrix([mat1.a mat2.a], [mat2.b mat2.b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = LowRankMatrix(rand(10,2), rand(10,2))\n",
    "B = LowRankMatrix(rand(10,2), rand(10,2))\n",
    "C = A+B\n",
    "rank(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily interact with the outside world and express that a low rank matrix times any other matrix is always a low rank matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: *\n",
    "\n",
    "*(mat1::LowRankMatrix, mat2::AbstractMatrix) = LowRankMatrix(mat1.a, (mat1.b'*mat2)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then we can create a more specific way to multiply low rank matrices. Though the more general routine above would work for this case too, the definition below is *slightly* more efficient because it uses fewer transposes in total. Since we know this special case, we plug it in, and Julia will make sure it is used on every occasion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "*(mat1::LowRankMatrix, mat2::LowRankMatrix) = LowRankMatrix(mat1.a, mat1.b*mat2.a'*mat2.b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which A*rand(10,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which A*B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rank(A*B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Static versus dynamic multiple dispatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Dynamic means things can change\n",
    "\n",
    "A language like C++ heavily supports *generic programming* with templates and function overloading. Yet, all of that is static: all types are known at compile-time.\n",
    "\n",
    "Julia is dynamic, and this has some advantages. Here is an example, going back to integers modulo a number `N`. Remember we had an awkward constructor before. We had to fully qualify the type parameters and write `ModInt{7}(2)`. But we can also define this function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "modint(n, N) = ModInt{N}(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = modint(2, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `modint` creates a new integer modulo `N`. Though `N` is a type parameter for `ModInt`, it is just an argument to the `modint` function. For all we know, its value could come from user input. There is no way to statically determine what the full type of the `ModInt` is going to be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Strictly speaking, this makes `modint` type-unstable (if that is even a word). The type of the output depends on a value of the input, not just on the input types. Julia will not know at compile-time what the type of the output will be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@code_typed modint(Int64,Int64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the last line, where it says `end::ModInt{Int64}`. The compiler can deduce that `N` is going to be an integer, but it does not know which one it is. Compare this to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@code_typed ModInt{7}(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the compiler has full knowledge about the type at the end."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `modint` function is convenient: we can create `ModInt`'s on the fly for any `N` the user would like. But type-inference has a problem tracking the flow of types through the code. I've warned against type instabilities. Is this going to be a problem? Well, yes and no. *One type instability does not a disaster make.* Read on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The effect of type-instability is always local: dynamic dispatch helps type inference recover\n",
    "\n",
    "That is one complicated title, but here is the thing, and it is rather fundamental. Perhaps Julia can not infer the type of one of your variables. Pity, expect things to slow down. Yet, as soon as you invoke a function, dynamic dispatch guarantees to call the most specific method available. Julia checks the types of the arguments and compiles a version of the method specifically for those types. If necessary, this happens at runtime. **Aha!** But that means that **inside the function we've called, the compiler knows all the types again** during compilation! Type inference can continue on its merry way.\n",
    "\n",
    "This definitely calls for an example, because you don't have to take my word for it. Let's have a type with an untyped field, which from its name we expect to be an array of some kind."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type VerySlowType\n",
    "    some_array\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps we want to sum over the elements in `some_array`. Here is one way to do it. Admittedly it is a strange way (see also our earlier discussion of sums in the previous chapter) but that is not the point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Base: sum\n",
    "\n",
    "function sum(v::VerySlowType)\n",
    "    z = v.some_array[1]\n",
    "    for i = 2:length(v.some_array)\n",
    "        z = z + v.some_array[i]\n",
    "    end\n",
    "    z\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v = VerySlowType( rand(10000))\n",
    "typeof(v.some_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sum(v)\n",
    "@time sum(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is the summation so slow? Julia has JIT-compiled `sum` for arguments of type `VerySlowType`. It has no way of knowing what the type of `v.slow_array` will be, so it emits general and slow code that can deal with any type. That is what is going on in the for loop. It is not optimized for arrays of `Float64`'s.\n",
    "\n",
    "The Julia miracle is that one indirection completely solves this problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sum2(v::VerySlowType) = sum_some_array(v.some_array)\n",
    "\n",
    "function sum_some_array(some_array)\n",
    "    z = some_array[1]\n",
    "    for i = 2:length(some_array)\n",
    "        z = z + some_array[i]\n",
    "    end\n",
    "    z\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v = VerySlowType( rand(10000))\n",
    "sum2(v)\n",
    "@time sum2(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There you are: blazingly fast and zero allocations! With almost exactly the same source code. Why? **Because `sum_some_array` knows the type of `some_array`, even though `sum2` does not!** The call to `sum_some_array` from `sum2` incurs a dynamic performance hit because the type of `v.some_array` is unknown. But this type will be determined and Julia sees it is an array of `Float64`'s. It will call a version of `sum_some_array` that is optimized for arrays of `Float64`'s.\n",
    "\n",
    "There is a one-time performance hit for compiling `sum_some_array`. There is also a recurring performance hit for looking up the type of `some_array` at runtime, each time we invoke `sum2`, because we call another function with `some_array` as an argument. But that is limited to just one function call, and the time it takes to resolve is dwarfed by the execution time of the sum. The important thing to optimize here is the summation itself, and that is what with a single indirection we have achieved.\n",
    "\n",
    "This is perhaps a subtle point to make, but it is a very important one. Go back to the example until you've convinced yourself that you've completely understood. It brings you that much closer to appreciating where the speed of Julia comes from and how you can write efficient code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 4. Some design patterns with dispatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data or metadata: do I use a type parameter or an extra field?\n",
    "\n",
    "Remember the `ModInt` type I defined earlier in this chapter, for integers modulo `N`? It had the number `N` as a type parameter. It exists only in the compiler. I could also have used this definition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "immutable ModInt2\n",
    "    n :: Int64\n",
    "    N :: Int64\n",
    "    \n",
    "    ModInt2(n, N) = new(mod(n,N), N)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are advantages and disadvantages. With the second definition, all numbers would take twice as much memory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sizeof(ModInt2(5, 7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sizeof(ModInt{7}(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, I could efficiently use thousands of different values of `N` in one program. When `N` is a type parameter as in `ModInt{N}`, each different value of `N` results in a different type. This means all functions and methods that use `ModInt{N}` will be recompiled, adding to compilation overhead and size of code.\n",
    "\n",
    "Another consideration: when `N` is a type parameter, I can dispatch on the value of `N`. Not so when `N` is a field, because you can not dispatch on the value of a variable, only on its type. I would be forced to write things like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function +(a::ModInt2, b::ModInt2)\n",
    "    if a.N == b.N\n",
    "        ModInt2(a+b, a.N)\n",
    "    else\n",
    "        error(\"No can do.\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Such kinds of `if` statements at the start of a function are a clear indication that dispatch may be a more elegant alternative.\n",
    "\n",
    "There is a trade-off between speed of execution and compilation time. In the case of integers modulo `N`, having a type parameter `N` seems like the most defensible option. But a different choice was made for the sizes of vectors and matrices. Julia includes the dimension of an array as a type parameter, but not its actual size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = Array(Float64, 10,10)\n",
    "typeof(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "size(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine the amount of recompilation if each function using matrices depended on the size of the matrix! Usually, the cost associated with verifying correct dimensions is dwarfed by the cost of the actual operation on the matrix anyway. Removing the `if` statement will not make much difference.\n",
    "\n",
    "The overhead may be significant for small matrices though. And, sure enough, there is a packaged called [`FixedSizeArrays`](https://github.com/SimonDanisch/FixedSizeArrays.jl) which includes the size of an array in its type and leads to better performance for small arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Efficient method selection with singleton types\n",
    "\n",
    "You often have several options for solving a particular numerical problem. It is customary in other languages to provide options with integer arguments, or even strings. An example could be `dwt(x, L, \"db2\")` to compute the discrete wavelet transform of a vector `x` to `L` levels using the Daubechies 2 wavelet (represented by the string \"db2\").\n",
    "\n",
    "Singleton types carry no data, but their type is useful metadata. They can be used to select methods.\n",
    "\n",
    "Compare this definition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function dwt(x, L, name)\n",
    "    if name == \"db2\"\n",
    "        \"Computing a wavelet transform using Daubechies 2\"\n",
    "    elseif name == \"db4\"\n",
    "        \"Computing a wavelet transform using Daubechies 4\"\n",
    "    else\n",
    "        \"Using a generic algorithm for some general wavelet\"\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "to the following one, which uses singleton types to select the right transform:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "abstract DiscreteWavelet\n",
    "\n",
    "immutable Daubechies{N} <: DiscreteWavelet\n",
    "end\n",
    "\n",
    "immutable SomeOtherWavelet <: DiscreteWavelet\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dwt{T <: DiscreteWavelet}(x, L, ::Type{T}) = \"Using a generic algorithm for some general wavelet\"\n",
    "\n",
    "dwt{N}(x, L, ::Type{Daubechies{N}}) = \"Computing a wavelet transform using Daubechies $N\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first method of `dwt` illustrates how to write code that applies to any wavelet: the type parameter `T` is restricted to be a subtype of DiscreteWavelet.\n",
    "\n",
    "The second method may use a more specific algorithm for a Daubechies wavelet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dwt(rand(10), 4, SomeOtherWavelet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dwt(rand(10), 4, Daubechies{2})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One important thing to remark is that option selection using singleton types makes your code type-stable. Each different singleton type yields a different type signature of the method. So, it is okay if the output type is different, for example due to using a different algorithm. (That was not the case above because I simply returned strings, but you get the point.)\n",
    "\n",
    "Other advantages include:\n",
    "* increased probability of inlining: no list of conditional statements\n",
    "* ease of adding options: you don't have to modify the `if elseif` list in `dwt`, you simply define a new wavelet and add a method for it\n",
    "* improved performance: your code jumps immediately to the right function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Is value-based dispatch a good idea?\n",
    "\n",
    "Once you're used to dispatch based on types, you may become greedy and notice some of your code would become faster if you could also dispatch on values. It may allow you to remove even more `if` and `case` statements.\n",
    "\n",
    "Tracking values of variables is orders of magnitude harder than tracking types. In Julia, I don't think it will happen. The performance overhead for analyzing possible values a variable might have would be immense. The only advantage seems to be that you can more easily add special cases for special values, for example to define the limit of a function at a special point.\n",
    "\n",
    "One thing that is provided is the [`Val{T}`](http://docs.julialang.org/en/release-0.4/manual/types/#value-types) type. The type parameter `T` can be any [*bitstype*](http://docs.julialang.org/en/release-0.4/manual/types/#bits-types), including `Float64`, `Bool`, `Int64` and so on. Here, the compiler knows the numerical value of `T`, because it is part of the type `Val{T}`. You can do computations with `T` at compile-time. The trade-off is that each numerical value for `T` results in a new type and, hence, in more code and compilation time.\n",
    "\n",
    "Here is an example where `Val{T}` seems to make sense. Let's say you want a routine for computing the $l^p$ norm of a vector. You could write:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function norm(x, p)\n",
    "    if p == 1\n",
    "        sum(abs(x))\n",
    "    elseif p == Inf\n",
    "        maximum(abs(x))\n",
    "    elseif p == 0\n",
    "        countnz(x)  ## number of non-zero entries\n",
    "    else\n",
    "        sum(abs(x).^p)^(1/p)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conditional statements make this routine less readable, and more difficult to inline. However, you are unlikely to use many different values of $p$. The important ones are mainly 0, 1, 2 and $\\infty$. Here is a more elegant solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# general definition\n",
    "norm{P}(x, ::Type{Val{P}}) = sum(abs(x).^P).^(1/P)\n",
    "\n",
    "# l0 norm: count number of non-zero entries\n",
    "norm(x, ::Type{Val{0}}) = countnz(x)\n",
    "\n",
    "# l1 norm: sum of absolute values\n",
    "norm(x, ::Type{Val{1}}) = sum(abs(x))\n",
    "\n",
    "# l∞ norm: maximum of absolute value\n",
    "norm(x, ::Type{Val{Inf}}) = maximum(abs(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = rand(10)\n",
    "norm(a, Val{2})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "norm(a, Val{1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "norm(a, Val{Inf})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which norm(a, Val{1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which norm(a, Val{Inf})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which norm(a, Val{3})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.0",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}