08_BestPractice_Optimisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimisation\n",
    "\n",
    "Often the code we initiially write to solve a problem is primarily focused on accuracy but often once this is achieved we would like to find ways to make the code run faster.  Python isn't the best code to do this but we will we will use it to investigate the process for optimising code.  We will revisit this when we learn cython.\n",
    "\n",
    "In order to understand how to make code faster it can be useful to know the basics of how a computer works. Most computers have an architecture that looks a bit like this:\n",
    "\n",
    "![](Plots/Computer_architecture.png)\n",
    "\n",
    "Modern CPU's are so fast that actually most optimisatising for high performance computing comes from two approaches: 1. optimising memory usage to make sure the CPU is fed with enough data 2. utilising parallelisation which we will meet next.  There are three numbers listed for each section in the above plot: Latency, which is how long it takes to respond to a command; Bandwidth, which is how fast it can pass data up the chain; Size, which is how much data can be stored there.  The main problem we see is that memory latency is ~200 times longer than CPU latency (which is 1 cycle).  So, if you want to multiply two numbers it takes 200 times as long to get them to the CPU as it does to multiply them.  The response to this was to introduce multiple layers of faster and faster, but smaller and smaller, memory inbetween the CPU and main memory called 'cache'.  This helps keep data that we are working on close to the CPU so if we want to use it we can without delay.  This is a problem that is getting worse with time as we can see in the following:\n",
    "\n",
    "![](Plots/cpuvsmem.png)\n",
    "\n",
    "This means that there are two things we can do to speed up our code.\n",
    "\n",
    "- Keep data local, ie near the CPU, and \n",
    "- reuse it as much as possible while it is there.\n",
    "\n",
    "Obviously you have no direct control over the cache but there are some things you can do that help the computer optimise this.\n",
    "\n",
    "The best analogy for a computer is trying to learn something from textbooks.  Now the heirarchy looks like this:\n",
    "\n",
    "![](Plots/Library.png)\n",
    "\n",
    "And it is easy to get some basic lessons\n",
    "\n",
    "- Get all the books you need from the library at once because it's a pain to keep going back there\n",
    "- Constantly swapping books is annoying so try to finish working on the open ones before you look at the next set\n",
    "\n",
    "With this in mind lets look at what we can do to optimise our code.\n",
    "\n",
    "Firstly, in python we don't actually have much control over what the computer does at the instruction level.  This is a shame as often it is possible to speed up some calculations by orders of magnitude just by using simple directives to compilers in lower level languages like C and Fortran.  This is why it is so important to use packages in python where this type of optimisation should already have been done for you.  Another key point is that you should try to use packages what have been build for your CPU (check this with `lscpu`), for example making sure you install the intel distribution for intel chips with:\n",
    "\n",
    "`conda create -n idp intelpython3_full python=3`\n",
    "\n",
    "### Timing operations\n",
    "\n",
    "There are always many different ways to do the same calculation so you should get into the habit of checking which is fastest.  For simple commands and functions we can just use `%timeit` (single line) and `%%timeit` (entire cell) which are useful for comparing equvilent code.  The best bit is they run the code multiple times to give you an accurate estimate of time.  Here are some examples of simple small speedups you can do:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Squaring vs multiplication vs adding\n",
    "x = 0.5\n",
    "%timeit x**2\n",
    "%timeit x*x\n",
    "%timeit x/x\n",
    "\n",
    "%timeit x+x\n",
    "%timeit 2.0*x\n",
    "%timeit 2.0*5.0\n",
    "%timeit x/0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Strictly speaking on a CPU addition is faster than multiplication which is faster than division.  Exponentiation is not a basic operation for the CPU and instead calls a function which takes the log then multiplies.  However, in python we see that while exponentiation is slower than multiplying, multipling is the same as adding.  Multiplying by a number is faster than multiplying by a variable as we don't have to go and find the number it in memory and multipling two numbers much faster again.  We can roughly conclude that it takes twice as long to find a variable as it does to multiply it (you can read faster than find the right bit of a book). In python division and muliplication take the same time.  This is because most of the time is lost in finding the variables rather than the operation itself.  But if we do the operation alot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit [0.5*x for x in range(10000)]\n",
    "%timeit [x/2.0 for x in range(10000)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then small differences can emerge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Square root\n",
    "import math\n",
    "%timeit x**0.5\n",
    "%timeit math.sqrt(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python functions can be slower than basic operations.   I think here `sqrt` most likely has been tuned for higher accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 3a. Math vs Numpy\n",
    "import math\n",
    "import numpy as np\n",
    "%timeit math.sin(x)\n",
    "%timeit np.sin(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So numpy is slower for scalars but..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3b. Math vs Numpy for vector\n",
    "X = np.random.random(100)\n",
    "Y1 = np.zeros(100)\n",
    "Y2 = np.zeros(100)\n",
    "Y3 = np.zeros(100)\n",
    "%timeit Y1 = np.sin(X)\n",
    "%timeit Y2 = list(map(lambda x: math.sin(x),X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "for i in range(100):\n",
    "    Y3[i] = math.sin(X[i]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "faster for arrays, note the tiny scaling which shows that almost all the time is getting the function itself so again most of the time is in finding things, not calculation.  Here we have used `map` which is a way of applying a function to a list with the format:\n",
    "\n",
    "map(function,list) ie:  `list(map(lambda x: x**2, items)`\n",
    "\n",
    "There is also `filter` which selects based on a condition:\n",
    "\n",
    "filter(function,list) ie: `list(filter(lambda x: x < 0, items))`\n",
    "\n",
    "These are faster than loops but can't compete with numpy.   Here is something similar for matrix multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4a. Math vs Numpy multiplication\n",
    "x=2.3\n",
    "y=3.4\n",
    "%timeit np.dot(x,y)\n",
    "%timeit x*y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 4b. Math vs Numpy for dot product\n",
    "X = np.random.random(100)\n",
    "Y = np.random.random(100)\n",
    "%timeit np.dot(X,Y)\n",
    "%timeit sum(list(map(lambda x:x[0]*x[1],list(zip(X,Y)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the call time for numpy is all the time and much faster than `map` construct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5. Constructing lists\n",
    "%timeit [x*x for x in range(100)]\n",
    "%timeit [x**2 for x in range(100)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "list1=[]\n",
    "for x in range(100):\n",
    "    list1.append(x*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "list1=[]\n",
    "append = list1.append\n",
    "for x in range(100):\n",
    "    append(x*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "list1 = np.zeros((100))\n",
    "for x in range(100):\n",
    "    list1[x] = x*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "np.fromfunction(lambda i,:i*i,(100,),dtype=int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So list comprehensions are fastest. Also note that using `append = list1.append` saves us about 20% in time.  This is due to not having to looking up the function, the assignment has make it more 'local' in memory.  Do be careful with this as it can make code very hard to read so you should generally aviod it.  Also note that subsequent operations on list1 may be quicker for the last 2 methods as the list is more likely to be contiguous in memory (all in the same place) for very large lists so we can read it to cache faster.   Finally, note that using a list comprehension won't make you code fast if you use bad algorithms like `x**2` rather than `x*x` which costs ~5 times as much.  In general better algorithms will always win over better code.  This is why prototyping is such a good idea."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#7. Different indexing orders\n",
    "import numpy as np\n",
    "A = np.zeros((100,100))\n",
    "B = np.random.random((100,100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "for i in range(100):\n",
    "    for j in range(100):\n",
    "        A[i,j] = B[i,j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "for i in range(100):\n",
    "    for j in range(100):\n",
    "        A[j,i] = B[j,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In C or Fortran this would have mattered alot, but here the loops are so slow you don't see any difference due to how you are accessing memory.  This could still be something to try for very large arrays as it may eventually win.  The problem is that python stores arrays as a single line of all the rows in order so the `(i,j)` loop reads the data in order, so the next value is just next to the last one.  The `(j,i)` has to jump the length of the row to find the next value slowing it down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#8. Move stuff out of loops if possible\n",
    "import math\n",
    "\n",
    "def func1():\n",
    "    x=math.sqrt(2)\n",
    "    y=0\n",
    "    for i in range(100):\n",
    "        y*=x\n",
    "        \n",
    "def func2():\n",
    "    y=0\n",
    "    for i in range(100):\n",
    "        y*=math.sqrt(2)\n",
    "               \n",
    "def func3():\n",
    "    y=0\n",
    "    for i in range(100):\n",
    "        y*=1.414213\n",
    "        \n",
    "%timeit func1()\n",
    "%timeit func2()\n",
    "%timeit func3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no point in working out the square root of 2 100 times. Again try to make stuff used in loops 'local' and avoid any calculation if possible (this doesn't happen for for loops which only expand the `for` line once, `while` does it every iteration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#8a. Move stuff out of loops if possible\n",
    "import math\n",
    "list1 = [range(100)]\n",
    "\n",
    "def func1():\n",
    "    y=0\n",
    "    i=0\n",
    "    while i<len(list1):\n",
    "        y+=i\n",
    "        i+=1\n",
    "        \n",
    "def func2():\n",
    "    y=0\n",
    "    i=0\n",
    "    length = len(list1)\n",
    "    while i<length:\n",
    "        y+=i\n",
    "        i+=1\n",
    "        \n",
    "%timeit func1()\n",
    "%timeit func2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same goes for test conditons for while loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#9. imports inside Vs outside of funtions\n",
    "\n",
    "import math\n",
    "def func1():\n",
    "    math.sin(0.3)\n",
    "    \n",
    "def func2():\n",
    "    import math\n",
    "    math.sin(0.3)\n",
    "    \n",
    "    \n",
    "%timeit func1()\n",
    "%timeit func2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be tempting to use imports inside functions so it looks tidy but this has significant cost. Best leave them at the top of the module as functions may be used multiple times but modules are (generally) only loaded once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#10. Function overhead\n",
    "\n",
    "def func1(i):\n",
    "    return i*i\n",
    "\n",
    "def func2():\n",
    "    x = 0e0\n",
    "    for i in range(100):\n",
    "        x += func1(i)\n",
    "        \n",
    "def func3():\n",
    "    x = 0e0\n",
    "    for i in range(100):\n",
    "        x += i*i\n",
    "        \n",
    "%timeit func2()\n",
    "%timeit func3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modulisation is great but function calls can be expensive so don't go crazy (especially in loops!).  Again this is introducing an extra lookup which makes the data less likely to be local in cache"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#11. Multiple assignment rather than tempory variables\n",
    "\n",
    "def func1():\n",
    "    a=0\n",
    "    b=1\n",
    "    for i in range(1000):\n",
    "        a,b = b,a+b\n",
    "        \n",
    "def func2():\n",
    "    a=0\n",
    "    b=1\n",
    "    for i in range(1000):\n",
    "        c = a+b\n",
    "        a = b\n",
    "        b = c\n",
    "        \n",
    "%timeit func1()\n",
    "%timeit func2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A small saving (plus looks better)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 12. Finding elements depends on the data type\n",
    "list1 = [range(100)]\n",
    "set1 = set(list1)\n",
    "\n",
    "%timeit 5 in list1\n",
    "%timeit 5 in set1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sets and dictionaries are hash tables so are faster to search than lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 13. Read entire files at once if possible\n",
    "import numpy as np\n",
    "def read1():\n",
    "    data1 = np.loadtxt('Data/Period1.txt')\n",
    "    return data1\n",
    "\n",
    "def read2():\n",
    "    file = open('Data/Period1.txt','r')\n",
    "    data1 = []\n",
    "    for line in file:\n",
    "        tmp = line.split()\n",
    "        data1.append(float(tmp[0]))\n",
    "    data2 = np.array(data1)\n",
    "    return data2\n",
    "\n",
    "def read3():\n",
    "    file = open('Data/Period1.txt','r')\n",
    "    data1 = file.read()\n",
    "    data1 = data1.split('\\n')\n",
    "    data1 = [float(x) for x in data1[:-1]]\n",
    "    data1 = np.array(data1)\n",
    "    return data1\n",
    "\n",
    "%timeit read1()\n",
    "%timeit read2()\n",
    "%timeit read3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So numpy is the cleanest to code but the slowest to run and reading the file in one block is better than line by line (remember getting all library books at once rather than making multiple trips).  This is especially true on shared systems where you may compete for bandwidth for I/O.  There the system may alternate read statements from competing programmes eg: someone reads 1Tb, you read a line, someone reads 1Tb, you read one line etc.. which can be crippling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above will only give you small gains but can be worth considering.  Real gains almost always come finding better methods as:\n",
    "\n",
    "- Algorithm improvement virtually always beats optimisation\n",
    "\n",
    "If your code is slow you should first think about \"how else I could try to calculate this?\" rather than trying to find faster ways to do the same approach.\n",
    "Don't be tempted do do any of these if it sacrifices code readability.  There is a famous quote amongst computer scientsts:\n",
    "\n",
    "<b>\"Premature optimisation is the root of all evil\"</b>\n",
    "\n",
    "Your first goal should always be clarity (luckly in python clarity often coincides with speed).\n",
    "\n",
    "## Algorithmic Complexity\n",
    "\n",
    "When selecting algorithms you should think about how they scale with the size of the problem.  This is the computational or algorithimic complexity.  There are seveal main groups: `linear`, `log-linear`, `polynomial` and `exponential`. Unlike the above where simplicity makes the code fast, usually the simplest algorithm will scale the worst.  A common example is sorting.  The simplest algorithim is `selection sort` which works like this:\n",
    "\n",
    "1. Iteritivly find the ith smallest element and swap it's positon with the element at postion i\n",
    "\n",
    "This is fairly simple and the way most people would think about sorting.  It is an $O(n^2)$ algorithm as we have to find the smallest of the remaining objects, which is $O(n)$, $O(n)$ times.  Code to do it is fairly simple:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def selection_sort(list1):        \n",
    "    for i in range(len(list1)):\n",
    "        minimum = i\n",
    "        \n",
    "        for j in range(i + 1, len(list1)):\n",
    "            # Select the smallest value\n",
    "            if list1[j] < list1[minimum]:\n",
    "                minimum = j\n",
    "\n",
    "        # Switch the element with the one at position i \n",
    "        list1[minimum], list1[i] = list1[i], list1[minimum]\n",
    "            \n",
    "    return list1\n",
    "\n",
    "\n",
    "        \n",
    "list1 = [7,5,7,9,5,3,5,4,8]  \n",
    "\n",
    "print(selection_sort(list1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However much faster algorithms exist which scale much better.  If we use a divide an conquer approach this complexity can be reduced to $O(nlog(n))$.  This can be seen in the example of `merge sort` which works using the fact that merging sorted lists is $O(n)$.  Our algorithm looks like this:\n",
    "\n",
    "1. Make $n$ lists 1 element long.\n",
    "2. Merge list in pairs to make $n/2$ sorted lists, 2 elements long\n",
    "3. Itererate until only one sorted list remains.\n",
    "\n",
    "Here is the code to do it, it's much more complicated now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_sort(arr):\n",
    "    # The last array split\n",
    "    if len(arr) <= 1:\n",
    "        return arr\n",
    "    mid = len(arr) // 2\n",
    "    # Perform merge_sort recursively on both halves\n",
    "    left, right = merge_sort(arr[:mid]), merge_sort(arr[mid:])\n",
    "\n",
    "    # Merge each side together\n",
    "    return merge(left, right, arr.copy())\n",
    "\n",
    "\n",
    "def merge(left, right, merged):\n",
    "\n",
    "    left_cursor, right_cursor = 0, 0\n",
    "    while left_cursor < len(left) and right_cursor < len(right):\n",
    "      \n",
    "        # Sort each one and place into the result\n",
    "        if left[left_cursor] <= right[right_cursor]:\n",
    "            merged[left_cursor+right_cursor]=left[left_cursor]\n",
    "            left_cursor += 1\n",
    "        else:\n",
    "            merged[left_cursor + right_cursor] = right[right_cursor]\n",
    "            right_cursor += 1\n",
    "            \n",
    "    for left_cursor in range(left_cursor, len(left)):\n",
    "        merged[left_cursor + right_cursor] = left[left_cursor]\n",
    "        \n",
    "    for right_cursor in range(right_cursor, len(right)):\n",
    "        merged[left_cursor + right_cursor] = right[right_cursor]\n",
    "\n",
    "    return merged\n",
    "\n",
    "list1 = [7,5,7,9,5,3,5,4,8]  \n",
    "\n",
    "print(merge_sort(list1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "list1 = np.random.normal(0,5,(10))\n",
    "list2 = np.random.normal(0,5,(10))\n",
    "\n",
    "%timeit selection_sort(list1)\n",
    "%timeit merge_sort(list2)\n",
    "\n",
    "list1 = np.random.normal(0,5,(1000))\n",
    "list2 = np.random.normal(0,5,(1000))\n",
    "\n",
    "%timeit selection_sort(list1)\n",
    "%timeit merge_sort(list2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This type of optimisation is tricky and requires a bit of clever thinking.  For standard operations it is work serching a bit to check if there is a fancy way of doing it (and a package that implements it!) but for specilised cases you are best to try asking your local computer officers or computer science graduates.  \n",
    "\n",
    "The key to optimising your code is to note that usually in any code it's only 10% (or less!) of it that takes 90% of the time. So optimising 90% of the code will have almost no effect on runtime. This means we should initially always write for clarity then try to find which parts ate taking most of the time.\n",
    "\n",
    "How do we find the 10% of our code to focus on?  The answer is to **profile** your code.\n",
    "\n",
    "\n",
    "### Profiliing\n",
    "\n",
    "Profilers analyse your code and tell you what parts are taking the most time (and memory) to run.  There are a few you can use.  Let's look at them with our little prime number generator we created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def primenum(N):\n",
    "    primes = [2]\n",
    "    for n in range(3,N):\n",
    "        if all(n%p>0 for p in primes):\n",
    "            primes.append(n)\n",
    "    return primes \n",
    "\n",
    "%prun primenum(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This indicates that we spent most of our time in `<genexp>` on line 4 with the next most time spent on the built-in method `all`.  Which is what we would expect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we could install the `line_profiler` package (`conda install line_profiler`) and use this instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext line_profiler\n",
    "\n",
    "def primenum(N):\n",
    "    primes = [2]\n",
    "    for n in range(3,N):\n",
    "        if all(n%p>0 for p in primes):\n",
    "            primes.append(n)\n",
    "    return primes \n",
    "\n",
    "\n",
    "%lprun -f primenum primenum(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is a little easier to read. Here is an intersting case to see how memory reuse is important"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%load_ext line_profiler\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "def nbody(particles):\n",
    "    for item1 in particles:\n",
    "        acc1 = 0\n",
    "        acc2 = 0\n",
    "        acc3 = 0\n",
    "        for item2 in particles:\n",
    "            dx = item1[0] - item2[0]\n",
    "            dy = item1[1] - item2[1]\n",
    "            dz = item1[2] - item2[2]\n",
    "            \n",
    "            distsquare = dx*dx+dy*dy+dz*dz\n",
    "            if abs(distsquare) > 1e-9:\n",
    "                distinv = 1e0/math.sqrt(distsquare)\n",
    "                acc1 += dx * distinv*distinv*distinv\n",
    "                acc2 += dy * distinv*distinv*distinv\n",
    "                acc3 += dz * distinv*distinv*distinv\n",
    "\n",
    "particles = np.random.random((100,3))\n",
    "\n",
    "%lprun -f nbody nbody(particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So why do the 2nd and 3rd lines get faster for the blocks of three commands next to each other?  It's because the first moves the data to l1 cache so it's easier to use in the next line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one of our main problems is likely to be memory useage we can also profile memory with the `memory_profiler` package (`conda install memory_profiler`) but this only works on files so we need to save this to a file first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file Tools/primetools.py\n",
    "def primenum(N):\n",
    "    primes = [2]\n",
    "    for n in range(3,N):\n",
    "        if all(n%p>0 for p in primes):\n",
    "            primes.append(n)\n",
    "    return primes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%load_ext memory_profiler\n",
    "import sys\n",
    "sys.path.append('./Tools') \n",
    "import primetools as pts \n",
    "%mprun -f pts.primenum pts.primenum(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use `%run -p` to profile scripts that we run.  Finally there is also another build in profiler `cProfile`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cProfile\n",
    "def primenum(N):\n",
    "    primes = [2]\n",
    "    for n in range(3,N):\n",
    "        if all(n%p>0 for p in primes):\n",
    "            primes.append(n)\n",
    "    return primes \n",
    "\n",
    "cProfile.run('primenum(10000)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can also be run from the command line with `python -m cProfile [-o output_file] [-s sort_order] myscript.py` which is usefull\n",
    "\n",
    "### Exercise:\n",
    "\n",
    "Profile both your periodic signal code from lecture 6 and your 'game of life' code from 7 and see if you can get them to run any faster. Try the same with you previous solutions from earlier lectures.  You should use this in your protoyping stage to check the quality of modules you may want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}