1. Starting out.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: The Newton method for finding roots of functions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton method is an iterative method to solve equations of the form $f(x)=0$, i.e. to find *roots* or *zeros* $x^\\ast$ such that $f(x^\\ast) = 0$. Given an initial guess $x_0$, we repeat the iteration\n",
    "\n",
    "$$x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_n)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement the Newton algorithm in Julia. We start from an initial condition $x_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_0 = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia always returns a value from any expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use LaTeX notation and tab completion for Unicode, e.g. `x\\_0<TAB>`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x₀ = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x₀"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Values in Julia have associated **types**. We can find the type of a variable using the appropriately-named `typeof` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(x₀)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can guess that this means an integer with 64 bits. [This result will be `Int32` if you have a 32-bit machine.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple functions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to define a function whose roots we wish to find. Let's find square roots of two, for example. Julia provides a concise mathematical syntax for defining simple functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(x) = x^2 - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(x_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need the derivative function, $f'$. For the moment, let's just give it my hand. (Later, we will see a neat way to avoid this.) We may like to write `f'`, using the apostrophe, `'`), but the apostrophe turns out to be a special character in Julia, so we get an error if we try to define a variable or function named `f'`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f'(x) = 2x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, Unicode comes to our rescue: `f\\prime<TAB>`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f′(x) = 2x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can do one step of our algorithm; mathematical operations work like we expect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_1 = x_0 - f(x_0) / f'(x_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that division of integers using `/` gives a floating-point result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(x_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iteration "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to repeat such steps several times. Julia has `for` loops and `while` loops. As usual, we tend to use `for` loops when we know how many iterations we want, and `while` when we iterate until a certain condition is attained. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ranges and arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with a simple `for`. Blocks of code in Julia *always* end with `end`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in 1:5\n",
    "    println(i)  # print the value of i followed by a new line\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, a variable `i` is introduced that is *local* to the loop, i.e. it exists only inside the loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`i` takes each value in the *iterable collection* `1:5`. Let's ask Julia what this object `1:5` is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "1:5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, Julia returns a value, but in this case it is (at first glance) apparently not very helpful. What type is this object?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(1:5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that Julia has a special type (actually, several different types) to represent **ranges**, in which the elements are calculated each time a new element is required, rather than stored. We can see all the elements that will be produced using the `collect` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "v = collect(1:5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an object of a new type, an `Array`, in this case one whose elements are integers and that is of *dimension* 1. Note that `1` is *not* the number of elements in the array, which is called `length`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "length(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Array`s are also iterable, so we can iterate over an `Array` using a `for` loop. 1-dimensional arrays, also called `Vector`s, are constructed using square brackets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "w = [3, 4, 7]\n",
    "for i in w\n",
    "    println(2*i)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing the Newton method "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to implement the Newton method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_0 = 3\n",
    "x = x_0\n",
    "\n",
    "for i in 1:10\n",
    "    x_new = x - f(x) / f′(x)\n",
    "    println(i, \"\\t\", x_new)\n",
    "    x = x_new\n",
    "end\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we see that the method rapidly converges to one of the square roots of two. Which root it converges to depends on the initial condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_0 = -3\n",
    "x = x_0\n",
    "\n",
    "for i in 1:10\n",
    "    x_new = x - f(x) / f′(x)\n",
    "    println(i, \"\\t\", x_new)\n",
    "    x = x_new\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton method is, in fact, not guaranteed to converge to a root (although it always does so if started \"sufficiently close\" to a root, at a rate that is known). Furthermore, *which* root it converges to can depend sensitively on the initial condition. Let's calculate this for several initial conditions.\n",
    "\n",
    "First we create a set of initial conditions on the real line, say between -5 and 5. We now include a step size in the range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initial_conditions = -5:0.1:5\n",
    "collect(initial_conditions)   # use tab completion for long variable names!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This range type is different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(-5:0.1:5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The array is also a new type: it is now an array of 64-bit floating-point numbers. We can also see that the `{...}` syntax thus gives the **parameters** of the `Array` type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each of these initial conditions, we will run the Newton algorithm for a certain number of steps and store the resulting value. We thus need a new array in which to store the results. One way of creating an array is using the `similar` function, which, by default, creates an array of the same type and same size, but with (currently) uninitialized values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "roots = similar(initial_conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we do the work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (j, x_0) in enumerate(initial_conditions)\n",
    "    x = x_0\n",
    "\n",
    "    for i in 1:100\n",
    "        x = x - f(x) / f′(x)\n",
    "    end\n",
    "    \n",
    "    roots[j] = x\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `enumerate` iterates over `initial_conditions` but returns not only the value at each step, but also a counter. `(j, x_0`) is called a **tuple** (an ordered pair):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t = (3, 4)\n",
    "typeof(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NB: In Julia v0.4, tuples have been completely reworked, and the resulting type is now\n",
    "`Tuple{Int64,Int64}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia does not show all of the contents of an array by default. We can see everything using `showall`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "showall(roots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, apart from the `NaN` value, the results are not very exciting. Let's work harder with more initial conditions. We can find out how long the calculation takes using `@time` by wrapping the code in a `begin...end` block:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@time begin \n",
    "    initial_conditions = -100:0.01:100\n",
    "    roots = similar(initial_conditions)\n",
    "\n",
    "    for (j, x_0) in enumerate(initial_conditions)\n",
    "        x = x_0\n",
    "\n",
    "        for i in 1:1000\n",
    "            x = x - f(x) / f′(x)\n",
    "        end\n",
    "\n",
    "        roots[j] = x\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages and visualisation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are now many values stored in the array, so it is hopeless to examine them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "length(roots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead, we turn to **visualisation**. There are several plotting **packages** in Julia: [`Gadfly`]() is a native Julia library that produces beautiful plots; [`PyPlot`] is a Julian interface to the well-known `matplotlib` Python library. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with `PyPlot`. First we need to download the package. Julia provides a built-in package manager, called `Pkg`, that gracefully handles dependencies, etc. To tell Julia that we require the package, we do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Pkg.add(\"PyPlot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step is necessary only once. In each session where we need to use `PyPlot` we do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this process of loading a package currently can take a considerable time. Work is in progress to reduce this loading time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "plot(roots);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are used to the performance of C or Fortran, we might start to be unhappy with Julia's speed in this rather simple calculation. A close inspection of the output of the `@time` operation, however, gives us a very important clue: Julia apparently allocated over a gigabyte of memory to do a simple loop with some floating-point numbers! \n",
    "\n",
    "This is almost *always* a very strong signal that there is something very wrong in your Julia code! In our case, it is not at all clear what that could be. It turns out to be something very fundamental in Julia: \n",
    "\n",
    "[almost] **NEVER WORK WITH GLOBAL OBJECTS!**\n",
    "\n",
    "Due to technical details about the way that Julia works, it turns out that **GLOBALS ARE BAD**. What is the solution? **PUT EVERYTHING INTO A FUNCTION!** Let's try following this advice. We take *exactly* the same code and just plop it into a new function. For longer functions, Julia has an alternative syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function do_roots()\n",
    "    initial_conditions = -100:0.01:100\n",
    "    roots = similar(initial_conditions)\n",
    "\n",
    "    for (j, x_0) in enumerate(initial_conditions)\n",
    "        x = x_0\n",
    "\n",
    "        for i in 1:1000\n",
    "            x = x - f(x) / f′(x)\n",
    "        end\n",
    "\n",
    "        roots[j] = x\n",
    "    end\n",
    "    \n",
    "    roots\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the last line of the function. This will automatically *return* the value of the `roots` object as the output of the function. So we can call it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots = do_roots()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now how long did it take?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# a semi-colon suppresses output\n",
    "@time roots = do_roots();\n",
    "@time roots = do_roots();  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It allocates a million times less memory, and is 50 times faster! This is the first lesson about performance in Julia: *always* put everything in a function.\n",
    "\n",
    "Note that the first time we ran the function, it took longer. This is due to the fact that the first time a function is run with arguments of given types, the function is *compiled*. Subsequent runs with the same types of arguments reuse the previously-compiled code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Use a `while` loop with a suitable condition to improve the code for the Newton method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generic functions and methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our code currently is not very flexible. To make it more flexible, we would like to pass in arguments to the `do_roots` function. We can make a version which takes as arguments the functions `f` and `f'`, for example. Functions are \"first-class objects\" in Julia, so they can just be passed around by name.\n",
    "\n",
    "Let's redefine our function `do_roots` to accept these arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function do_roots(f, f′)\n",
    "    initial_conditions = -100:0.01:100\n",
    "    roots = similar(initial_conditions)\n",
    "\n",
    "    for (j, x_0) in enumerate(initial_conditions)\n",
    "        x = x_0\n",
    "\n",
    "        for i in 1:1000\n",
    "            x = x - f(x) / f′(x)\n",
    "        end\n",
    "\n",
    "        roots[j] = x\n",
    "    end\n",
    "    \n",
    "    roots\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the output that Julia returns: \"generic function with 2 methods\". This is a sign that something interesting is happening. In fact, we have not \"redefined\" the function `do_roots`; rather, we have defined a *new version* of `do_roots`, which accepts a *different set of arguments*. (The collection *and types* of the arguments that a function accepts are called its **type signature**.)\n",
    "\n",
    "Indeed, the function `do_roots` now has *two different methods* or versions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(do_roots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we call `do_roots` with no arguments, the first version will be used; calling it with two arguments will call the second version. The process of choosing which \"version\" of a function to call is called *dispatch*. The fundamental fact in Julia is that (almost) all functions are such \"generic functions\" with multiple version, i.e. Julia is one of very few languages that use **multiple dispatch**. This turns out to be very natural for many applications in scientific computing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arguments `f` and `f'` in the second method of `do_roots` are names that are local to the function. We have functions of the same name defined globally, so we can pass those in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@time do_roots(f, f′);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is faster than the first version of `do_roots`, but much slower than the good version. It turns out that Julia currently *cannot optimize* (inline) functions passed in this way. This is something to bear in mind -- there is (currently) a trade-off between user convenience and speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia also has *anonymous* functions, which allow us to pass in a function that we define \"in the moment\", without giving it a name. For example, let's do the exercise with a more interesting function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@time roots = do_roots(x->(x-1)*(x-2)*(x-3), x->3x^2-12x+11);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that anonymous functions are currently very slow. However, there are workarounds, e.g. the `FastAnon` package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results for this function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "plot(-100:0.01:100, roots)\n",
    "xlim()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(6,4))\n",
    "plot(-100:0.01:100, roots)\n",
    "xlim(1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complexifying Newton "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous result is still pretty boring. It turns out that the Newton method gets interesting if we look for roots of functions of *complex* numbers. [If you are not familiar with complex numbers, you can think of them as pairs of real numbers that have certain mathematical operations defined.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to use the Newton method starting from initial conditions distributed in the complex plane, i.e. pairs $a + bi$, where $i = \\sqrt{-1}$. First of all let's see how Julia handles complex numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sqrt(-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh dear, that didn't work very well. It turns out that Julia is carefully designed to respect, when possible, the type of the input argument. Indeed, let's ask Julia what it thinks `sqrt` means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sqrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `sqrt` is a generic function, with the following methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(sqrt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia gives us a list of the available methods, together with links direct to the source code on GitHub (in IJulia) or locally (in Juno).\n",
    "\n",
    "`sqrt()` acting on a `Float64` returns a `Float64` when it can, or throws a `DomainError` when its argument is negative. To get square roots in the complex plane, we must *start* with a complex number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The names of types in Julia start with capital letters, so let's try `Complex`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Complex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we will see later, types have functions with the same name that act as **constructors** to make objects of the type. Let's see the available functions with the name `Complex`. Note that output has changed rather a lot between Julia v0.3 and Julia v0.4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(Complex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try playing with `Complex`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = Complex(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "b = Complex(3, 4.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `Complex` is also parametrised by the type of its real and imaginary parts.\n",
    "\n",
    "We can also make complex numbers directly using `im`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "3.0 + 4.0im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Here, 4.0im is multiplication of 4.0 by `im`, which represents $i$, the imaginary unit.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do complex arithmetic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is happening here? Julia knows how to do `*` for complex numbers. Let's ask Julia what `*` is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, mathematical operators *are generic functions too*! We can list all the ways to do `*`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of these are defined in Julia itself. (Although the definitions for basic types like `Int` are only shallow wrappers around underlying C code.) We see that generic functions can be a complicated \"patchwork\" made of different methods for different types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the exact method used for a given operation using `@which`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which a * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial conditions: matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to think about how to generate a grid of initial conditions of the form $a+bi$ in the complex plane, $\\mathbb{C}$. Firstly, we could just iterate over the initial conditions in two repeated `for`s, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in -2:1\n",
    "    for j in -2:1\n",
    "        println(\"($i, $j)\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have used **string interpolation**: the *value* of the variable `i` is substituted into the string instead of the sequence `$i$`. [Note that this is not recommended for performance-critical applications.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we still require somewhere to store the results. It is natural to use a **matrix**. A simple way of generating a matrix is the `zeros` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeros(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that with a single element, we generate a *vector* of zeros, while"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeros(3, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "gives a *matrix*, i.e. a 2-dimensional `Array`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiple dispatch allows Julia to provide convenience versions of functions like this. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeros(-3:2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "creates a vector of the same length as the range!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this does not work for two different ranges:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zeros(-3:2, -3:2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `length` for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "linear_initial_conditions = -5:0.1:5\n",
    "roots = zeros(length(linear_initial_conditions), length(linear_initial_conditions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if we try to store a complex number in this matrix, we find a problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots[1, 1] = 3+4im"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An `InexactError` is a sign that we are trying to put a value into a type that it \"doesn't fit into\", for example a `Float64` into an `Int64`, or, in this case, a complex number into a float. We must instead create the matrix to hold complex numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "linear_initial_conditions = -5:0.1:5\n",
    "roots = zeros(Complex128, length(linear_initial_conditions), length(linear_initial_conditions))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `Complex128` is just an alias (an alternative name) for `Complex{Float64}`, so called because two 64-bit `Float64`s require `128` bits of storage in total."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can insert complex values into the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots[1, 1] = 3+4im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing Newton for complex functions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to make a version of Newton for complex functions. We will try to find cube roots of $1$ in the complex plane, by finding zeros of the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(z) = z^3 - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f′(z) = 3z^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function do_complex_roots(range=-5:0.1:5)  # default value\n",
    "\n",
    "    L = length(range)\n",
    "    \n",
    "    roots = zeros(Complex128, L, L)\n",
    "\n",
    "    for (i, x) in enumerate(range)\n",
    "        for (j, y) in enumerate(range)\n",
    "            \n",
    "            z = x + y*im\n",
    "            \n",
    "            for k in 1:1000\n",
    "                z = z - f(z) / f′(z)\n",
    "            end\n",
    "\n",
    "            roots[i,j] = z\n",
    "        end\n",
    "        \n",
    "    end\n",
    "    \n",
    "    roots\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots = do_complex_roots(-5:0.1:5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use `PyPlot` to plot the result. `PyPlot` only understands floating-point matrices, so we'll take the imaginary part:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow(imag(roots))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia uses \"column-major\" storage, whereas Python uses \"row-major\", so in fact we need to flip $x$ and $y$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function do_complex_roots(range=-5:0.1:5)  # default value\n",
    "\n",
    "    L = length(range)\n",
    "    \n",
    "    roots = zeros(Complex128, L, L)\n",
    "\n",
    "    for (i, x) in enumerate(range)\n",
    "        for (j, y) in enumerate(range)\n",
    "            \n",
    "            z = y + x*im\n",
    "            \n",
    "            for k in 1:1000\n",
    "                z = z - f(z) / f′(z)\n",
    "            end\n",
    "\n",
    "            roots[i,j] = z\n",
    "        end\n",
    "        \n",
    "    end\n",
    "    \n",
    "    roots\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow(imag(do_complex_roots(-3:0.01:3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Array comprehensions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia has a neat syntax for constructing arrays from iterables that is very similar to mathematical notation.\n",
    "For example, the squares of the numbers from 1 to 10 is\n",
    "\n",
    "$$\\{x^2: x \\in \\{1,\\ldots,10\\} \\},$$\n",
    "\n",
    "i.e. \"the set of $x^2$ for $x$ from $1$ to $10$. In Julia we can write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "squares = [x^2 for x in 1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a Newton function by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function newton(x0, N=100)\n",
    "    x = x0\n",
    "    \n",
    "    for i in 1:N\n",
    "        x = x - f(x) / f′(x)\n",
    "    end\n",
    "    \n",
    "    x\n",
    "end\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then our Newton fractal can be written very concisely as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "methods(newton)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the effect of a default argument is simply to create an additional method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function newton_fractal(range)\n",
    "    [newton(b+a*im) for a in range, b in range]\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add labels using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "?text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow(imag(newton_fractal(-3:0.01:3)), extent=(-3, 3, -3, 3))\n",
    "\n",
    "text(1, 0, L\"1\")\n",
    "text(reim(exp(2π*im/3))..., L\"e^{2\\frac{\\pi}{3}}\")\n",
    "text(reim(exp(-2π*im/3))..., L\"e^{-2\\frac{\\pi}{3}}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have used Julia's `reim` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reim(exp(2π*im/3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It returns a tuple. The `...`, or *splat*, operator, unpacks the tuple into two arguments.\n",
    "The `L\"...\"` notation is a special string macro available in the `LateXStrings` package used by `PyPlot`,\n",
    "that makes a LaTeX string."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise: Make a version that accepts functions and experiment with other complex polynomials. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introspection and iteration protocol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does Julia know how to iterate using `for` through a vector or range? Let's look at a Unicode string:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = \"aαbβ\"  # use `\\alpha<TAB>`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia provides access to several layers between the high-level code we write and the low-level machine code that is finally produced by the compilation process. The first of those is a \"lowered\" version of the code, in which high-level syntax is transformed to Julia code at a lower level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@code_lowered iterate(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there are three important functions: `start`, `next` and `done`.\n",
    "For example, iterating through a Unicode `UTF8String` is complicated, since characters have different lengths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonetheless, we can iterate through `s`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "function string_iterate(s)\n",
    "    for c in s\n",
    "        println(c)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "string_iterate(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, we can extract a list of the characters in `s` with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "chars = [c for c in s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "chars[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typeof(chars[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in Julia, strings are written with `\"` and characters with `'` (as in C)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interface that allows us to iterate over an object using `for` is provided by three functions `start`, `next` and `done` that must be defined for that type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "start(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "next(s, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "next(s, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "done(s, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which start(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "For more details about introspection, check out Leah Hanson's [blog post](http://blog.leahhanson.us/julia-introspects.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@which next(s, 1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.10-pre",
   "language": "julia",
   "name": "julia-0.3"
  },
  "language_info": {
   "name": "julia",
   "version": "0.3.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}