9_Cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cython\n",
    "\n",
    "One of the best ways to sppe up python code is to convert it into compiled C code. Luckly this is fairly easy we can actually do it in a jupyter notebook which we can use for testing things (to use it in scripts we will need some extra steps).  First we install cython using `conda install cython` then we need to load the extension to the notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare a normal python function and it's cython-ized version.  We will use the example if a function that calculates the Nth Fibonacci number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fib1(N):\n",
    "    a,b = 0,1\n",
    "    for i in range(N):\n",
    "        a,b = b,a+b\n",
    "    return a\n",
    "\n",
    "%timeit fib1(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def fib2(N):\n",
    "    a,b = 0,1\n",
    "    for i in range(N):\n",
    "        a,b = b,a+b\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit fib2(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def fib3(int N):\n",
    "    cdef int i\n",
    "    cdef int a=0,b=1\n",
    "    for i in range(N):\n",
    "        a,b = b,a+b\n",
    "    return a\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit fib3(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So just adding `%%cython` gives us a factor ~2x speedup.  But if we simply add types to our variables with `cdef` this increases to a ~240x speed up!\n",
    "\n",
    "This is because the function is dominated by the loop which C can do much better.  The `%%cython` magic actually does something tricky in the background. It takes the cell and converts it to C code then compiles it and stores the resulting executable in a temporary location.  We can see the actuall C code generated using the annotate option by adding `-a` after the `%%cython`.  This gives us a window to how the code has been converted to C with highlighting to show how much python interaction is left for each line.   If we click the little '+' on the line number it shows you what this line has been converted to in C and the stronger the yellow the more python interation remains.\n",
    "\n",
    "We will come back to compilation later but let's look at the difference between writing cython and python code:\n",
    "\n",
    "1. We don't have to do anything to cython-ise most python code.  We can put almost any python code through the cython compiler and it will work fine and usually run faster.\n",
    "\n",
    "2. To access performance of C with Cython we usually only have to declare types using `cdef` and sometimes switch the default behaviour of some operations using simple flags.\n",
    "\n",
    "3. In cython we can now use all C libraries and easily access threaded parallelisim by avoiding the GIL.\n",
    "\n",
    "So we see that there are very few differences.  Cython is a superset of python so we don't have to change anything if we don't want to.  As cython is effectivly an optimisation tool we should profile the code and only cythonise the slowest parts.  This is the main advantage. If you wanted to access the speed of C you would otherwise have to re-write all your code in C where lots of things can be significantly more difficult.  Instead we can use the convenience of python for most of the code and only invoke C in the sections where performance is most important.\n",
    "\n",
    "## Types\n",
    "Using cython is it's basic form is pretty easy.  Let's look at the cdef statement a bit more. Here are the following basic cdef types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "cdef char i=1           # Oddly an 8 bit integer (-128 to 127) (it's enough to label all charcaters so can be used for strings)\n",
    "cdef short j=2          # 16 bit integer (-32,768 to 32,767)\n",
    "cdef int k=3            # 32 bit integer (-2,147,483,648 to 2,147,483,647)\n",
    "cdef unsigned int l=4   # 32 bit +ve integer (0 to 4,294,967,295), \"unsigned\" can go infront of all numeric types\n",
    "cdef long int m=5       # 64 bit integer (-9,223,372,036,854,775,808 to 9,223,372,036,854,775,807)\n",
    "cdef float x=0.0        # 32 bit float (6 decimal places, max exponent 38)\n",
    "cdef double y = 0e0     # 64 bit float (12 decimal places, max exponent 1023)\n",
    "cdef list list1 = [1,2,3]       # just a normal list (not much performance gain)\n",
    "cdef dict dict1 = {'a':1,'b':2} # just a normal dict (not much performance gain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here once we define a type we have to stick with it unlike python which dynamically changes types to accuratly store any number you give it.  This means we are now in danger of overflow errors.  This is when you assign `cdef short j` then write `j = 200**2` and get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "cdef short j\n",
    "j = 200**2\n",
    "print(j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is because 40,000 is larger than 32,767 so we wrap around to the negative part.  Similarly if we try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "cdef unsigned int j\n",
    "j = -1\n",
    "print(j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so we have to be a bit carefull with our variables to avoid strange results.\n",
    "\n",
    "Strings are stored completly differently in C so there is no `cdef` just for them.  Instead they are just a array of `char`.  The `char*` means that it is an address to the point in memory where the string begins. Also python and C encode strings differently so you have to `encode` and `decode` for them to be able to talk to each other.  It's best just to keep strings as python variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def test(input):\n",
    "    input_byte = input.encode('utf-8')\n",
    "    cdef char* c_string = input_byte\n",
    "    cdef bytes py_string_byte = c_string\n",
    "    output = py_string_byte.decode('utf-8')\n",
    "    print(output)\n",
    "test('Hello')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " We can also use any of the standard C math libraries with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "from libc.math cimport sin\n",
    "\n",
    "def sin_c(double x):\n",
    "    return sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "x = 0.5\n",
    "%timeit math.sin(x)\n",
    "%timeit sin_c(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which are a bit faster.  To use numpy arrays we have to do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "cimport numpy as cnp\n",
    "\n",
    "cdef cnp.ndarray array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but this won't give us all the speed improvement possible as C doesn't know how to allocate the memory as it doesn't know the shape and datatype of the array.  Instead it is better to do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "cimport numpy as cnp\n",
    "\n",
    "def matrix_dot(cnp.ndarray[cnp.int_t, ndim=2] array1, cnp.ndarray[cnp.int_t, ndim=2] array2):\n",
    "    cdef cnp.ndarray[cnp.int_t, ndim=2] array3\n",
    "    array3 = np.dot(array1,array2)\n",
    "    return array3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "array1 = np.ones((100,100),dtype=np.int)\n",
    "array2 = np.ones((100,100),dtype=np.int)\n",
    "%timeit array3 = matrix_dot(array1,array2)\n",
    "%timeit array4 = np.dot(array1,array2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: we had to put this declaration in a function, this is so cython knows how long the memory needs to be allocated for as it's local to the function. Also you can `cimport numpy as np` I did it to a different name so you could see which is doing what.  Also numpy is already in C so as expected wraping it in cython doesn't help.\n",
    "\n",
    "We can also optimise the function call by specifying the return type.  For functions we have three choices: `def`, `cdef` and `cpdef`.  The first says it's callable in python or cython, the second cython only with optimised call, the third is callable in python and cython but optimised in the second case. If you use cdef or cpdef you need to add the type for the return variable like below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "cimport numpy as cnp\n",
    "\n",
    "cpdef cnp.ndarray[cnp.int_t, ndim=2] matrix_dot2(cnp.ndarray[cnp.int_t, ndim=2] array1, cnp.ndarray[cnp.int_t, ndim=2] array2):\n",
    "    cdef cnp.ndarray[cnp.int_t, ndim=2] array3\n",
    "    array3 = np.dot(array1,array2)\n",
    "    return array3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython for scripts\n",
    "\n",
    "So using the `%%cython` magic is pretty cool but we can't write a code using it.  So how do we use cython in our normal python code?  It's a four step process (two more than normal):\n",
    "\n",
    "1. Put your cython code in a file with extension `.pyx` like `cython_module.pyx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file cython_module.pyx\n",
    "\"\"\"\n",
    "Cython code for fibonnaci numbers\n",
    "\"\"\"\n",
    "cpdef int fibonacci(int N):\n",
    "    cdef int i\n",
    "    cdef int a=0,b=1\n",
    "    for i in range(N):\n",
    "        a,b = b,a+b\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Create a file called setup.py with the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file setup.py\n",
    "from distutils.core import setup\n",
    "from Cython.Build import cythonize\n",
    "\n",
    "setup(\n",
    "    ext_modules = cythonize(\"cython_module.pyx\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Now compile the code on the command line with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "python3 setup.py build_ext --inplace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Use the new cython functions with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cython_module as cym\n",
    "cym.fibonacci(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are now free to use the functions in `cython_modules` in python.\n",
    "\n",
    "**Example**\n",
    "Do this with other simple functions we have created.\n",
    "\n",
    "\n",
    "If we look in the directory we see two new files `cython_modules.c` and `cython_modules.so`  The `.c` is the transliteration of our cython code into C and the `.so` file is the compiled version of it.  If you open the `.c` file you will see that it is now about 2600 lines long.  Mostly it's definitions with the actual calculation appearing around line 1070 and lasting about 80 lines.  It is clear from the `.c` code that the code is doing a lot of checks which python does in the background which can slow down operation of the code.  Again we can see how well we are doing by using the annotate option in our `setup.py` file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file setup.py\n",
    "from distutils.core import setup\n",
    "from Cython.Build import cythonize\n",
    "\n",
    "setup(\n",
    "    ext_modules = cythonize(\"cython_module.pyx\", annotate=True)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can build it again but we will have to delete the `.so` files to make it run (otherwise it doesn't think anything has changed). This generates a `.html` file which shows us how much of our code has been converted to C.  It should have highlighted two lines the `def fib3()` line and the `return a` line.  This is because we haven't specified what type the function should return.  We can correct this by changing the definition to: `cpdef int fib3()`.  Now when we re-compile the `return` line is white and the `def` line is paler yellow.  This can't be changed as we want the function to be available in python so it must interact with it.\n",
    "\n",
    "## Extensions\n",
    "Now we have access to all of the functionallity of C and C++.  This is a massive topic and I couldn't begin to address it here.  There are however a couple of options I will flag up for you to think about in future\n",
    "\n",
    "Here is a link to compiler directives that can be specified in the setup file for all code or using decorators (which we haven't discussed but are just lines above a function begining with an @) for specific functions:\n",
    "https://cython.readthedocs.io/en/latest/src/userguide/source_files_and_compilation.html#compiler-directives\n",
    "Some common decorators are:\n",
    "- @cython.boundscheck(False)  Remove checks that you are accessing valid array entries\n",
    "- @cython.wraparound(False)   Remove the ability to use negative indexing in arrays\n",
    "- @cython.cdivision(False)    Use C's version of division rather than pythons so no more divide by zero errors\n",
    "\n",
    "- @cython.profile(True)  This is nessecary if you wan to profile using cProfile\n",
    "\n",
    "Turning these off and on can help you access more of the C speed by removing python style checks.  If you turn these off you code will usually just produce nonsense or explode when you do something wrong (like in C!) rather than raise an error (like in python).  These can buy some speed but are only really important if they are blocking a loop from being converted to C (where it would vectorise) or this particular loop contains only this type of calculation but this is quite hard to set up. You probably don't need to worry about them much\n",
    "\n",
    "The second is that now we can access task parallelisim both through the cython `prange` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "cimport numpy as cnp\n",
    "from cython.parallel import prange\n",
    "import cython\n",
    "\n",
    "@cython.cdivision(True)\n",
    "@cython.boundscheck(False)\n",
    "cpdef cnp.ndarray[cnp.int_t, ndim=1] func1(cnp.ndarray[cnp.double_t, ndim=1] Xin):\n",
    "    cdef int i\n",
    "    cdef int N = Xin.shape[0]\n",
    "    cdef cnp.ndarray[cnp.double_t, ndim=1] Xout = np.empty_like(Xin)\n",
    "    \n",
    "    for i in range(N):\n",
    "        Xout[i] = 1e0/Xin[i]\n",
    "        \n",
    "    return Xout\n",
    "\n",
    "@cython.cdivision(True)\n",
    "@cython.boundscheck(False)\n",
    "cpdef cnp.ndarray[cnp.int_t, ndim=1] func2(cnp.ndarray[cnp.double_t, ndim=1] Xin):\n",
    "    cdef int i\n",
    "    cdef int N = Xin.shape[0]\n",
    "    cdef cnp.ndarray[cnp.double_t, ndim=1] Xout = np.empty_like(Xin)\n",
    "    \n",
    "    for i in prange(N, nogil=True):\n",
    "        Xout[i] = 1e0/Xin[i]\n",
    "        \n",
    "    return Xout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "Xin = np.random.random((10000))+1e0\n",
    "\n",
    "%timeit func1(Xin)\n",
    "%timeit func2(Xin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can easily run into issues if you want to sum numbers as they all threads have to access the same variable. Cython does make sure the answer is right (unlike in C) but the code becomes effectively serial so the code will run slower due to the overheads for creating the threads in the first place. Still, this can be an easy way to paralleise simple loops.  Note that this will not happen if there is any python inside the loop, it has to be all cython."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython --compile-args=-fopenmp --link-args=-fopenmp\n",
    "from cython.parallel import prange\n",
    "\n",
    "cdef int i\n",
    "cdef int n = 30\n",
    "cdef int sum = 0\n",
    "\n",
    "# for i in range(n):\n",
    "for i in prange(n, nogil=True):\n",
    "    sum += i\n",
    "\n",
    "print(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrapping C with Cython\n",
    "\n",
    "ONe useful advanced topic is learning to wrap pure c code to be used in python.  This is fairly easy but a little fiddly for people new to C.  Here is a simple example to show you how it works.\n",
    "\n",
    "Suppose we have a function writtein in C that we want to use in python.  The C code is as follows:\n",
    "\n",
    "cexample.c :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <omp.h>\n",
    "#include \"cexample.h\"\n",
    "\n",
    "// The top two lines are the C version of loading modules\n",
    "// This is done via \"header\" files which list the functions to load.\n",
    "\n",
    "int fibonacci(int n){\n",
    "\t\n",
    "\tint i,a,b,tmp;\n",
    "\ta=0;\n",
    "\tb=1;\n",
    "\tfor (i=0;i<(n-1);i++){\n",
    "\t\ttmp = a+b;\n",
    "\t\ta = b;\n",
    "\t\tb = tmp;\n",
    "\t}\n",
    "\treturn b;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with header file cexample.h:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ifndef C_EXAMPLE_H\n",
    "#define C_EXAMPLE_H\n",
    "\n",
    "// The bit at the top is to check if the module has already been loaded\n",
    "// It asks if it is already defined and only loads it if not.\n",
    "\n",
    "int fibonacci(int n);\n",
    "\n",
    "#endif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must compile this code into a library which we will then load in cython.  To do this we run the two commands, the first compiles the `cexamples.c` into an object file.  The second creates the library **ar**chive from the object files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd ../Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "cd lib\n",
    "gcc -c cexample.c\n",
    "ar rcs libcexample.a cexample.o\n",
    "cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have created the C library we need to use cython to wrap it for use in python.  All we need to do is create the `.pyx` file for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file wrapcexamples.pyx\n",
    "cdef extern from \"cexample.h\":\n",
    "    int fibonacci(int n)\n",
    "\n",
    "def cfib(n):\n",
    "    cdef int m\n",
    "    m = fibonacci(n)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then create the `setup.py` to create the python module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file steup.py\n",
    "from distutils.core import setup\n",
    "from distutils.extension import Extension\n",
    "from Cython.Build import cythonize\n",
    "\n",
    "cexample_extension = Extension(\n",
    "    name=\"wrapcexample\",            \n",
    "    sources=[\"wrapcexample.pyx\"],\n",
    "    libraries=[\"cexample\"],\n",
    "    library_dirs=[\"lib\"],\n",
    "    include_dirs=[\"lib\"]\n",
    ")\n",
    "\n",
    "setup(\n",
    "    name=\"wrapcexample\",\n",
    "    ext_modules=cythonize([cexample_extension])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now if we run the build command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "python3 setup.py build_ext --inplace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can access the c function in pure python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import wrapcexample as cexmpl\n",
    "\n",
    "cexmpl.cfib(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a useful way to access thread level parallelisim as we can now use OpenMP (OMP) in native C where it's much easier.  An example is in the directory CythonOMP. The C code is below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include \"ompexample.h\"\n",
    "#include <omp.h>\n",
    "\n",
    "int pointlesssum(int n){\n",
    "\t\n",
    "\tint i,sum;\n",
    "\tsum = 0;\n",
    "\t\n",
    "\t#pragma omp parallel for default(none) private(i) shared(n) reduction(+:sum)\n",
    "\tfor(i=0;i<n;i++){\n",
    "\t\tsum += i;\n",
    "\t}\n",
    "\treturn sum;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now thread parallelisim is accessed via simple `pragmas` avoiding the trouble of having to turn off lots of python flags to get true parallelisim like for `prange`.  The above paralleises the `for` loop where `i` is private meaning each thread keeps it's own copy, `n` is public which means that all threads see the same variable and `sum` is a reduction variable meaning that each thread gets a local copy which are merged at the end of the loop by `+`. The `default(none)` means that parameters are neither private or shared by default, you have to specify them. Altenativly you can set the default to `private` then only specify variables that are shared or, more dangerously, set the default to shared then specify which are private. This is a pretty common way to use it.  You can create a parallel region with:|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pragma omp parallel default(none) private(...) shared(...)\n",
    "{\n",
    "\tt = omp_get_thread_num();\n",
    "    \n",
    "    ...\n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which you can programme like a normal mpi template but now communication is done by shared variables and reduction by critical regions which only one thread can access at a time eg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pragma omp parallel default(none) private(t,i,thrsum) shared(totsum)\n",
    "{\n",
    "\tt = omp_get_thread_num();\n",
    "\tint thrsum = 0;\n",
    "\tint totsum = 0;\n",
    "\tfor(i=50*t;i<50*(t+1);i++){\n",
    "\t\ttrdsum += i;\n",
    "\t}\n",
    "\t\n",
    "\t#pragma omp critical\n",
    "\t{\n",
    "\t\ttotsum += thrsum;\n",
    "\t}    \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason they can't all write to a shared variable at the same time is due to race conditions, where one thread reads the variable then preforms the operation but before it can write the variable back another thread reads it.  Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 2\n",
    "#pragma omp parallel default(none) private() shareda)\n",
    "{\n",
    "\ta+=1;   \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. thread1: read a=2\n",
    "2. thread1: add one\n",
    "3. thread2: read a=2\n",
    "4. thread1: write a=3\n",
    "5. thread2: add one\n",
    "6. thread2: write a=3\n",
    "\n",
    "so a=3 not 4 as it should."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Makefiles\n",
    "\n",
    "Compiling the C code and converting it into python modules used a few commands which can be a pain to remember especillt if the C code has multiple files.  Generqallt when writing C this is automated with `Makefiles`.  Here is a short introduction to them.\n",
    "\n",
    "Makefiles consist of a list of rules, usually for updating files when files they depend on change, but they can be used more generally.  The rules take the form (note we need <tab> rather than 4 spaces) of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target ... : dependancy ...\n",
    "\tcommands\n",
    "\t...\n",
    "\t..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `target` is usually the file to be updated or rule name to pass to `make`, `dependancy` are the files/rules that the target depends on and `commands` are the bash commands to run in order to update the target.  A very simple example would be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file Makefile\n",
    "\n",
    "test :\n",
    "\techo \"Hello!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "make test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has no dependencies so when we `make test` it finds the rule for `test` and runs the commands below. `make` always does the first rule if nothing is specified (unless you have set a default) so equvelently we could have typed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "make"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can add a dependancy so in order to do the first rule we need to do the second first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file Makefile\n",
    "\n",
    "test1 : test2\n",
    "\techo \"Hello two!\"\n",
    "    \n",
    "test2 :\n",
    "\techo \"Hello one!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "make"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "both `target` and `dependancy` can be files.  Now make checks to see if any files in the \"rule tree\" have been updated then runs the commands from there up.  A simple make file for the compiling and wrapping the c functions would be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file Makefile\n",
    "\n",
    "wrap : lib/libcexample.a\n",
    "\tpython setup.py build_ext --inplace\n",
    "    \n",
    "lib/libcexample.a : lib/cexample.o\n",
    "\tar rcs lib/libcexample.a lib/cexample.o\n",
    "    \n",
    "lib/cexample.o : lib/cexample.c\n",
    "\tgcc -c lib/cexample.c -o lib/cexample.o"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "make wrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is OK but a bit verbose.  Just like when writing python we are better to use variables to make the code simpler to understand.  In make files variables are created with `=` and `:=` signs. The values are accessed by $() construct so if `var1 = 2` the $(var) evalueate to 2. The first assignment `=` is `implicit` which means that it doesn't expand the rhs immediatly, the second, `:=` is `explicit` in that it does expand it before assignment.  The difference can be seen in the example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var1 = $(var2)\n",
    "var2 = \"hello\"\n",
    "echo $(var1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we don't expand `var1` until we get to `echo $(var1)` so it doesn't matter that `var2` isn't defined when we assign `var1`.  With `:=` this would matter as we would try to expand `$(var2)` when creating var1 and it wouldn't exist.  Conversely in this example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var1 = \"hello\"\n",
    "var1 = $(var1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the second assignment is implicit this creates an infinite loop that can't be expanded.  Here `:=` would work fine as we would expand it before assignment.\n",
    "\n",
    "\n",
    "There are also `?=` which assigns the variable only if it has not previously been assigned and `+=` which will add another element to a list, ie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var1 = one two three\n",
    "var1 += four"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also create pattern specific variables using `%`.  This will match any non-empty string and can be used in any string object but only once. With the automaticly defined variables `$@`, `$<`, and `$^` which are the filename of the target, the first filename of the dependency and all the dependancies this allows us to set up generic rules for all files of a specific type like object files which are always created from their c files, ie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file Makefile\n",
    "# compilers, flags and libraries \n",
    "CC = gcc\n",
    "CFLAGS := -g -O3 -xHost \n",
    "\n",
    "# librarys\n",
    "LIBS := \n",
    "    \n",
    "# Directories code objects and librarys     \n",
    "LIBDIR := lib\n",
    "\n",
    "# source file(s) without suffix \n",
    "CFILES = cexample\n",
    "\n",
    "#This says don't look for a file called \"wrap\"\n",
    ".PHONY : wrap\n",
    "    \n",
    "# dependany says wrap depends on all files in LIBDIR\n",
    "# matching anything in CFILES but with 'lib' on the front and '.a' on the end\n",
    "wrap : $(LIBDIR)/$(CFILES:%=lib%.a)\n",
    "\tpython setup.py build_ext --inplace\n",
    "\n",
    "# for anything in CFILES but with 'lib' on the front and '.a' on the end\n",
    "# dependent on anything in CFILES with '.o' on the end\n",
    "$(LIBDIR)/$(CFILES:%=lib%.a) : $(LIBDIR)/$(CFILES:%=%.o)\n",
    "\tar rcs $@ $^\n",
    "\n",
    "# everything that ends in '.o' should be made from the same file with '.c' instead\n",
    "$(LIBDIR)/%.o : $(LIBDIR)/%.c\n",
    "\t$(CC) -c $< -o $@\n",
    "    \n",
    "clean : \n",
    "\trm *.c *.so $(LIBDIR)/*.a $(LIBDIR)/*.o"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "make clean\n",
    "make wrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fortran\n",
    "\n",
    "As an aside we note that something similar exists for Fortran code, the f2py package (which can also handle C).  It's only for wrapping code where as Cython alows you to mix it together and in fortran strings and arrays need a little more work to pass (see: https://docs.scipy.org/doc/numpy/f2py/).  F2py can be used on the command line a bit like a compiler, first you create some Fortran code like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MODULE funcs\n",
    "  \n",
    "  CONTAINS\n",
    "    \n",
    "  SUBROUTINE fibonacci(n,m)\n",
    "    IMPLICIT none\n",
    "    INTEGER, INTENT(in) :: n\n",
    "    INTEGER, INTENT(out) :: m\n",
    "    INTEGER :: i,a,b,tmp\n",
    "    a = 0\n",
    "    b = 1\n",
    "    DO i=1,n\n",
    "      tmp = a+b\n",
    "      a = b\n",
    "      b = tmp\n",
    "    END DO\n",
    "    m = b\n",
    "  END SUBROUTINE fibonacci\n",
    "  \n",
    "END MODULE forfunc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then to create a loadable module (`.so` file) you would then compile with f2py:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2py3 -c fibonacci.f90 -m fortran"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you should be able to use the function after importing the module in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd Fortran/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fortran as fort\n",
    "\n",
    "fort.funcs.fibonacci(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otherwise you can use it via python by creating signiture files.  I'll leave you to explore further with the documentation (https://docs.scipy.org/doc/numpy/f2py/getting-started.html#the-quick-and-smart-way)"
   ]
  },
  {
   "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
}