03_higherperf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higher-performance computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are various ways to speed up the execution of code in a Jupyter notebook.  \n",
    "We will start scratching the surface by looking at\n",
    "\n",
    "- Numba\n",
    "- Fortran magic\n",
    "- Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other ways to increase performance include"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Numexpr** - a package that speeds up complex array operations (sometimes slow in numpy due to creation of many temporary arrays)\n",
    "- wrapping C libraries with **ctypes**\n",
    "- **CFFI** (C Foreign Function Interface for Python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy offers fast vector computations of arrays, but some algorithms can't be vectorized and require explicit loops.  \n",
    "However, Python loops are slow. `Numba` can speed up Python through a *just-in-time* (JIT) compiler, compiling Python code directly to machine code. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random walk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First without Numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def step():\n",
    "    return 1. if random.random() > .5 else -1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def walk(n):\n",
    "    x = np.zeros(n)\n",
    "    dx = 1. / n\n",
    "    for i in range(n - 1):\n",
    "        x_new = x[i] + dx * step()\n",
    "        if x_new > 5e-3:\n",
    "            x[i + 1] = 0.\n",
    "        else:\n",
    "            x[i + 1] = x_new\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 100000\n",
    "x = walk(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x118007b10>]"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAECCAYAAAAsBKpOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcU9X5/9+ZlWUWYBj2fTsgCArIjix1wQWli1X71bYo\n4lZttdaq/dnaxVZrpYqtrQu4t7Zqba1URAGVTUAUWYTDALIvDsPAzDDMPr8/kkxuknuT3DvJTCZ5\n3q8XL3LvOefm5Exyn3vOeZ7P46qvr0cQBEEQYkFKc3dAEARBSFzEyAiCIAgxQ4yMIAiCEDPEyAiC\nIAgxQ4yMIAiCEDPEyAiCIAgxI81JI6WUC3gSGAFUAHO01rsN5TOB+4Fq4Dmt9bNWbZRSZwBPeZoW\neM7XKaVuAOZ6rvGg1nqRo08oCIIgNBtOZzKzgEyt9QTgXmCet0ApleY5Pg+YCsxVSuWHaPMgcI/W\nejLgAmYqpToDtwHjgRnA75RS6Q77KgiCIDQTTo3MJGAxgNZ6LTDaUDYEKNBal2itq4EVwJQQbb6h\ntV6llMoAugAngTHASq11jda6BPcMZ7jDvgqCIAjNhFMjk4PbGHipUUqlWJSVAblAtlkbrXW9UqoX\nsAXIAz4PcQ1BEAShBeHUyJTgNhoN19Fa1xnKcgxl2UBxqDZa631a60G492b+iNvABF7jhMO+CoIg\nCM2Eo41/YBVwKfC6UmocsNlQtg0YoJRqB5QDk4FHPGVBbZRS/wF+rLXeCZQCtcB64EHPElprYDDu\nmU5I6uvr610ul8OPJAiCkLTE7MbpciKQafAU8+6TzAZGAW09nmSXAL/A3fEFWuu/mrXRWu/wGJw/\nAJW4jdIcrfVRpdT1wI2eazyotf53BF2rLywstf15EpH8/GxkLNzIWPiQsfAhY+EjPz87voxMHCNG\nxoP8gHzIWPiQsfAhY+EjlkZGgjEFQRCEmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCEmCFG\nRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCE\nmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZYmQEQRCEmCFGRhAEQYgZaU4a\nKaVcwJPACKACmKO13m0onwncD1QDz2mtn7Vqo5Q6C5gP1ACVwHe11oVKqceAiUCp57KXa629rwVB\nEIQWgNOZzCwgU2s9AbgXmOctUEqleY7PA6YCc5VS+SHaPAbcqrWeDrwJ/NRzfhRwodZ6uuefGBhB\nEIQWhlMjMwlYDKC1XguMNpQNAQq01iVa62pgBTDFpM0oT/0rtdabPa/TgArPrGcg8LRSaqVSarbD\nfgqCICQ0ldW1fPT5ISqqapq7K6Y4NTI5wEnDcY1SKsWirAzIBbIDztcqpVK01kcBlFITgFuBPwJt\ncS+hXQPMAG5RSg1z2FdBEISE5d8rdvP8O9t5denO5u6KKU6NTAluo9FwHa11naEsx1CWDRSHaqOU\nuhL3fs3FWusioByYr7Wu0FqXActw7+UIgiAkNTW1dWzceYzaujoqq2p5d91+APZ/FZ87Co42/oFV\nwKXA60qpccBmQ9k2YIBSqh1uYzEZeMRTFtRGKXUNMBeYqrU+4ak3CPiHxykgDfdS2/ORdCw/Pzt8\npSRBxsKHjIUPGQsfLXEsZv74PwB069iWWVP6N5xPTU2Jy8/j1Mi8CZyvlFrlOZ6tlLoaaOvxJLsT\nWAK4gAVa68NKqcA23/cssT0O7AXeVErVAx9qrX+plHoRWAtUAS9orbdF0rHCwvi05k1Nfn62jIUH\nGQsfMhY+WvpYHDp2ihMnTzccV1XXOv48sTROjoyM1roeuDng9A5D+SJgUQRtAPIs3uNR4FEn/RME\nQUgG6urqG17vO1rWjD2xRoIxBUEQWiiVNXXhKzUzYmQEQRBaKG9+tDt8pWZGjIwgCEIC0Ll96+bu\ngiliZARBEJKM+vr68JWihFPvMkEQBKGJqa2z3oMpO10d0TVeWqJZ/ulBABbeMz0q/QqFzGQEQRBa\nCAsXbbcsO1VREzRDqa6p5Xcvb2DdtqPsPHCS6pq6BgPjbhOZYWoMMpMRBEFoIXyx93jI8o07j3H2\nwPyG4w07Cik4cJKCA25FrwnDuvjVv+2xFTGfzchMRhAEoYVwsqwqZPnLSxrCFVm64QBPv/WFX/nq\nLUdi0q9QiJERBEFIEIpLKxtev/LejhA1mw4xMoIgCC2EVhmpzd0F24iREQRBaCFUVUcW4X+wMH4k\nZsTICIIgtBDqIohvKTtdzf0L1jVBbyJDjIwgCEKcEiouxorjJRUx6IlzxIVZEAQhDik5VcWPnlgJ\n2AuafOC59bHqkiNkJiMIghCH/GOZL53ypl1FfmWDeuQ2dXccI0ZGEAQhDtm+r7jh9WOvfc7itfsa\nji+f1Jff3zS+4bh1pnOvs1jrmImREQRBiENq6/xv/v9c7pvZuFwuOrbzqS736ZLj+H0KT8Z2D0f2\nZARBEOKQklPW0f29u7jTJd97zUjSUlOoqKpl295iy/qhOB5jIyMzGUEQhBZEZnoqrTPd84OBPdrR\nt2sOQ3q3Z3CvdhFfY/bFgxte7/8qtjE1YmQEQRBaECkWd+3Dx8tNz9942dCgc8P7d2x43TWvTVT6\nZYWj5TKllAt4EhgBVABztNa7DeUzgfuBauA5rfWzVm2UUmcB84EaoBL4rta6UCl1AzDXc40HtdaL\nnH5IQRCElkaXDm04YmI4TlfWmtZPcblMz3++61jQuTrDfs/zi7czbWwfZ52MAKczmVlAptZ6AnAv\nMM9boJRK8xyfB0wF5iql8kO0eQy4VWs9HXgT+KlSqjNwGzAemAH8TimV7rCvgiAILY7pI7vbqt/N\nZEbywOxz+Hjr0aDzuW0zGq5/vKQyqDyaODUyk4DFAFrrtcBoQ9kQoEBrXaK1rgZWAFNM2ozy1L9S\na73Z8zoN9yxnDLBSa12jtS4BCoDhDvsqCElLSXkVtz++Ar3P2aaw0HKYMbZ30LmUFP/ZzfwfTmbh\nPdNJSXEFlcUKp0YmBzhpOK5RSqVYlJUBuUB2wPlapVSK1voogFJqAnAr8McQ1xAEwQY/mr+SstPV\nPPy3z5q7K4JN6myGr6SlBhuNDtmZjB/qS1RmVHFObSIj49SFuQS30fCSorWuM5QZnbazgeJQbZRS\nV+JeQrtYa12klDK7xolIOpafnx2+UpIgY+FDxsI3BjIWPuJ5LFq3zrAsM+v3UcOy15XnDeLrUwfQ\ntnU6nfLaNpzv2sX3rJ7VNjNKPQ2NUyOzCrgUeF0pNQ7YbCjbBgxQSrUDyoHJwCOesqA2SqlrcG/w\nT9Vaew3JOuA3SqkMoDUwGNgSSccKC0sdfqTEIj8/W8bCg4yFm8LCUhkLA/E+FsUnzb3FwPw+d7zY\nV/+CUd0pL6ugvKyCcUPy+c9Hu4LaVVRUR7G31jg1Mm8C5yulVnmOZyulrgbaejzJ7gSWAC5ggdb6\nsFIqsM33PUtsjwN7gTeVUvXAh1rrXyql5gMrPde4T2sdOu+oIAhCAvHvFV8CkJmRSmWVuUeZEWMa\nAJfB0ywvp5Vp/RiryTTgyMhoreuBmwNO7zCULwIWRdAGIM/iPRYAC5z0TxAEIVGorKrliqn9ee2D\nXSHr9e7sXkI7b3QPv/NWey9rvzgSnQ6GQWRlBCFJ2HnwJCcraslt1fJS+CYzV04fQM9OWWHrZbVO\n59m7pwV5jbks4mdKTzfNcplE/AtCkvDblzZwx2MfNnc3hAj52kj3jOSMPh34bEdwQKUZdtySI03l\n3FjEyAhCAqN6ButZHS461Qw9Eeyydps7iDI1xcWWL4vC1I5fxMgIQgKj9wd7/v/smbXN0BPBSF19\nPTW1oWcSZZ7lrNLyKmZfPKQpuhUTxMgIgiA0MQ8sXMfcRz6IqG5Kiot22b6Ylu/NUDHqVWyQjX9B\nEIQmoraujh//eXVDrpi6+npLYUsvqSkppBn2WiYP72b7fa+7eAiZGf4OH20y0yivrLF9LbvITEYQ\nBCHGVNfUsftQCSWnqv2SkVVX17F0wwGefHNzyDTIaWm+W7UTzbFJw7tyzuBOfud+c8NYAK69YJDt\n69lBZjKCkKDUNVW0nRCW59/ZzpqtR5g5oY/f+ROnKnnlPXeIYWV1La0ygm/JPfLbcqoi+jOOdlmZ\nLLxnetSvG4jMZAQhQdl9sMSy7O/vF1BSLiIaTcWare7Ax/+u3uN33mhUXJjPUDLSU/2ELVsaMpMR\nhARlw46vLMve+2Q/nxUU8vubJzRhj4RAjEtk9fheV1S5Zy7efZTWmWnccOkZdIlxFstYIEZGEBKU\ncMmojp2saKKeJA8l5VX8aP5KgIbo+1WbD1vWr6z2aZLV1LqNTHVNHbfM+8hdbtAsGz+sCy0RWS4T\nhARl/XbrmYwQGzYW+CLzvZ5bCxZts6z/yN99eX6OFLlVlBMtWFaMjCAIQpSoNQmwzMuxzttinG12\n8NRLNH8NMTKCkAQ0hReRANW1Pgvx/if7Aejf3V5SX+PeTCIgRkYQBCFKnCzzzUzeWrWH+vp61m2L\nbNmy2NNWZjKCILRIbv/m8ObuQsKz76h/xsodJtpxVjz44gYAausSy8qIkRGEBGfuzDMAOKNP+2bu\nSeKT3Taj4fX5o3vy1YnTIetfOKZn0LnM9JYbE2OGGBlBSHD6dcsBID1Nfu6x5uOtRxteD+yRy7ov\njoaoDZ9sLww6V2eYyVwxtX/0OtdMyLdOEBKcqhq3x1NghsSOuea534XoUFtXT9vW6SHrFJUExyoZ\nl9guGtc76v1qasTICEKCU3rKXD5GgjFjS1V1bcSb/kY27Aie3bRkHEX8K6VcwJPACKACmKO13m0o\nnwncD1QDz2mtn42gzTxgu9b6ac/xY8BEwLuTdrnW2n9XTRCEsPTrZs+FVogO/175pWXZd2coXlys\nTcsSTdjUqazMLCBTaz1BKTUWmOc5h1IqzXM8CjgNrFJK/QeYZNZGKdUReBEYCGw3vMco4EKt9XGH\nfRQEAYLyiAhNQ3GptazPmMGdTI3MyVNVDTPPMUM6BZW3RJwul00CFgNordcCow1lQ4ACrXWJ1roa\nWAFMMWkzylM/C/gF8JL3Ap5Zz0DgaaXUSqXUbIf9FATBQG5WRvhKQkzp2zXH0vDf8cRKjha7PdIO\nFCaGvIxTI5MDnDQc1yilUizKyoBcIDvgfK1SKkVrvUdrvR78dK7bAvOBa4AZwC1KqWEO+yoIScei\nNXtMz58sE3n/5uTcEV25/3ujSU0Jf+s9dCwxjIzT5bIS3EbDS4rWus5QlmMoywaKw7QJpByYr7Wu\nAFBKLcO9l7MlXMfy87PDVUkaZCx8JNtYvPFhw3ZnyM+e264NGQkWl2GHpvpeDB/QkV/dOIFUG1kt\nr79saEJ8b50amVXApcDrSqlxwGZD2TZggFKqHW5jMRl4xFNm1SaQQcA/lFJnefo4CXg+ko4VFopv\nALh/PDIWbpJ9LEJ99k+2HGJgj3aNuv6rSwtYsn4/v7p+DD3ysxp1raakKb8XlZU1HC8qs9Xm+InT\nTda/WBozp8tlbwKVSqlVwKPAHUqpq5VSc7TWNcCdwBLcxmiB1vqwWZuAaza4VGitt+N2BlgLLAde\n0Fpb62ULgtDASQuXZXA/URsx5jNxypL1biHI3728odHXSlS8AbF2SLMx64lnHM1ktNb1wM0Bp3cY\nyhcBiyJoYyz/VcDxo7iNkSAINrjjiZWWZT+bPYYrf/a/huO8nOgFZJ6ubLzBSgT6d8th1yH/1Nez\nJvezfR07S2vxjARjCkIS0aZVup/s/+6Am6HgnLc8cTFHjpcHlTmR9MlIENdzMTKCkEBU10Q2m/BK\nyoTK2ijYwxt8eaqiJqL6D8w+J2R5iktmMoIgxBk//vPqiOoN758Xk/ef98+NrNxkndM+Gbn2QmV6\nvlfnlu85FgliZAQhgSg7Xe13POfSIab1tI08J6F44Ll1fsdbdh9n4f+Sb3a094i1F9i0s7s7umai\nyMuIkRGEBKHe5KY0fmgX07oHTaLJ1207ynUPLeNURTU1tXW8urSAw0WhAwL3HbXnlpsIrP3iKB99\nfqjh+K1VX/LL59f71clIj+zW+uc7zuWWWcP4y51TGNzL35W8PkGSlzmNkxEEIc4we/INlPcPxV//\nsxWAF97ZTnpaKmu2HmHJ+v3ccOkZjB9mbqySkafeco/TuSO6AfDvFcFCmHk5rThcFOwAEEjrzDRG\nD3ZrlE0b2YPt+3wzzCF9OkSju82OzGQEIUGoqPLf9P/2tAGWdUMt4ew9WsqarUcajp95+wtKTGJv\nEmU5Jxb88IoRttuMHOSLYbr6awPp0qFNNLvUbIiREYQE4bXluyKuW2hICxy4zFZ4IjjPTE1tsALU\nmi1Hgs4lE6E8+Tq1a237eqkpKQzp7U6R3btL4jgFyHKZICQIxn0CgPwQN7qjxb6lnBNlVbTPzvQr\nz83K8BPTrAvYH6iprUt69+cfPbGSP98xJej8Pf83EoDf3zSeWpuzvdu/OZwDhWX07544OYBkJiMk\nDMYb4eGiU7zx4S5q66w0WBMf4/JLIO2yMi3LAHp28tcgS031v1W8tco6IVeyYKVw0MczC+nYrjWd\n29tb8srMSE0oAwNiZIQE4Z6n1jDn98s58JXb2+lnz6xl0Zq9vPCOefbBZCDUpv/si32uzWaGuLY2\n9BN44P5PIvDpjsKQrshg7sEXSGpqYgRRRgsxMkJC8JUn0dMn2j+n+srNEhhohnFT+a//2Rp089y2\nt9jvOHC5LFE2pb3U19fzp39tDnJFDkTvCx9flCiR+tFCjIyQUOw+nLxaXGkOn6B3HyqhNkxMRuBs\nJ9FupOE+v5cdB8IbGTtu48mAGBkhodiy+3hzd6FZePOj3dQYlrjsKPj27pwddnlsY8Exv+P/rPTt\nycwY2yvi94pXIjUy67Z9Fb6S4IcYGaFFs3rLYTbtKmo4zm6T7ld+3ugeTd2lJqemto7/rt7TcNwh\nJ5OnfzI1bLtLJ/QG4KJxvagJ4yARKENjzFkzsEcu83842a88cHkt3om0v2YpkTvkhHaiSHbEyAgt\nmmff3sZjr33ecHzJ+D5+5e9/cqCJe9T0BO5F3/7N4REt2XjVgl9bvstUZsZIqCf9tNQUslr7G/fb\nHl8R9v3jiUhnMr0DRC036EKOl1TGoksJgxgZocViFnHeLiujGXrSfFRV17LlyyK/cykRLpUt//Qg\nAEUlFX4R/mYERvx7UwWAuSfa6crI5O7jhUiNzEXj/JcG//ymfxb5K6b2j1qfEgUxMkKL5eFXPg06\nV18Pyz9N/NmLl5eX7OCJN/xvdE4yKn648VDI8j0Brr3HTvpUAea/scn2+8UbxuWyUG7KoZbVslqn\nc9G43lHtVyIgEf9Ci6XgwMmgc0eLy00FCxMVMxftSGcyTnlt+c6YXr85MHrP1dbVW3rqWY1tr05Z\nPHDdmJj0raUjMxkhoUgmA3PUJM0vQGHxadPzgYwZ0snR+76zdp/f8fWXmOesaUnsOeybqYXytEtN\nMb9lHjRxCBDcOJrJKKVcwJPACKACmKO13m0onwncD1QDz2mtn42gzTxgu9b6ac/xDcBczzUe1Fov\nctJXIbnITE+lsjrxotHNMJvJAUTq1zW0T4cgl9yLxvYKMiLhyEyPfS76L/YcJ6dNBj0C5G6ihXFP\nav9XZQzoYS7tYqU8HemeTjLidCYzC8jUWk8A7gXmeQuUUmme4/OAqcBcpVS+VRulVEel1P+AmYZr\ndAZuA8YDM4DfKaX83VeEpMbqx95SDUxVdS3XPbSM9z/ZH3EbKwcy46Z8KMyWfgLVfwOTnpWWB0v+\nnzXQWiMtWvzh1Y38fOG68BUd8pkhDui3L2+wrFfVQr9fzYlTIzMJWAygtV4LjDaUDQEKtNYlWutq\nYAUwxaTNKE/9LOAXwEuGa4wBVmqta7TWJUABMNxhX4UExKgQnAgs83h6/e39gkZfK1BR2YpO7YNV\nmgOdBgKX1ErK/dM7g9uFGeCH3xpOZoZ7VtMjv21EfYiEWOetue6hZRHXTXblaSc4NTI5gHGuXqOU\nSrEoKwNygeyA87VKqRSt9R6t9XrA+O22uoYgANZP8WZEImrY3DhxvS41ueGD9b5BIAN7tAs6F7js\nk5Hmu5beV0zbVtYr7CMGdGTerRMB6JhrP5+KFadOm39Ouzz/znaue2iZX24cO9+Ng4XJl2o6Gjj1\nLivBbTS8pGit6wxlOYaybKA4TBuz6wdeI7xoEJCfnzjJfhpLIo9FWmbkq6c/+MNy/vyT6THsTePp\nmOe7gUX6d/unhZdXt66hn8dCXf+1D/wTn+Xl+fZAHv7bZ1ww1t9Fd1j/PL/rVVS542NS0lKi9v17\n+e8+V/Wy6jr6dnP2vOnNt/PjP6/mb7++CIAOeeZ7PGZ93xpGHDORf2+NwamRWQVcCryulBoHGB31\ntwEDlFLtgHJgMvCIp8yqTSDrgN8opTKA1sBgYEskHSssDC3VnSzk52cn9FgUl0YeZb3vSGncj8XJ\nkz5Pscb0tXVmasj24b4XFQFBlK0CJkVL1u71Ox7YPdfvet5ZQkVFddTG/Lghi+ftj37Awnsa98DQ\nLiuDwsJS8vOzOXLUXFD1lUVbuWCMf+BlaWlwxlAj8f4dC0UsDaTT5bI3gUql1CrgUeAOpdTVSqk5\nWusa4E5gCW5jtEBrfdisTcA1G+atWuujwHxgJfA+cJ/WOrEW4ZOMzwoK2XMkegrJm3cXha8UR7y3\nfj9/+tdm0+WZr06c9tMeA/cyzpL1+zlVYb1U1K1j8L5HY3O0/eAbZzJuaOeG45y2oZfxxp7R2e/Y\n60wQTe2yjTuPha9kg/1flTX8HWpN0koDvLoseJYoHmTOcDST0VrXAzcHnN5hKF8ELIqgjbH8VwHH\nC4AFTvonxBfVNXUNUelWT6FvfrSb9tmZTD27e0TXfP6d7VHrX1Pw96XuDf0d+0+gerX3K7vnr2v8\njuvr6/nbewUs/fQAry4tsByzLJP9kUE9g/dZ7NAqI43rLh7Cx1uPhq17yfjeQbnsU1wuXMT/Dbns\ndDWdgM02VLuNm/63zBrGk//2La4EeuEJPiQYU4g5VTXh3T7/u3oPL74bPotlTW1dozbyj5dU8OrS\nAspDzBBiycN/+yxsnSPHy1kagTTODpM4mcbO8FJTXaSlpnDzrGF857yBIet+c4qFTpcLv831xtK3\na074SjYpr6zhxf99wfLPDjpqP3qwz+vu2bunccPMM6LVtYRDZGWEmBMupa0d5j7yQaPaP/XW1oYg\nxqu+Fvom2lw0Z9Krgv0n6JGfxTmDnakBgFs/7svDpWzdc5yhfTo0qj/HSyr40kEiul0HT1JcWuln\nDIys3HSYRWv2mpbZJdYyPi0dmckIMWfXodA3CaNuVHVN8BPwvz7azR1PrOSAhQvpZRP7cO81IyPq\nS1GJe/M21F5HtKirq+fuv6zm1aX+sS97j5Q2GN7KquBZ3qOvhp/tGJk1qa/zTgbQ2SStcp8uzjaF\nX1myI3ylAEpOVXHTox/wwUb3DOOPhjQOdnjwpQ1+y1mBSPbKpkOMjBBzRhoiws2WUbyBiGAet/D2\n6j2cPFXFzxeYR3zPmtzPNOYjkKKTFQ25P5woFdulvLKGYycrWLLeP4r/l8+vb8glX20yHkUR5Cd5\nYbH/ntSsyW5D841z+zntLgBDercPOvfDbzmLgz5ioa0Wio07j1FVXceLizXVNbVh89wY+ceyAt74\ncJdfmoFHX/2MT7YHZ7N8O8DRAoJjr1pCfFVLQJbLhCaltraetACpK2NmS7u/629Oifym+pO/rG54\n3RQ56iOJVHfqhWWU5p9ydneyWqdxRu8O9O3WOFdUsyf81NSmexY1pnm+8Q8fmtY5eaqKXBOvt3fX\nuY250eNt655itu4p5rKJfcK+d9tW6ZQZAj+Pl1TSISeTeswVJub9YGKLy5vTHMhMRog5Rk+jWhMf\n236GjV27EiJf7Ck2PT9hWBceuXlCw3HgU+nhIvtP2XYJ91FOVVRH5IV1559Whnyqzm2bQWpKCgN6\n5EYc7W8HM4PcNS94WS0aROKufMcTKxsCK80wm/G+tWpP2Otee6HyO/7JX1bzxBubef6d7fz4z6uC\n6rfLyqRrXvTkcxIVMTJCzDEajmoTGfV1huUMu66vZ/QJXt4BmHPpGeQZhCKNSbagaaTZw81S6urq\nI/LCOlFWZam43CqjcQrIxmVDK2FNs6XFu646u1Hv21gCXdhPlDUuBfKZ/fI4Z3CnoLifjTuPsXKT\nf86eO749olHvlWyIkRFijtFwnDCJ1DfmRbG7fDRjbK/wlYDte/1nPGOHdLaoGT3MZm1GKqpqeeqt\nrRFdy5icrOCAT96kwsRxwA7P3D2NS8a7pWIundDHtI6Z91SkIpx22B3GQSQUd/4peKZhhwHd3bPp\nuRG4IgcqVQuhkT0ZIeYYDccn+quQP1K7Mxmz5SEzb6ute/yD7lItMh9Gk3Cf5acBQZihqPF43dXX\n1/O7l4PTTjeGb5zbj8kjupFvYyYTC0pO2RP1OFhYRvf86OSX6efRQ4vE6yynjX0x02RGZjJCzDEa\nmXBxE4EzmfUmnkFWZLV2i2b26RpsxAKTczWFB2uoDIt2+fgLdwR+LCLpXS4Xndq1trzB2o0DuWr6\nAEf9sOvNdcKmUQrFoJ4i8h4rxMgkGKXlVbz50W7KTJJLeSk5VcW/PtrVJLEi4H9jXG3IQGhG4Mb/\nX0LEOgzt62+wfj1nLHNnnsGZ/fIazqlGyqw0BrtBhMP6hQ9cjAe5lnD7QOef0xOAQRbZJa0ot+mp\nleYxfnc92bilMoD0QJdHIWqIkUkwXnpX89/Ve3h+0ReWdX70xEreXr2X2x5bwa5D5hvK0eR0pW/f\nIHATNZDGCCvmts1g3NAufk/k/bub3+iMsTmxYttec883KyLRv3rlPfsBjtHAGD8TbtPfO/5msjeh\nCEwIFs6DraKqlu17ixtin4T4RIxMguH1ojpmkEcPhTHeIla8/oF53hMIzkpo50ndLBVwIFbLYmbK\nAtHGzmfKXogtAAAgAElEQVT5/kWDw+59rNlyJKyRjhXnje7R8Lpft+hriZkRzni8+K7mVIXEqcQ7\nYmQSDO9NNdLbm1k0dLQor6hm96ESjhZHZvDALcMeKfuOOs9U2D8GN8q6unq2fnmc8ooajhwvZ+0X\n4ZWMvbTLygwbIPrM29az01hz1oCO/PQ7Z/P47ZNi9h7TR/orcFdWh/acKy6tpLNJCmm73HjZ0EZf\nQ7BGvMsSDveNqtwkZW1FVQ0Z6akB5xrnAhuK2x5fYTuC/6m3tgbFKjQGKxHEPl2ib2Te+GgX73y8\nz7L8jD7tLYNH9xwpoYdNT6nWmWk8dOM4W22c4nK5glIUhG9j7z2ceIqdbOTm/73XjAySJJp98WCe\n+595KolYuG4nOjKTSTC8exrb9xb7PUmfrqzhlnkfMefh5U3Sj5NllZYGpixKOdsbw76vop/FcFWY\npSwrAwNw9PjpoJnMj688K+T1WmWkkh2n7rR5OZl0yDZ3iY4mLzdyj8ps/CYP72ZZf2YE8jSCP2Jk\nEoy9R303z5WbfPstx0tCp46NNqGCDJduCJ8rJdZYRdA3hpJy58azrr7ez1V4eP88hvbtwM++O8qy\njZ0U1E1NSoorIomg8opqfv+3TzlRVmnp9HHBOT356XfO5uJxvYPKGhvpb9c9u5vIyNhGjEwCU27w\n6rJKzpSXE5vp//Z9JyzL/rPyy4iuYRYB/uANYx33yYy6unq27S2OapItJ4xWnSiv9Bkpr2ho/24t\nM36j8EQFxaWV6H2hPex+8NgKtu87wZ1/WsWKTeZOKN+eNgDVq72pGKpZqgTweaYZU0mD23gbSbVY\n07vaImHbQJtu2YIYmYTGGKdhJdlhdx8gmlRbZMx85O+fsfaLo/zmxU/8zt/69WG2BQm754eu//6G\nAzzy9894bfkuwD3jqwqz4RwLRql81mwJ7yjQKQob3U1JJJlAvVg5cnhnGy6Xi8dvn0RmeuiYlvx2\nrXjwhnEsvGc6c2f6b+qfNbAjAwxu7VYzmfNH9+RbU4Mzf0oeGvuIkUkS9lhkp2wqd1Qzaiwi4rft\nLQ5abvv9TeMZpexnaxwTRqPMq2n23if7qaiq4a4nV1vmrWksU86yXusH980xHF/Z8NRLRLLbZDTk\nzjEjLTWF780YbFrWv3sO5w7vxs6DvqXSUG7jZstzgn0ceZcppVzAk8AIoAKYo7XebSifCdwPVAPP\naa2ftWqjlOoPPA/UAVu01rd6rvEYMBHw3h0v11pHf7c2gdjpYJ9h064iZk6MXmZFiEwe5OCxU7TP\ninzTumM7Z0/wk87sypsf7bYsN0rLex0SvoowxsguZob+rAEdue2bZwIwfliXJgkSbWq2fFnEsL55\n4StGSH6I78Ljt0+idab/be2q6QN4ddlOvn/RkKCZi6ROjj1OZzKzgEyt9QTgXmCet0ApleY5Pg+Y\nCsxVSuWHaDMPuE9rPQVIUUpd7jk/CrhQaz3d808MTBh++/KGoHP19fUhY2HCpUZ2QiRuy/c/u5by\nJgiky0g3/4oHStIANFaxxaudZkWxwfni2Z9O49mfTuP2bw1vWIKx2viOJOFWPLN6S2gpIbucZci0\nGkiaSYK1C8b04pm7p9K9Y/DSaVMkr0t2nBqZScBiAK31WmC0oWwIUKC1LtFaVwMrgCkmbbxuM6O0\n1is8r98BzvPMegYCTyulViqlZjvsZ9Jz4x8+DJnrHEInCntpieafy60j9u1ez8iLS7St63p56MZx\nDO3bgT/eFj4wMM0iiZfZDf10I42eWbZGIy7DU3OKyxV0gyu0mEFZLduMGpRvs4fNgw7hBGJGuH2n\nUIYhPc38722VzM2OwvRVXzN3BhBC49TI5ADGtZkapVSKRVkZkAtkB5yvVUql4o0edFPqqdsGmA9c\nA8wAblFKDXPY16QmEq+p4hDyHcs/PcjitdYBhmZYyan06uTvZLBlt09+/57/Gxnx9Tu1b8OPrzwr\n7E0dIC3N/CZyyrM0NvVsX5S5cenMSZCfmZPBzbN8X9sLRvcM2d7KNlst6VgZpXjDrqv1b+c2TYAp\nRL5cdvfVZ3PBOaH/foI5TiP+S3AbDS8pWus6Q5lxNzkbKLZoU6uUqguoewIoB+ZrrSsAlFLLcO/l\nhH4kB/LzkzOhUG0jXHC7dskhNyvYlbnCoIprZ1w/1ebLcxNGdGffe+azl4kje8Ir1nlSovV3nTSi\nGys/P8S+r8rIz88mzaC+28qw3HXHEyt56YEZ5GZlROxRFJhOAKBjB5/huXLGECpr6zl/bG/Tz9Oj\nqy9524ThXRvqWC2jTT+nl6NxaY7fiJ337NzJ3xnFTlu7n61Tp2zTJTYvLz0wg50HTjC6CZLcJSpO\njcwq4FLgdaXUOGCzoWwbMEAp1Q63sZgMPOIpM2vzqVLqXK31R8BFwDJAAf9QSp3l6eMk3M4BYSks\nbNlbN+u2HSUvp5WlerAVjdEgKywspeq0/5P7/q/K+MXCdX51ImXN5+ab1yP6tufVEH0I18fG0qtz\nFv933kBWevLDFxaW8v563ywtP0Ay5NoHFgPuZZKvjepuueRixdzLzqCVwdv2eFEZszz7K2afp3t7\nn3fZpKFdwn7mutpa2+OSn5/d5L+RjPQUW+8ZWNes7YDuuX5eYqHqBnL/90bz6xfc7vHHi8rCPkT0\n7timxd9XwhHLBw+ny2VvApVKqVXAo8AdSqmrlVJztNY1wJ3AEtzGaIHW+rBZG8+17gJ+5TmfDryu\ntd4OvAisBZYDL2it/XXAE5S//mcrD74UvIEfDrsaZNltfE/tZstbnxuWjuxitQbvcujJ0ydK6W4f\nmD3GL8ZiQYDg5J/+tTmwCQCvLi3gl899YloWig7ZrWiV4XuOCzchMu4PDO4dXiesX9fmcz+3QziH\niFBYyf1HmnbbjF6dfcu2EvcSexzNZLTW9cDNAad3GMoXAYsiaIPWugC3F1rg+UdxG6OkoTG5VOz+\nVn7x/XO468nVANQY3ndjwTGWrN8XJBoI7pz181/fzLihnUPmPrFKo+vUW/RHV4xw1tDDdy8eQrVJ\nQqxVNryeDhTaV3zOaZvh584dbiZk54Y37wcTaWeyxBmPlJVXs2brET7eepTbv3UmqSkpHCgs44XF\n5iKUAK0zUzldWcv93xttWm43VbMRMSxNi6gwxxG1dc73VQITPllx7QWD2PLlcT812a27i+g00p0v\nZP4bmwCC3ItPV9ZwtLiczbuL2Ly7iDP6dDDdeK+vr6fIQifNqbtoTgQb/KG44muD/JY7unVsy6Fj\npxp1zUjokJ1JqU09swdvGEubzPA/y3g3ML+ZM5b/9+xaAKpq6njmv+5Z44GvTtG7SzYLFm1jr0nc\n0DSP3P/DN03gRGml30zQyMQzu/Liu868E1NcLu789gi6xkCJWwhGIv7jCKsI+Ggx9ezuTBvZg9u+\nOdzvaa7IxLtsX0Bel92HS3AZHAHveGJlUO6X2ro6DhT6bt59u+bwf+cPAiAzPTXIe8o7GxoXRWn/\nSGgKAwOQkZ5KmoVLrRVd89qaOmE09Rg1lm4mMSkA9Z5MR1az9mzP0lpW63R6dLKWPEpLdQUcp5jK\nwFgxrF8eg3uHT3ctNB4xMnFEpPElgUSSIRLggwCRTO9Nfmif8Ov/9QEqwQCbdvn2bQ58VcYNv/+A\nVwyxL1PO6sbUs7tx5fQBPHjDWLLa+K/Ne2dul0zoE1H/452ueW2C9h9y22Yw+6LBPDD7nEZde/bF\nQxrVvjn4+rnBgpZLNxzgdGWNreR0ZgQueT39k6kiAxOnyHJZHFHrcCbz9H+dZUzs4tlU9T5Ubg+R\nk766ui5oT+WND3dzyfg+ADz1X7fWmDGve2qKi9SUFC4cY75J63X5TUsQaY/aunpSU138+Y5z/RwM\nJo8IrVkWCampLW+MLhnfO0jSZ9XmIyH3pt5atYdZk4ONk9BykZlMHGEnJ7wRs7XtSNj6pTsY8tMd\nhZSUV/H7v1sr5q7ZGnqT/GBh8BJUpHswibIPW3jiNKdO19A6My3qmlgpLhe/njOWYX07NHpW1FRY\n/f0/+txc0l9ITMTIxBFW0vfhcJppcsd+t6vx8s8Ohoz6B/d+UbmJh5YXMy2w7WFyiXhZs9Va4v7O\nbzfOs6wpqa+PTGHBKd07tuXOK8+iV+fkDDgWWiZiZOKI99Y3X8bI6jA3x407j/Huuv2mZf9YVtAw\nKzKyaXdR0Lm5M88IOnd2CMHDYf2ip97bWH7/N2tFgsa4nwvOaGnOEMmKGJk4ootF4JkTxprEsXgl\n5c2IRDGgvCJ4xnT9w8ssjc+ZJvLu44Z2YZpBLwxo8idzK8+ncGzfd8LSzbwx7ueCMyTepWUgRqaZ\nWbP1CNc9tIyDx05FLYHYhGFduO/7YxqSYA3p3Z4n7zyXswf6q/Zee8GghtdL1psbinCEcoizMppW\n+xU/umIEl0+Kbm4bM266bGj4ShZYfd5Yu58nC+cMjjwxnU2VH6GZkD9TM+MNUrv/2bWOvcuMXDS2\nF9df4k7OVHjCHRS5bW+xaVBbqORPZkwa3tVWfSuvMasN4eH985rEyISKvwB3KmQrrJbFvE4boZb+\nkpEeYdJfB/Kd8yKX05dcMC0DMTJxhHHJxalsxuDe7SNfRrD5G81IC51bPZA2rcw1q+zk8IgV114w\niN6ds3n6J1ODyqaEcDm28gD0Gp94+GzxxAETr0MjgeoGdoJX7X4fheZBjEwcYbyBbXQoUDnMxMvL\nik07gzfmQ2E3WPTM/uab9idO2csvEgumjezBL2afYyrzPrRvB+ZcOoRHb53Iwnum+5VZeY95/3ap\nIWTjhWB+df0Yv2M7RnrK2e6HAcnzEt/ILyKOMC7F5Oe2ClHTnIduHGdrMzSSusY1cjseVD+5yjqp\n2MchXJaNdO4QPUeISJk+sjsul4sJw7r66bt52bon2IsOfDPPtV9E9tkEN4FGxU46hR75WSz46TTJ\nWBnnJLSROVBYxrEWkj0Q3JIbXnJsCCD26pxFq4xUOrW3d1MeMyT8Jmtvg8y+nZlMenrjlzJCLVvF\nimMnzcU9vZScMo9JKj3tXBU4mQlM/2B3M188zOKfhDYyP1+wjrv/uqa5uxExn++yt3zlpba23nTZ\n56wB7k1ou5uvRozX9W4Zzb54cNh20ch1EiiCGCuM45MXZgbZtpW5ElP7bHe7icOsUyAIwWQHaL3J\nZn7ikbBGpt6h2GS8UG9jaep0VY2pNMvXz+1HblYG116oTNtV1diL7fDOZCJZN2+MrMpNlw8lLyeT\nMU0UbGf0aOsZxvPsg43mWT+9S4mtIpDpF9wz5F/PGRs0E5GZSeKRsEZm96GS5u5Co/jMxsb/8ZJK\n07wlPTtl8ccfTDJNQOYtN/LtaQOC6hh/8t7N7cY+bfbtGjr4csyQzjxyy0Ry2jQuj0ykDOrpG59w\njhO7Dpp/r8S7zB4De+TS3RMUm5Huvg3dMmtYc3ZJiBEJa2SWfWr+xBlPVFZba5UFqtfGgkBZ+nAp\nbV/yJIkKdEs1W6oLxZeHfYKe8bA8kt0mg2fvnsazd0+jY25w7JCZA0AgNZ61xGgLY7Z0bv26W2Xi\nN3PG+p0fbvA8/MX3z+GyiX0YOcg6PklouSSUkTlkSJFbYsixsmCRMyn8WKP3nWhU+/1flfHAwnVR\n6o2bQMkXM975eK/fcZtWaXRqby+w08uAHrmO2kWblBSXpYGIZOlVZjLmjFL5LLxnOt06tuVHVwyn\nW8e2jBnSiTP6+GaMXfPaMmtyPzHQCUpCGZkbH1ra8Lpnvm8paNXmI7y2fGdzdCkk67c1zt31FwvX\nBWWwbCzXXDCI6y/xJcgy++Gnp/t/bTpkZ9qSn7/3mpENr1vCbSWS3bEyz3KlGBlrhvfvyG/mjOWm\ny4fFxQxWaBoc7VIqpVzAk8AIoAKYo7XebSifCdwPVAPPaa2ftWqjlOoPPA/UAVu01rd6rnEDMNdz\njQe11osi6VtNbR0pLheL1+3zO//O2n1cYbLn0JyEkzdpDlwuFxPPdMvHvLq0gNGDO/HKezv86lw6\nvg95ua0aJHHuuupsWmWk0TG3VVgXYICBPdoxsEcuBQdOtoxcMhFYmSf+tRlo/OxUEBINp64ws4BM\nrfUEpdRYYJ7nHEqpNM/xKOA0sEop9R9gkkWbecB9WusVSqm/KKUuBz4GbgNGAm2AlUqpJVrrsIlT\nbpn3EdltzOVMnvvftrhKY5uZEb+yGBPP7NpgbALp3SWbYX07UF1Tx+Be7Wjjceu1E+2e4wnU7GhT\nPy3e0fvFyAiCEafLZZOAxQBa67XAaEPZEKBAa13iMQorgCkmbUZ56o/SWq/wvH4HOB8YA6zUWtdo\nrUuAAmB4JB2rqa2juNRctmTFpsOUlsdP0FzvFpp8qkd+Fi6Xi3NHdPMLAD16vDzia1xz/iBmjOll\n6tEWb7SxiI0RBCE8To1MDnDScFyjlEqxKCsDcoHsgPO1SqlU/JflSz3tA+t6r9FoIlnOaSrsaoHF\nAq/u03dnmMfSmGHlbeX1VsvLCS+Jk5uVybenDwjycItHbvn6mSGVmQVBsMbpI1oJbkPgJUVrXWco\nM4Z7ZwPFFm1qlVJ1AXVPWFwjKusQ/Xt3oENOK558YxP9u+cyY3yfaFzWNm+t2MXnO6xjYfp1zyU/\n395MJ7B+JO1vvuIsLhjfh4E920fs3RPuukUlFbb7Hmsa05/8/GzOGtKFb9+3iNOVNUHXMipmz5zc\nL+4+eyDx3r+mRMYi9jg1MquAS4HXlVLjgM2Gsm3AAKVUO6AcmAw84ikza/OpUupcrfVHwEXAMmA9\n8KBSKgNoDQwGtjjsqx9FRacoLi5n8Zo9AJzZux0ZUdDZsssz/w79cepq6ygsLLUsD3SrvWR8b7/6\n+fnZIdsb6dAmnaKiyL3UrK5bdro6bJ3mwM5YhOJ0ZQ0ABw+d8PvO3PiHDxpezxzXK64+eyDRGotE\nQMbCRyyNrdPlsjeBSqXUKuBR4A6l1NVKqTla6xrgTmAJbmO0QGt92KyN51p3Ab/ynE8HXtdaHwXm\nAyuB93E7BkRlM6W2rs7vBn2oKHS+i+YgLTUlrOJxwYGTfscDujd/vIk3s2dmMxjtpsRrbLxUG+R5\n7AamCkKi42gmo7WuB24OOL3DUL4IWBRBG7TWBcBUk/MLgAVO+meka14bDhf5NqSrquv8XFILDpyk\nT5fopD2OFhlpKdSG2a+pClALiIf4DK+UTyglg0TAKnGZIAjBJPxj18+uHeV3/OsXPuGp/25tOC6K\nI0cAL2mprvC5WwJtShPZmP/33dHhKyURpyrCetULQlKT8EYmMAlSZXWtn8jhBl3Y1F0K61VWUl7N\n4aJyW0rSrhhamafumgLA1yf3bVgSM6OjRyY/3UYK3ZbE0D7t/Y7r6uq57bEVFrUFQYAkMDLhAh5r\n6+zJ3UeDSDNM/nf1HsuyTwONYwxnMulpqSy8ZzozJ/YNWW/uzKEA3Pr1xFTTbedx3a7x/P3e+2S/\nX/nF43o3eZ8EId5JeCMTjv7NsGEeqZF5/5MDlmUfbDzkdxwPf8gBPXJZ8NNpDO/fsbm7EhNWbT4C\nwBHPHl9gKuamSrImCC2JeLg3xZw+Xazd85pjuSwWG8d2Uy/HimRIOvXs2+aq3h9/0TjBU0FIRBLa\nyJw90P1EHW62csSGHEo0CNxr+dqoHvw6IN8G+LvGhuL+740OmzZYiB7eZG+B+2BmieMEIdlJaCNz\n4Zhenv97hqx339Mf89aqL1m8dl/IetHi811FfscHC8vomueeiRg3zSN1Be7bNb5csBOV80b1AHzf\nq9y2/pk7wyV9E4RkJKGNjPeG3TG3NUPDpNX994ov+efynbY8uiIl8JpeiXwvOw+WkOJysfCe6Tx1\n11Rb187LCZ+1UYgOuw65A2A/3OjOuhqYyXHycHPVakFIZhLayBgDFM8Z3CmiNjsPngxfyQaV1bVc\n//ByXlqiLeuMHRJZ38yItxw5iYw3bfRnBW7NucDtp3ZZYvAFIZCENjIHj/kkYwb3ahdRm2Mnohuc\n+VXxaQCWf3rQss50zzKMEwLjgITYcdV0t0GfOaEPQJMtrwpCSyah71DG/CaRpnutqrEviXLrHz/k\nuoeWmZZFsvwWak/lZ898zP4QKZZHDMgL30EhKngTrHljryRBmSCEJ6GMTF5uKy4ybL4O7OmbvUQq\nY+8k9/jpSmvDdKLMX9dz+WfWMxozDheV84uF6/jycAmvvLcjKHhUBBmbjuMl7lnu6x/sauaeCELL\nIaHuUM///EK+NbV/w3HfEPExVrTOtNYMDScHs/dIsGy4cWP+aHE5L71rvTfjpUd+VtC5X7/wCUs3\nHGBTgGea0HSESnhn/N4JguAjoYwM+AcDGg1GpFkoreJNrn94GXMeXm65LAbwy+fXs3GnfyIyY977\ne5/6OKjNzbOCJVi+OaWf5XvU1Lo/R9tWafTIb2tZT4g+xlnj5wF/Z5GUEQRzEs7IANx0+VAum9jH\nz+BEKiBpFY0fqWfz/Nc3UVruWyILJyFj5vUWKoma91PU1tVHvAQoRAejbMzjr29qeO2NcRIEIRin\nmTHjmjFDOgedqycyK1Fba08w07tOb8QYqf/Ox3ttXQ9gUE9rhQKv3ayrq3e0fyQ4Z/zQLry1ak/Q\n+d+YqDUIguAmIWcyZrTPjiyG4e3Veygtr+KV93ZQXFoZtr5ZnWqPodq06xirthyxbHvWAHMhyVBu\nyacq3FkZ6+rr4yJRWTJhlcIgGfTaBMEpSWNkjDfu1plpzP/hZH5/0/igelv3FPPq0gKWbjjA4699\nDsDGAv/1d6Nbstky2r1PfczSDQfYvOt4cKGB6y4ZYucjAPD8O9sB995MKNdmIfrI8qQg2CdpjIyR\n6y4eTFbr9Ia4h0DWbHWr6e7z3MTfXrPHr3zb3uKw7/HKezuoC7NEl9U6PXxnTfA6H1RFKKApRAeJ\n6BcE+ySlkbGLN3e9lz+8urHhdSivNZlpCIKQ7Dja+FdKtQJeBjoBJcD3tNZFAXVuAOYC1cCDWutF\nVu2UUuOAxzx139Na/8pzjX8DeZ7zp7XWlzjpr5e8nFYUlVTQIceeLH5+u1YUWsjNfLzVes9l54Ho\n6qAJzc83p/TjjQ93N3c3BKHF4HQmczOwSWt9LvAScL+xUCnVGbgNGA/MAH6nlEoP0e4vwFVa68nA\nWKXUCM/5gVrryVrr6Y01MODOu3L7t4b7ybhEsnk+wmKDHoIzVAqJTc9O9gN8BSGZcWpkJgGLPa/f\nAc4LKB8DrNRa12itS4ACYIRJu68ppbKBDK31Hs/5d4HzlFKdgHZKqbeUUh8ppRptZHLaZgR5dF1z\nwaCw7YpCRHrHip9cfXaTv6cQnsAgTEEQQhN2uUwpdR1wBzTsYruAI4B3LagUCFR4zDGUA5QBuUB2\nQDvvOeOmRynQF0gH/gA8jnvJbJVSaq3WOqq/8knDu3L0+GlcKfDOx8GquoeLTjVIuwfy5eES0/PR\nYEjv9jG7tuCcTu3NnUUEQTAnrJHRWi8EFhrPKaXewG0c8PwfKEdbgr/hyQaKPecD2wUaKe/5I8BT\nWus6oFAp9RmggJBGJj/f/nLGrVeezfIN+02NzM+eWWv5PqEkZiIhXF/nzjqTZZ/sC7m3E+oaTsYi\nUYnWWHznojP4x7KdMbl2U9HS+htLZCxij9OI/1XAxcAnnv9XBJSvA36jlMoAWgODgS3A6sB2WutS\npVSlUqovsAe4EHgAOB/3vs4lSqksYCiwLVzHCguDRSojobTU3pLYgUP2ZN4H9silIMBYhOvruMH5\njOzfgUf/sZFjJytMAz+trpGfn+14LBKNWI9FSxpn+V74kLHwEUtj63RP5i/AMKXUCmAO8EsApdQd\nSqlLtdZHgfnASuB94D6tdZVVO+Am4G/Ax8CnWuv1WuvFgFZKrcG9j3Ov1jp0dGMjCJeeOZA/vbEp\nfCUDP/3OSL/jnID88FZkpKdy7zWjGD+0i633E2LH0z+ZyuBe7eiRn8X8H05u7u4IQlzjaCajtT4N\nfNvk/B8NrxcACyJstw63J1rg+Tud9M8JOW3C3/QvGd+bRWvcWmRb94QPyDRijBZPS3Vx37WjbLU3\nk8XxZmgUmpa01BTuDnhoEATBHAnGNJDfLnT8zOThXS3Lxg0NFuX04vVo+9pId5rl2785nE4WagN2\n+Pq51ikBBEEQ4gExMgamnt09ZHl+CMMwd+ZQ2rYynxjeePlQAK4+fyAP3TSeYf3sp0zu3tE/d8yV\nnnzzgiAI8YwYGQP9ugZ6YvsTTm3XSlU505MfJsXlcjyDGRzg0jxqUL6j6wiCIDQlYmQMqF6Ni035\nzvnhAzujRZtWzsQ1BUEQmhIxMlHEmO45Fgzv71tmk1wygiC0BMTItCBUr3bN3QVBEARbiJFpQbz/\nyYGG1xnp8qcTBCH+kTuVTS4a2yvo3LwfTGyS9zZG/EvKX0EQWgJiZGxyxTR/1+GRg/L9MiYO759H\nWqoMqyAIAoiRCcnCe6Yza3JfAO666qyG813z2jS8DtyA/+G3hvPXu6bEpD/eYE5BEISWghiZAFpn\numNa+nd3x8xcNrEv8384mTP6+LTNZl80pOF1epr/ELpcLlJcLub9YCKDeuTyyM0Tota30YPzPf93\nito1BUEQYklsfW5bIF3z2rL7UAmtMnxDk9XaPyalQ45veax3F3P10nZZmdxzjT19snCoXu351fVj\n6NKhTfjKgiAIcYDMZAK4avpAslqn840QumAZngh+gDFDrDXLYkGP/CzZ8xEEocUgM5kABvTIDSvf\nbpzZSFCkIAiCNfJI3EjEyAiCIFgjRqaR1NXXN3cXBEEQ4hYxMo0kI02GUBAEwQrZk3HIwzeN53hJ\nBelpqeErC4IgJCliZByS3651yCRmgiAIgkMjo5RqBbwMdAJKgO9prYsC6twAzAWqgQe11otCtVNK\npQKvAs9orZd4zv0cuMRzjTu01uud9FcQBEFoHpxuKNwMbNJanwu8BNxvLFRKdQZuA8YDM4DfKaXS\nrSKppPIAAASySURBVNoppfoBHwKjDdc4GzhXaz0WuBr4s8O+CoIgCM2EUyMzCVjsef0OcF5A+Rhg\npda6RmtdAhQAI0K0ywKuB5YHvMcSAK31fiBVKZWHIAiC0GIIu1ymlLoOuAPw+uq6gCPASc9xKZAT\n0CzHUA5QBuQC2WbttNabPO/lCrjGMZNr+C3LCYIgCPFLWCOjtV4ILDSeU0q9gdtg4Pn/RECzEvwN\nTzZQ7Dkfql3gNYzCYOHqC4IgCHGGU++yVcDFwCee/1cElK8DfqOUygBaA4OBLcDqMO0C3+NhpdSj\nQE/ApbU+HqZfrvx8c8HKZETGwoeMhQ8ZCx8yFrHHqZH5C/CCUmoFUAl8B0ApdQdQoLV+Wyk1H1iJ\ne3ntPq11lVLKtJ2BhvB5rfWnnnprPNe41WFfBUEQhGbCVS+yKIIgCEKMEE0UQRAEIWaIkREEQRBi\nhhgZQRAEIWaIkREEQRBiRkIIZHqCOJ/ErSpQAczRWu9u3l5FB6VUGu44pT5ABvAg8AXwPFAHbNFa\n3+qpG7FenFJqHPCYp+57WutfNeXnagxKqU643eDPA2pJ0rFQSt0DXAak4/7+f0QSjoXnN/IC7t9I\nDXADSfi9UEqNBR7SWk9TSvUnRp/frqZkosxkZgGZWusJwL3AvGbuTzS5Bjjm0XubAfwJ9+e7T2s9\nBUhRSl1uVy8Otxv6VVrrycBYpdSIJv1UDvHcUP4KlHtOJeVYKKWmAOM93/mpQC+SdCxwx9ylaq0n\nAr8GfkuSjYVS6ifAM0Cm51RMPr8TTclEMTINmmha67UYhDYTgH/i+6On4n5SG6m19gayvgOcT+R6\ncV9TSmUDGVrrPZ7z7xKsPxev/AH3l/8Q7vipZB2LC4EtSql/A28Bb5O8Y7EDSPOsaOTifsJOtrHY\nCXzdcDwqBp//fBxoSiaKkQnUSqtRSiXEZ9Nal2utT3n+6K8BP8N9c/Xi1YAz6sKBtV6c91xJwDVy\nY/IBoohS6vvAV1rr9/CNgfHvnDRjAXQERgHfwv0k+grJOxZlQF9gO/AUMJ8k+41ord/E/QDqJVaf\n3+oaliTEjZhgnbMUrXVdc3Um2iilegLLgBe01q/iXmf14tV0s6MXFyhq2lJ04WYD5yulluN+AnsR\nyDeUJ9NYFAHvep5Kd+DeizT+2JNpLO4AFmutFb7vRYahPJnGwkss7hGBdY31LUkUI+PVUsOzWbW5\nebsTPTzrqO8Cd2utX/Cc/kwpda7n9UW4NeDWA5OUUhlKqVyC9eLw/L9Ca10KVCql+nqWGC4ktI5c\nXKC1nqK1nqa1ngZsBK4F3knGscAt2TQDQCnVDWgLLPXs1UByjcVxfE/XJ3A7NH2WpGPh5dMY/S5W\nAxcqpVxKqV5EoCmZEN5lwJu4n3BXeY5nN2dnosy9QDvgfo9XRz3wQ+AJz6bdNuB1rXW9Tb24m4C/\n4X7QWNKCs47eBTyTbGPh8QqarJRah/sz3gzsAZ5NtrHA7QG1UCn1EW5Pu3uADSTnWHiJ2e/Crqak\naJcJgiAIMSNRlssEQRCEOESMjCAIghAzxMgIgiAIMUOMjCAIghAzxMgIgiAIMUOMjCAIghAzxMgI\ngiAIMUOMjCAIghAz/j+VoBUTy6D8NAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118001090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 90.6 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "walk(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"blue\"> Demo, *just-in-time* compilation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import the Numba package, and then add the decorator `@jit` immediately before the function definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `nopython=True` argument activates nopython mode, which means that the code is translated directly to machine code, bypassing the CPython interpreter. The nopython mode is faster than python mode, but is more limited and for example lists and dictionaries are not supported. For optimal performance, try to stick with nopython!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing in Fortran and C "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* One can compile external functions in Fortran or C, and write python wrappers...\n",
    "* but it's even simpler to directly use Fortran or Cython, hiding away the wrapping! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!pip install cython fortran-magic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fortran"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The %%fortran cell magic allows us to write Fortran code into a cell, which gets compiled and imported using `f2py`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%fortran?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The fortranmagic extension is already loaded. To reload it, use:\n",
      "  %reload_ext fortranmagic\n"
     ]
    }
   ],
   "source": [
    "%load_ext fortranmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running...\n",
      "   /Users/ktw/anaconda2/bin/python -m numpy.f2py -m _fortran_magic_67c6d4756b44a5d5ce607e55199e220e -c /Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90\n",
      "running build\n",
      "running config_cc\n",
      "unifing config_cc, config, build_clib, build_ext, build commands --compiler options\n",
      "running config_fc\n",
      "unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\n",
      "running build_src\n",
      "build_src\n",
      "building extension \"_fortran_magic_67c6d4756b44a5d5ce607e55199e220e\" sources\n",
      "f2py options: []\n",
      "f2py:> /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220emodule.c\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7\n",
      "Reading fortran codes...\n",
      "\tReading file '/Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90' (format:free)\n",
      "Post-processing...\n",
      "\tBlock: _fortran_magic_67c6d4756b44a5d5ce607e55199e220e\n",
      "\t\t\tBlock: my_function\n",
      "Post-processing (stage 2)...\n",
      "Building modules...\n",
      "\tBuilding module \"_fortran_magic_67c6d4756b44a5d5ce607e55199e220e\"...\n",
      "\t\tCreating wrapper for Fortran subroutine \"my_function\"(\"my_function\")...\n",
      "\t\tConstructing wrapper function \"my_function\"...\n",
      "\t\t  z = my_function(x,y)\n",
      "\tWrote C/API module \"_fortran_magic_67c6d4756b44a5d5ce607e55199e220e\" to file \"/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220emodule.c\"\n",
      "\tFortran 77 wrappers are saved to \"/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e-f2pywrappers.f\"\n",
      "  adding '/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.c' to sources.\n",
      "  adding '/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7' to include_dirs.\n",
      "copying /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.c -> /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7\n",
      "copying /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/f2py/src/fortranobject.h -> /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7\n",
      "  adding '/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e-f2pywrappers.f' to sources.\n",
      "build_src: building npy-pkg config files\n",
      "running build_ext\n",
      "customize UnixCCompiler\n",
      "customize UnixCCompiler using build_ext\n",
      "customize Gnu95FCompiler\n",
      "Found executable /usr/local/bin/gfortran\n",
      "customize Gnu95FCompiler\n",
      "customize Gnu95FCompiler using build_ext\n",
      "building '_fortran_magic_67c6d4756b44a5d5ce607e55199e220e' extension\n",
      "compiling C sources\n",
      "C compiler: gcc -fno-strict-aliasing -I/Users/ktw/anaconda2/include -arch x86_64 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes\n",
      "\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7\n",
      "compile options: '-I/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7 -I/Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include -I/Users/ktw/anaconda2/include/python2.7 -c'\n",
      "gcc: /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220emodule.c\n",
      "In file included from /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220emodule.c:19:\n",
      "In file included from /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.h:13:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:18:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1777:\n",
      "/Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: \"Using deprecated NumPy API, disable it by \"          \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-W#warnings]\n",
      "#warning \"Using deprecated NumPy API, disable it by \" \\\n",
      " ^\n",
      "1 warning generated.\n",
      "gcc: /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.c\n",
      "In file included from /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.c:2:\n",
      "In file included from /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.h:13:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:18:\n",
      "In file included from /Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1777:\n",
      "/Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: \"Using deprecated NumPy API, disable it by \"          \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-W#warnings]\n",
      "#warning \"Using deprecated NumPy API, disable it by \" \\\n",
      " ^\n",
      "1 warning generated.\n",
      "compiling Fortran sources\n",
      "Fortran f77 compiler: /usr/local/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -arch x86_64 -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /usr/local/bin/gfortran -Wall -g -fno-second-underscore -arch x86_64 -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /usr/local/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -arch x86_64 -fPIC -O3 -funroll-loops\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/Users\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/Users/ktw\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/Users/ktw/.ipython\n",
      "creating /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/Users/ktw/.ipython/fortran\n",
      "compile options: '-I/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7 -I/Users/ktw/anaconda2/lib/python2.7/site-packages/numpy/core/include -I/Users/ktw/anaconda2/include/python2.7 -c'\n",
      "gfortran:f90: /Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90\n",
      "/Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90:1:2:\n",
      "\n",
      " # -vvv for verbose output for what is happening under the hood\n",
      "  1\n",
      "Warning: Illegal preprocessor directive\n",
      "/Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90:2:2:\n",
      "\n",
      " # One can also give compiler flags, e.g. for OpenMP parallellization\n",
      "  1\n",
      "Warning: Illegal preprocessor directive\n",
      "/Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.f90:3:2:\n",
      "\n",
      " #%%fortran -vvv --f90flags='-fopenmp' --extra='-lgomp' # for OpenMP support\n",
      "  1\n",
      "Warning: Illegal preprocessor directive\n",
      "gfortran:f77: /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e-f2pywrappers.f\n",
      "/usr/local/bin/gfortran -Wall -g -arch x86_64 -Wall -g -undefined dynamic_lookup -bundle /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220emodule.o /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/fortranobject.o /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/Users/ktw/.ipython/fortran/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.o /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN/src.macosx-10.6-x86_64-2.7/_fortran_magic_67c6d4756b44a5d5ce607e55199e220e-f2pywrappers.o -L/usr/local/Cellar/gcc/6.2.0/lib/gcc/6/gcc/x86_64-apple-darwin15.6.0/6.2.0 -L/Users/ktw/anaconda2/lib -lgfortran -o ./_fortran_magic_67c6d4756b44a5d5ce607e55199e220e.so\n",
      "Removing build directory /var/folders/2y/083kpwjs5zlc_3db05ll22180000gn/T/tmp9PZPHN\n",
      "\n",
      "Ok. The following fortran objects are ready to use: my_function\n"
     ]
    }
   ],
   "source": [
    "%%fortran -vvv\n",
    "# -vvv for verbose output for what is happening under the hood\n",
    "# One can also give compiler flags, e.g. for OpenMP parallellization\n",
    "#%%fortran -vvv --f90flags='-fopenmp' --extra='-lgomp' # for OpenMP support\n",
    "subroutine my_function(x, y, z)\n",
    "    real, intent(in) :: x(:), y(:)\n",
    "    real, intent(out) :: z(size(x))\n",
    "    ! using vector operations  \n",
    "    z(:) = sin(x(:) + y(:))\n",
    "end subroutine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x = np.random.normal(size=100)\n",
    "y = np.random.normal(size=100)\n",
    "z = my_function(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cython is a superset of Python which supports calling C functions and declaring C types on variables and class attributes.\n",
    "Cython allows you to\n",
    "- wrap C/C++ libraries into Python\n",
    "- optimize Python code by statically compiling with C\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Cython extension is already loaded. To reload it, use:\n",
      "  %reload_ext Cython\n"
     ]
    }
   ],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"blue\"> Demo: the Mandelbrot fractal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "size = 200\n",
    "iterations = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pure python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mandelbrot_python(m, size, iterations):\n",
    "    for i in range(size):\n",
    "        for j in range(size):\n",
    "            c = -2 + 3./size*j + 1j*(1.5-3./size*i)\n",
    "            z = 0\n",
    "            for n in range(iterations):\n",
    "                if np.abs(z) <= 10:\n",
    "                    z = z*z + c\n",
    "                    m[i, j] = n\n",
    "                else:\n",
    "                    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 1: 2.35 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit -n1 -r1 \n",
    "m = np.zeros((size, size))\n",
    "mandelbrot_python(m, size, iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First cython attempt\n",
    "\n",
    "First just add the cython magic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Second attempt\n",
    "\n",
    "Now add type information, use *memory views* for NumPy arrays, and replace `np.dot()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color=\"red\"> *Exercise: calculate primes*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Take the following Python code which returns prime numbers\n",
    "- Run it for `kmax=100` and time it\n",
    "- Add the most simple cythonization\n",
    "- Run again and time, compare to pure Python\n",
    "- Cythonize it properly! (Hint: replace `p=[]` also)\n",
    "- Try just-in-time compilation using `Numba` and compare to `Cython`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def primes(kmax):  \n",
    "    p = []\n",
    "    result = []  \n",
    "    if kmax > 1000:\n",
    "        kmax = 1000\n",
    "    k = 0\n",
    "    n = 2\n",
    "    while k < kmax:\n",
    "        i = 0\n",
    "        while i < k and n % p[i] != 0:\n",
    "            i = i + 1\n",
    "        if i == k:\n",
    "            p.append(n)\n",
    "            k = k + 1\n",
    "            result.append(n)\n",
    "        n = n + 1\n",
    "    return result"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}