Lecture-6A-Fortran-and-C.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Fortran and C code with Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J.R. Johansson (jrjohansson at gmail.com)\n",
    "\n",
    "The latest version of this [IPython notebook](http://ipython.org/notebook.html) lecture is available at [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).\n",
    "\n",
    "The other notebooks in this lecture series are indexed at [http://jrjohansson.github.io](http://jrjohansson.github.io)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The advantage of Python is that it is flexible and easy to program. The time it takes to setup a new calulation is therefore short. But for certain types of calculations Python (and any other interpreted language) can be very slow. It is particularly iterations over large arrays that is difficult to do efficiently.\n",
    "\n",
    "Such calculations may be implemented in a compiled language such as C or Fortran. In Python it is relatively easy to call out to libraries with compiled C or Fortran code. In this lecture we will look at how to do that.\n",
    "\n",
    "But before we go ahead and work on optimizing anything, it is always worthwhile to ask.... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAEgCAYAAABchszxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdclXX7B/DP9xw4bMUByhFEwgUOVBy5QkX0lyMtZ5qr\n1BwVj6OsnifBtByplZapmJjkyp0jUBFHuMUJiQuVIYgsRTbn+v1xd249MkQ4cESv9+t1Xh3u9b3u\nmxPn8jsFEREYY4wxxl4hCkMHwBhjjDFW0TgBYowxxtgrhxMgxhhjjL1yOAFijLEKkp6ejgcPHhg6\nDMYYOAFijL2AsrOz8ejRI0OHoTcajQZxcXFo0aIFXF1dcefOHQBAcnIy5syZg0uXLhk4QsZePZwA\nMfYKS0hIgJeXF65du2boUAAAWVlZiI+PR//+/eHo6IgtW7YYNJ4bN26gTp06MDc3x4EDB0p1jXnz\n5sHIyAj29vaYOHEivLy8UK9ePXzwwQfYs2cPZs6cidGjR+s3cMbYM3ECxFgFe/jwIQYMGIDmzZsj\nPDxcZ9+6deugUCgKvNRqNQICAp67rJ07d0KhUGDjxo3YvHkzPv30U3kfEWHgwIEIDg6WE434+HiY\nmZlh4sSJSExMRPXq1QuNx9zcHGfPni2yrKioKCiVSqxcubJATKNHj4ZCocDRo0cBAJmZmZg1axaU\nSiXMzc2hVqvx4MEDPHr0CLt37y7RfXbq1AmfffZZofsyMjIwefJknD9/HnFxcWjatGmhxz148AC3\nbt2Sf46Pj0ePHj2QlpaGvLw8eHt7IyMjo0TxPOmzzz7D3LlzsX//fkybNg1Tp04F8Pj516hRA+fP\nn39hklDGXhVGhg6AsVdJZGQkRo0ahVOnTgEAunbtiv3798PNzQ0A5C9BR0dHTJgwATY2NsjJycG3\n336LUaNGwdbWFl26dEFKSkqx5VhaWsLS0hIHDx6Uy01JScGSJUugVCoxb948/Pzzzzh16hQUCgW2\nbNmCL774AsePH0d2djYiIyNhbW2NpUuXIjs7W75ubm4u5syZg9jYWCxevBjr1q2T92nLunr1Klq2\nbAkiwv/+9z94eHigUaNGBWJ0dHRETk4O+vbti4MHD6JOnTqYPHkyWrVqhUaNGuH+/ftwcnIq0XP1\n8vKCv78/FixYUGDfb7/9hlWrVuHdd9+FtbV1kU1rYWFhmDt3LoKCggAAs2bNQlJSEq5cuYKTJ09i\n0KBBWLp0KWbMmAEASElJ0Xk2hbGxsYFSqZTPAYDVq1cDAIyNjWFmZobJkyfj66+/xvHjx9GgQYMS\n3S9jTA+IMVYhkpKSqGXLltSiRQvatGkT+fj4kLW1NVWvXp1Onz5NREQ+Pj4khKCTJ0/qnPvHH3+Q\nEIJ++eUXWr58OQkhSAhBZmZm8vsnX3369CEior1795IQgnx9fek///kPCSGoVq1aFBYWRmq1mhYv\nXkzz5s0jIQQREWVkZJAQgrp27VroPdy8eZMcHBzI2tqazp49q7NPW9asWbPoypUrciy9evXSOW7U\nqFGkVCopISGB1q5dS0IIeu+99ygnJ6fUzzYoKIiEEHTmzJkC+0aNGkVNmjQhIqL9+/eTk5NTodcI\nCQmhLl26yD+7ubnR4MGD5Z/d3d3Jzs5O/vn1118nIQQZGRmRsbFxob+HLVu2yMdnZmbK92tsbEyn\nTp0iIiJ/f38SQtBvv/1W6vtnjD0/bgJjrIJ8/vnniI+Px9GjRzF48GD4+vri2rVrqFmzJoYPH478\n/HwAUq1B27Zt5fNycnKwZMkS2NnZYdiwYfjwww9x7do1hIaGIi0tDY0bN0aDBg2we/duBAYGYv/+\n/di+fTsAYMeOHfJ16tatCwC4d+8e3N3dUbt2bUyYMAG1a9eWj3ny+KcdOXIETZo0QUpKCvbs2YNW\nrVrp7H/y3NzcXPl9YGAggoOD5Z/z8vLQoUMH2Nra4saNG3LznrGxMTQaDeLj4xEfHy8/j5Jo2LAh\nAMi1O3/99RcyMzORlZWFHTt2oFevXgCAv//+Gz169Hjm9WJiYnDx4kU4ODjI20aOHIn4+HicPHkS\nAHD8+HEcOXIEN27ckGu/xo8fj8DAQAQGBuL8+fMYMGAAAOD3339Ho0aNMGrUKCiVSqxZswZt2rQB\nAFhZWQEAqlSpUuL7ZYyVHSdAjFWA+Ph4BAQE4Msvv4SlpaW8vWbNmhg3bhyuXbuGbdu2AQASExMx\nc+ZMxMXF4fTp03jzzTdx4sQJrFmzRv6SdHZ2Rvv27WFsbIxatWqhY8eO6NWrF3r06AFPT08YGUmt\n29ommqpVq8oJi7e3N0xMTLB8+XKYmZmhY8eOcjza45/+Ml68eDF69+4Ne3t77Nq1Cx06dChwj0+e\nu3v3blhZWeHChQtwc3ODt7c30tPTAQDbtm2T46tXrx7S09ORmpoKANi6dSvUajXUajVcXFwwYsQI\naDSaZz5fc3NzKJVKhIeHY9euXejduzdmz56NkydP4sGDB7C3twcACCFga2srn5ecnIz4+HicO3cO\nCxYsQExMDACpOQwAPD09AUj9dfr27QtAak7U6tSpE+rWrStff9iwYejRowd69OiB5s2bAwDOnTuH\nkSNHIjo6Gg4ODggLC8OwYcPkawwYMABbtmxB//79n3mfjDH94QSIsQpw9+5dAMA777xTYN/06dNh\nZWUFFxcXeducOXNgb2+Pdu3aITQ0FBs2bICXl1epy+/fvz/o32X/pkyZgtDQULkGQqEo+GfgyS/j\nxYsXY/r06ahTpw5CQ0PRpUuXZ5aVlZWF6tWro1mzZliyZAmuXr2Kfv36ITw8HFlZWRBCAACGDh0K\ne3t7fPnllwCkZODAgQMIDAxEnz59sG7dOvz555/PvD9bW1u0atUKhw4dgp+fHwAgNDQU//zzDxQK\nBd5++20AUiITFBSEuLg4nDt3DrVq1YJarYa7uzsCAwPlGh9tAnTq1Clcv34dnp6ecHZ2BiAls88j\nLS0NgJR8/f7772jWrFmBYwr7XDDGyhcnQIxVgD179qB27dpQq9WF7lcoFHINCQC4urpi7dq1CAwM\nRFhYGAYOHFiqcpOSkuTrP+nJ5ivtaKzCjvfz88P06dNRr149BAcHw8bGptiyhBBycqPVqVMnBAQE\n4PDhw/KXf7du3QAApqamWLx4MVatWoXAwEAoFAp069YNPXr0wIIFC2BmZobo6OgS3auFhQU2b96M\n3bt3w9jYGEePHsWSJUvQsmVLnaas06dPw97eHu7u7tBoNJgyZQr++usvucbmSbNmzULDhg1hZGSE\nrVu3QgiBQ4cOlSgeLUdHRxgZGYGIcPny5ec6lzFWfngUGGMVwNHRETExMThw4AC6d++us2/VqlVI\nS0uTawoAwN/fX66hKYujR48WSEiepq0Z0h7/pLlz5wIAli5dijp16iA9PV0nUQOkPiwWFhYFzn3S\nkCFDsG/fPvj7+0MIoZPQ9ejRA++88w5mzJgBDw8PmJmZAQA+/vhjZGVloWvXriW61z59+uDQoUNo\n0KABQkJC8Nprr+HKlSuYPn26znE1a9bEjBkz0KxZMwgh4OHhAZVKhfXr1+skW0ZGRvDz80Pbtm3l\n2rkuXbrgzJkzyMvLk5vxnsXJyUlOsgprOgQezzdkampaomsyxsqOEyDGKsCQIUOwevVqTJo0CZ9+\n+qmclISFhcHf3x/e3t7o2bMnjh8/DrVa/VzJz5MJzJMiIiLw4MEDCCFgbGwMa2vrIq+hrfE5duwY\nAGmINiB1yL516xaGDh2Kvn37IjAwUO6vozV06FB89dVXOmUV5qeffsKDBw/wxhtvFBgWv2rVKnTp\n0gWtWrXC4sWLceTIEaxYsQI+Pj5FztvzNBMTE5iYmGDJkiVQq9UYOnQoAgICMHz4cPkYIQQmTZqE\nadOmFTj/6URRrVZj1KhROtsaNWqEkJAQZGRkPFen5fnz52P+/PmF7tu+fTsGDBiAWrVqYefOnTod\n4Blj5YcTIMYqgEqlwt69ezFp0iRMnDhR7tjbokUL/Pbbbxg8eLDOsSWVn5+PuLg4vPbaawX23bt3\nDwBQp04d2NnZwc7ODvXr15drWLRef/11eZLF+/fvQwiBdu3aAZDm0BkxYgTCwsKwadMmAJCTnI8/\n/hheXl5wd3eXl3KoU6dOkc18ZmZm2Lx5c6H7LC0t8ddff2HcuHHo3bs3TExMsHDhQnnSwJJ4//33\n0blzZ7kpS5twaudY0iqqRqxdu3aoWrUqAGkuI20z3ZPeeustLF++vMD2kjbTFaZq1apQKpVQKBRy\n+Yyx8scJEGMVxMzMDP7+/vD390dOTg7OnTsnJxpaXbp00RlC/ixZWVl48OBBoV/WWk/WyFy9erXA\n/saNG6Nx48aFntO4cWOcPn26xPFoz3tWs1thbGxssH37doSEhKBWrVpo0qTJc51vZmam049HoVAU\nqE0xNzeXpwN42sSJE+X3UVFRaN++fYFjevTogTZt2sDc3Fxne1xcHGxsbNCiRYtCr3379m04OjoW\nuq9bt27P9TtnjOmHoKLqzxljldrx48fRuXNnLFq0CN7e3iU6x9TUFG+99RY2bdr0XEnM02WdOnUK\nM2fORGBgYGnDf2kcPXoUHh4e+OWXX/Dhhx8aOhzG2L+4BqiUEhIS8M0338De3r7INYgYM6T27dsj\nLy/vuc7JysrSS1lt27bl5OdfO3bsgEKhQJ06dQwdCmPsCS9FDRARlarKvSyWLFkCb29vGBkZITY2\nVmdyNcYY03rw4AFu3LiBli1bGjoUxtgTXooESEuj0UChUMijYsorKcrLy0OzZs1w5coVAEDv3r2x\na9euCk/CGGOMMVY6lToB+vvvv3H58mXcuHEDAODi4oK+ffsWO1mbPqxatQrjxo2DEEJOtn7++WdM\nmjSpXMtljDHGmH5U2gTozJkzOiM8zMzMkJmZiZo1a+LLL7+El5dXsaNISttslp6ejgYNGiA+Ph6N\nGjWS1wUyNTVFWFiYznIGjDHGGHsxVdpO0CqVCj4+PmjRogVUKhXy8/MRExODffv2wcfHB7t27cKy\nZcvQoEEDnen5g4ODERcXhxEjRpSq3EWLFiE+Ph5t2rRBp06dEBkZCTc3N1y4cAHDhw/HiRMnnmse\nF8YYY4xVvEpbAwQUXouTnp6OQ4cO4YcffkBKSgoWL14MDw8Peer6QYMGISMjA5s3by4wl8ez3L17\nF/Xr10dGRgYOHz6M+Ph4DBkyBF5eXrh+/TqioqLw2WefFTnjK2OMMcZeDJV6MVQhBPLy8pCXl4f8\n/HwA0oyyffr0wYIFC2BlZYWxY8fi+vXr8ro9p0+fhoODg9x35+lhv8UNA/bx8UFGRgb69euHN954\nQ26CCwsLQ0BAABQKBb777juEhISUx+0yxhhjTE8qdQIESAsWGhkZQalUApBGgmk0GrRq1QqHDh1C\nlSpVsHjxYuTk5CAvLw937tyBu7s7LCwskJubi++//15OeogI8+fPR2JiYoFyoqOjsXr1aiiVSrmG\nx9HRETY2NkhKSoKdnR2++uorEBHGjh373POvMMYYY6ziVNoEKCYmBn5+fjhy5Aju378vJzEKhUJe\n2BEAhg8fjk2bNkGlUuHkyZMQQsjrJoWGhmLu3LnyCsxXr17FrFmzCl3Xx8HBAUeOHMGiRYvkhRyf\nXDPp1KlT+N///oeRI0diy5YtJV4pmjHGGGMVr9J+S+/evRuTJk2CmZkZWrdujRYtWsDNzQ2NGzeG\nlZUVHj58iHPnzuG7776TR4MdOnQI9vb28lpAp06dQrVq1XDv3j3Y2toiLCwMVapUgYWFRaFldujQ\nAR06dNDZ5u7ujt27dyMsLAxDhw7Fb7/9Vr43zhhjjLEyq5QJkEajwYQJE/B///d/uHTpErZt24Yd\nO3Zg6dKlMDY2hrm5OdLS0gAAAwYMwJw5cwAAZ8+eRXJyMkxMTAAAR44cgZeXF6pXrw5AqhHKz8/H\nvXv30KhRoxINldfO7nru3Lnyul3GGGOM6VmlTIAUCgXy8/Oxfv16mJmZ4ddff4VCoUBOTg5Onz6N\nQ4cOITk5GYMHD0bLli3lYenVqlXDo0ePMGjQILRp0wYHDx7Erl275OaquLg41K5dGzVr1gRQsrmC\nOAFijDHGKp9KOQz+/PnzeO+99xAREQFTU1N4enpi5cqVsLOzK/a8mzdvYvDgwbC2tsadO3fw4MED\nKJVKjBgxAg0bNsSCBQtQv359rF+/HlWqVClRLEQEU1NT5OTkIDMzU+5PxBhjjLEXV6VLgE6fPo0P\nP/wQubm5WLhwIZKTk/HJJ5/gzTffxNq1a5GTkyN3hH6yM7R2mLxCoYAQAo8ePUJoaCgWLlyI2NhY\n/PPPPwCALl264ODBg88Vk52dHeLj4xEbGwu1Wq2/m2WMMcZYuah0TWA///wzoqOjsW7dOvTo0QMA\nkJCQgK+//hoHDhxA9+7dCzRd3bp1CwEBATh//jzu3LkDOzs7dO/eHd26dcO+ffuQlZWF0NBQbNiw\nQa750U6cWBLVq1dHfHw8kpOTOQFijDHGKoFKNwz++PHjeOedd9CpUyd526BBg6BSqRAcHAyg4Crw\nN27cwLx587B9+3bY2toiPT0dy5cvR4cOHaBUKtGhQwds3LgRH330ERYvXgwAzzWMXduJOjk5uay3\nxxhjjLEKUOkSoJs3b0Kj0UClUkGj0YCIkJubi4yMDNSuXVtu6nqSp6cnJkyYIL/fvXs3IiIikJqa\nit9++w1t2rTB8ePHsWnTJjx48OC5Y+IEiDHGGKtcKlUT2Pnz55Gfnw8ikmtoNBoNfHx8YGVlhb59\n+8ozQj9t0aJFqFq1KhYvXgwzMzNMmDABCoUCw4cPx4ABA5CWlgYhRIk7Pz+pRo0aAICkpKTS3xxj\njDHGKkylSoC0tTPh4eEICAhAZGQkbt++je3bt2PIkCHyDM9FmTlzJvLy8jB58mRYWFhgxIgREELA\nzMwMZmZmpY6La4AYY4yxyqVSJUCtW7fG9OnT8csvv2D8+PHIzs4GAFhYWMDf3x9HjhxB06ZN0bJl\nS4wcORL16tUDICVOf/75J9atW4fw8HAIIRAREYH8/HwYGRlBo9EAgM6osefBCRBjjDFWuVSqBMjc\n3BwLFizAggULEBkZifDwcKSlpSE2NhZXrlzBtWvXcOnSJezatQuOjo5yArRu3Tp89dVXcoJib2+P\nWrVq4erVq7C3ty9Vs9eTOAFijDHGKpdKkwARESIiIpCeno527dqhUaNG8qKkgNQX6MKFC5g3bx7+\n+OMPNGvWTN7n6uqKiRMnIj4+HhqNBsnJyZg7dy6mTZsGS0tLNGrUCJ06dcLXX38NKyur545NmwBx\nHyDGGGOscqg0CdDFixfx1ltvITs7G6GhoXB2dkZeXh6EEFAqlVAoFIiKisLmzZvxzTffwNjYGAAQ\nFRWFtWvXombNmmjcuDFq1qyJ1q1bw9bWFjdv3sSFCxcQERGBR48elbofkLYTdLnWAIWEAGvWAP7+\nQCmb6hhjjDEmqTQJ0Pr16xEdHY2dO3fC2dkZQMG5epo2bYrmzZtjzJgx2Lt3L6ysrBAeHg5/f3+o\nVCrk5eWhbt26cHd3x9ixY1G/fn2MHDmyzMtXlHsTWE4OMHIkEBMDeHpK7xljjDFWapWmKuHOnTuo\nUaMGHBwcAKDQ+X4aNmyI6dOnIzQ0FNu2bQMgzfuzYcMG+Pn5Yfr06cjLy8OuXbvw5ptvomHDhmjV\nqhU+/vhj7Nq1C/fv3y9VbOWeAKlUwDffSO+//BLIyCifchhjjLFXRKVJgNq3b4+kpCRs3LgRAKBU\nKuXRW0969913Ubt2bWiXODMzM8OQIUOwe/duDB06FHv27EFAQAC+/PJL1KlTB1euXMHPP/+Mfv36\noX///rh9+/Zzx6ZtAktMTES5La323ntAy5ZAbCzw/fflUwZjjDH2iqg0i6HeunULHTt2xN27d7F5\n82YMGDAAAOR1v7T/vXbtGpo0aYKIiAjUr19fPve1117DiRMn0LZtW/ma+fn52LZtGw4cOICIiAiE\nhobi4sWLaNq06XPHZ2VlhfT0dKSkpMDa2lo/N/20gwelJjBLS+D6daBWrfIphzHGGHvJVZoaoHr1\n6mHu3LlQqVTw9vaGv78/cnJy5HW/hBBITk7GqFGjYGFhISc/AHDu3DmYmZnBxMQEGo0GOTk50Gg0\nUCqVGDRoEBYvXow//vgD//zzT6mSH0BaER4A7t69W/abLUq3bkCfPkB6OuDrW37lMMYYYy+5SpMA\nAcDIkSOxZs0apKWl4YMPPsBXX32FK1euIDw8HIsWLcLrr7+OEydO4NNPPwUAuTnq2LFjcHZ2houL\nCxQKBVQqFRQKBYgIGo0GFhYWsLOz0xlW/7y0q8DHxcWV/UaLs2ABoFQCK1cCERHlWxZjjDH2kqo0\no8AAKaEZOnQo6tevj82bN2Pv3r1YsWIF0tPTodFooFar8eOPP2Ls2LE6512+fBmRkZGYO3cu2rdv\nj7p168LR0RFmZmYFVo4vrQqpAQIAFxdg/Hjgl1+Azz4Ddu8u3/IYY4yxl1Cl6QNUmBs3biAsLAzG\nxsawsbGBra0tGjRooHNMTk4OHBwckJWVBbVaDRMTE1hZWcHOzg5OTk5wcnKCp6cn6tevX6ZkaNq0\naVi8eDHmz5+Pzz77rKy3Vrx794D69YGHD4EDB6R+QYwxxhgrsUpVA/Q0Z2dneU6goly7dg2ZmZmI\ni4vD5cuXcfHiRVy6dAnXrl3DgQMHcPfuXWzYsAF79+6FhYVFqWOpsBogALC1Bb74QhoSP2UKEBYG\nGFXqXyVjjDFWoV7ab03tqLDQ0FCYmJgAAF5//XW8/vrr8jFxcXE4cuQIQkNDy5T8AI8ToHLvA6T1\nn/9I/YAuXQJWrAAmT66YchljjLGXQKXqBF0aly9fRpcuXWBpaQkiQl5enjx/kFqtxtChQ7F06dIy\nl6PtBF0hNUAAYGYGLF4svf/qK4DXIWOMMcZK7KVNgLT9ecaNG4eFCxfK242MjKD4dy0tIkJ+fr5e\nJi+s8BogAOjfH+jeHUhJkZIgxhhjjJVIpe4E/SK5f/8+bGxsUK1atfJdFPVpERFA8+YAEXD2LNCi\nRcWVzRhjjFVSL20NUEXTzv6cmppa6BId5cbVFfjoI0CjAT75REqEGGOMMVYsToD0xMjICFWqVAER\nIS0trWIL9/UFatYEjh4F/l0rjTHGGGNF4wRIj5RKJQCU34KoRbG2BubOld5PnQqkplZs+Ywxxlgl\nwwmQHuXm5gJ4nAhVqPffB9q3B+Ljgc8/r/jyGWOMsUqEO0HrSVZWFszMzGBkZITs7Gx5pFmFCg8H\nWrYEcnOBI0eAzp0rPgbGGGOsEuAaID3RDn+3s7MzTPIDAE2aPK79GT8eyM42TByMMcbYC44TID0J\nDw8HgAJrkVW4L78EGjUCrlwBvv3WsLEwxhhjLyhOgPTk3LlzAICWLVsaNhBTU2mJDEBKgM6cMWw8\njDHG2AuIEyA9OXz4MACgbdu2Bo4EwBtvSHMC5eUB770HZGQYOiLGGGPshcKdoPUgMzMT1apVQ3Z2\nNu7duwcbGxtDhwRkZgKtW0szRU+aBPz8s6EjYowxxl4YXAOkB8eOHUN2djZatGjxYiQ/gLRY6rp1\ngLExsGwZsHevoSNijDHGXhicAOnBgQMHAACenp4GjuQpLVoAc+ZI799/H0hMNGw8jDHG2AuCEyA9\n2L9/PwDAy8vLwJEUYto0wMMDSEgAxo3jtcIYY4wxcB+gMktKSoKNjQ2MjY2RkpICc3NzQ4dU0J07\n0orxaWnSCLFx4wwdEWOMMWZQXANURgcPHgQRoVOnTi9m8gMAdetK/YAA4D//Aa5dM2w8jDHGmIFx\nAlRGx44dAwB07drVwJE8w7BhwLvvSkPihw+XlstgjDHGXlGcAJWRdgLEVq1aGTiSEli2DHBwAE6f\nBmbPNnQ0jDHGmMFwH6AysrOzQ3x8PG7fvo26desaOpxnO3wY6NoVEAI4ehTo0MHQETHGGGMVjmuA\nyiA3NxcJCQkQQsDOzs7Q4ZSMhwfw2WeARgN88AEvmMoYY+yVxAlQGaSmpoKIUK1aNRgbGxs6nJLz\n9QUaNpQWTF2wwNDRMMYYYxWOE6AyyM/PBwAYGRkZOJLnZGoKLF8uvf/mG+DqVcPGwxhjjFUwToDK\nQKlUAgDy8vIMHEkpdO0KjBolNYFNnMgTJDLGGHulcCfoMsjLy4NKpQIA5OTkVL6aoPv3gcaNgaQk\nICBAWjmeMcYYewVwDVAZGBkZwdbWFkSEu3fvGjqc51ezJvDdd9L7qVOBlBTDxsMYY4xVEE6AysjF\nxQUAcPnyZQNHUkqjRwNvvCEtlPrFF4aOhjHGGKsQnACVUbNmzQAAly5dMnAkpSQE8MsvgLExsGIF\ncPy4oSNijDHGyh0nQGVU6RMgAHB1BT79VHo/fjyQlWXYeBhjjLFyxglQGWkToIsXLxo4kjL673+B\n+vWBy5eBGTMMHQ1jjDFWrngUWBmlp6fDysoKxsbGePToUeWaEPFpZ84A7dsDeXnAn38CffsaOiLG\nGGOsXHANUBlZWlrC2dkZubm5iIyMNHQ4ZdO6NTB3rvR+zBggKsqw8TDGGGPlhBMgPWjevDmAl6AZ\nDJCGw/fqJc0N1Lcv8OCBoSNijDHG9I4TID3Q9gOqtEPhn6RQAOvXAy4uQHg48O67wL9LfjDGGGMv\nC06A9OC1114DANy5c8fAkehJ1arArl1A9erA3r3AtGmGjogxxhjTK06A9ECtVgMAYmNjDRyJHjk7\nA9u2SfMD/fij9GKMMcZeEpwA6YE2AYqLizNwJHrm4QH4+0vvp0wBtm83bDyMMcaYnnACpAdVqlQB\nIA2Jf+kMHw7MmSOtFj9sGHDihKEjYowxxsqMEyA9MDExASCtCP9S+vJLYOxYaYbovn2BGzcMHRFj\njDFWJpwA6YE2AcrOzjZwJOVECGDZMqBnT+D+feDNN6Vh8owxxlglxQmQHqhUKgAvcQIESJ2hN28G\n3NyAa9eAfv14zTDGGGOVFidAevBkE9hLvbKIlRWwZw9gbw+EhgIjRwIajaGjYowxxp4bJ0B6oFAo\nYGRkBAD9okB1AAAgAElEQVTIzc01cDTlrE4daW6gKlWkGqHPPzd0RIwxxthz4wRIT7TNYC9tR+gn\nNWsGbN0KGBkB330n9Q9ijDHGKhFOgPREuwr8S18DpNW9O+DnJ73/+GPgwAHDxsMYY4w9B06A9EQI\nAQAvdx+gp40eLQ2R12iAwYOB69cNHRFjjDFWIpwA6YlCIT3KVyoBAoDZs6W5gVJSpJFhDx8aOiLG\nGGPsmTgB0hNtDZDmVRsVpVAAv/8urR4fEQGMGMEjwxhjjL3wOAHSk1eyCUyrShVg507A2lr6r6+v\noSNijDHGisUJkJ7k5eUBAJRKpYEjMZAGDYBNm6QaodmzpSHyjDHG2AtK0CtZZaFfRAQjIyNoNBrk\n5ubKcwK9khYvBqZNA8zNgWPHpJmjGWOMsRcM1wDpQXp6OjQaDSwsLF7t5AcApkyR+gFlZEidou/f\nN3REjDHGWAGcAOlBQkICAKBGjRoGjuQFIASwciXQpg1w+zbw9ttAZqaho2KMMcZ0cAKkB1evXgUA\nNGjQwMCRvCBMTYHt2wG1Gvj7b2DoUODfPlKMMcbYi4ATID2IjIwEADRs2NDAkbxA6tQB9u0DqlUD\n/vwTGDcO4O5mjDHGXhCcAOmBtgaoUaNGBo7kBdOkibRwqrk5sGYN8OmnnAQxxhh7IXACpAfaGiBO\ngArx+uvAtm2AsTGwaBHw1VecBDHGGDM4ToD0gPsAPUPPnsC6dYBSCXzzjbR+GCdBjDHGDIjnASqj\nrKwsmJmZQalUIisri4fBF2fLFuDdd6UO0dOnAwsWSKPGGGOMsQrGNUBldOvWLQCAo6MjJz/PMnAg\n8McfgJERsHAhMHYsjw5jjDFmEJwAldHt27cBAPXq1TNsIJXF228DO3YAZmbA6tXAO+9IkyYyxhhj\nFYgToDJ69OgRAKBq1aoGjqQS6d0bCA4GqlcHdu0CuncHkpMNFs7du9L6rf9W5jHGGHsFcAJURtnZ\n2QAAlUpl4EgqmfbtpUkSHRyA48eBTp0MloHMmAHMmgUsW2aQ4hljjBkAJ0BlZG5uDuBxTRB7Di4u\n0oKpTZoA//wDtG0r/VzBPv5Y+q+fH8C/Rsaez31e769Syc7ONsj31Yv4OeEEqIyqVasGAEhKSjJw\nJJWUvb1UE+TlBSQmAl27Ar//XqEhtGkjTVeUmiqN1meMFY+IcO/ePaxZswa2trYYM2YM9DmgODIy\nEuPHj0d4eLjerllWsbGx8Pb2RkhISKmvUdL70mg0GDZsGBQKBcaOHVvq8p6UlZWF+Ph49O/fH46O\njtiyZUupr7Vz506oVCo4OTkhNja2yOPK+3NSZsTKJCoqigBQ7dq1DR1K5ZabSzRpEpE0QxDRf/9L\nlJ9fYcWvXy8V26QJkUZTYcUyVuFiY2OpSZMmhe5LS0ujqKgoIiLKyMggtVpN//d//yfvDw4OJnt7\nexJCyK+uXbuSpaUlPXr0iIiIevbsSSdOnCj0+ufOnStRjJ6eniSEoKVLlxZ5zLx588jKykonFu3L\nwsKCYmJi6MKFC4XGMnXqVPLx8SEiotWrV5OtrW2h11EqlRQWFkZERB988AEJIWjatGkluofS3hcR\n0aRJk0gIQVWrViUhBB08eFDed/PmTVIoFLRixYoC540aNYqEEHTkyBEikn6Hvr6+pFAo5Hvq0KED\nmZqa0qhRo4qNYcKECTR58uQC2w8fPkympqZybCNHjixwTEk+JydOnKBmzZoV+tyFEPTDDz8UG58+\ncAJURnl5eWRiYkIA6OHDh4YOp/JbupRIoZCykYEDif79n6W8ZWcT2dlJxQYHV0iRjBlEVFQU1atX\nr9B9ISEh1KNHDyIiCgsLIyEEOTk5ERHRqlWrSAhBZmZmNGzYMNqzZw8dO3aM4uLi6OLFi/I1unTp\nQocPHy70+k5OTnT37t1i40tNTaWqVauSQqGg2NjYQo85efKknOgsWrSIfv31V51XaGgoERGNHj2a\nWrRoUeD8N954g8aMGUMJCQkkhCAjIyPy9fUtcJ2goCD5HEdHR1IoFPK1n1dJ7otISgyEEPTjjz9S\nYmIiOTk5kaurK2n+/ZfZP//8Q0IIsrGxoStXruicq02Abt++TdnZ2XLCZW9vT3PnzqWgoCC6desW\nnTlzhpKSkuTzcnNzKSsrS+dakyZNIhMTE53fl0ajIVdXV/Lw8KD09HT66KOPSAihk2SW9HPi4uJC\nQgh67733Cjz3gIAA+X7LEydAeuDq6koASvyvG/YMgYFEVapI2Yi7O1Exfyz06euvpSL79auQ4hgz\niP3798tJzdNCQkKoS5cu8s+2trbk5OREDx48IGtra3JycqLw8PAir52dnU3Ozs5yDcTT6tWrR7dv\n3y42vqioKDkpKUpOTg5ZW1vTzJkzi73W6NGjSQhBiYmJ8rb79++TqakpLVy4kIiImjZtWmgtxtOE\nEKRQKCgmJuaZxxamJPdFRPT999+TiYkJpaSkEBHR0qVLSQhBe/fuJaLHCZAQgnr16qVz7qhRo0ip\nVFJCQgKtXbtWTjBycnKKLbN///7k5uams+3ChQukUqno888/l7elpKSQEIKWLVtGRNKzNDExocGD\nBxORVINYks8JEdHkyZOLTMQrCvcB0gPtEhjXrl0zcCQviZ49pZFhTk7A2bNSJ50zZ8q92PHjAZVK\nWrw+Kqrci2Ps2TIz9X7Jv//+Gz169HjmccePH0diYiIAIDk5GWlpaVi5ciVcXV3lbfHx8UhPT5fP\niYuLQ0pKCtzd3Usdn0ajAQC0a9euyGOMjY1hamqKq1evIj4+Xn6lpKToHNe5c2cAwOHDh+VtDx8+\nRHZ2Ntq0aQMAsLCwQFRUFOLi4uTrPN1hVxuTWq1GnTp1yu2+AGDPnj2wtLSEtbU1AKBPnz6oVq0a\n/vzzTwBA3hOTxwYGBiI4OFj+OS8vDx06dICtrS1u3LgBtVqNgIAAGBsbQ6PRyPeXn5+vE9elS5cQ\nExOjE0fz5s0xfPhw+Pn5yfv27NkDAHBwcAAA1KhRA7169cLBgweRmZlZ4s8JAFhaWiItLQ1XrlyR\n40pISKjQPkKcAOlB/fr1AXACpFeursDJk9Lw+Lg4oHNnYNOmci2yVi1g6FCpE9KSJeVaFHvVbdwo\ndfjv2hXo1k2aC8vDA3BzAxwdAVtbwNISsLPTe9FCCNja2so/a7+gzp07hwULFshfdtopPqpUqYLq\n1avDysoK0dHR8nkODg5Qq9VQq9Xw9PREWFgYACmh0I6O1Xa8jY2Nxbx58xAdHY3o6GgQEdasWYOL\nFy8iJiYGMTExuHfvHgDg9OnTAIqfWoSIkJOTg61bt8LZ2RlqtRr29vbo2bOnzginkSNHwsLCAgsW\nLJC3hYWFQaVSwcXFBQCQk5ODY8eOoVGjRvL9dOrUSadzrzYmY2Pj533cBa5R3H1lZmYiIiICnp6e\n8raqVauiY8eO8qLbu3fvhpWVFS5cuAA3Nzd4e3vLycW2bdvkFQnq1auH9PR0pKamAgC2bt0q35+L\niwtGjBgBjUaD69ev4+bNm+jXr1+BePr27Yvk5GSMGTMGAHD27FkIIdC1a1cA0u93wIABSEpKQmJi\n4nN9TnJycpCamop27drJ+5s0aYL9+/eX7gGXhkHrn14Sy5cvJwA0ZswYQ4fy8snOJvrgg8edo//3\nv3LtHH3unFSMhQXRE03kjOnXggWPP9PFvUxMiDIy9Fq0j48PtW3blmJjYyksLIyMjIwKdFYlkprD\nhBByZ+FZs2ZR9erV5b4i4eHhFBQURL/++iup1Wrq168f3bp1i5RKJR0+fJju3r1LXbt21bm2QqGg\n27dvk7+/f5EdjlevXk1CCBowYECR93D06FESQlBAQABlZWXR/v376caNG4UeO2zYMBJCyPs9PT3J\nwcGBiIiio6NJCEGzZ88mjUZDwcHBdPny5QLXOHjwIAkhyN3dvdTPvST3pW0ma968OcXExND8+fNJ\nqVSSEIKaNm1KRNLvT9t0dPToUTI2NqZu3brR5cuXSQhB3bp1IyKizMxMatKkCU2cOJGIiPLz8yk4\nOJiCgoJoypQpJISg7du3k5+fX4GO1lopKSmkUqmoZcuWRETk4eFBQgjat28fHTp0iJycnOTf36lT\np4ioZJ8TIiJnZ2fy9PQkIqLTp08X2WxanjgB0oNdu3YRgALtsUxPNBqiH3543Dn67beJyrHDeY8e\nUjGzZ5dbEexVd/s20cGD0is4mGjfPqKQEKKwMKKoKKK7d4lSU8tlSKKPj0+BpGTq1KkUGBhIbm5u\ncgK0ZcsWEkLQrFmziIgoKyuLnJ2d6d133y1wzbFjx1KrVq3o1q1bBRKbZs2aUVBQEE2fPl3uoJuZ\nmUn79++noKAg+RUREUFEjzvy+vv7F3kP2uTsWf2JiB4nL9q+LG5ubvLoJm3CUVSn7aefmTYZLI2S\n3Jc2Hu3LysqKtm7dSp06dSIhBCUkJOgkQEREGzdulJMkIQTNmTNH3hcUFETGxsb0119/6ZSTm5tL\n5ubmtGTJEvneCnsG2s7m2iRKmwBpX9OmTaNZs2aREILmz59PRCX7nBBJ/cG0ny1D4dU79UDbVvt0\n+zPTEyEAb2+gcWNgyBBg+3agY0eps46jo96LmzED2LdPagabNk1atowxvapbV3oZSM2aNTFjxgw0\na9YMQgh4eHhApVJh/fr1cvPF0aNHdc4xMTHB0qVL0adPH0yaNAmdOnUCAISGhmLt2rX4WDujKIAR\nI0Zg6NChUCqVaNCgAZycnGBsbIxFixYBAExNTdG9e/dCY9M22axevRrdu3fXWWRaoVDA1tZW7m5w\n//591K5dG8lPLKWjUqlQvXp1+WcPDw9069YNfn5+mDx5MtLT09G7d28A0LlOfn6+3OcJAJRKJWxs\nbHRi2rRpE9577z1YWlrKxwkhUKtWrWc88ZLdl9bIkSMxevRoeHh4QAiBjIwMhIaG4lghE8UOGTIE\n+/btg7+/P4QQGDhwoLyvR48eeOeddzBjxgx4eHjA7N8/Zh9//DGysrLQtWvXYucD2rRpE+zt7TF7\n9mx5m6OjI3766Se0bt0atWrVQlZWFhYuXIjjx48DKNnnJDU1FYmJifL8eYmJiTr9kmxtbaFQVEAP\nHYOmXy+JkydPEgBq3bq1oUN5+V25QtSggVRFY2ND9Pffei9CoyFq3Voq4uef9X55xgzK19e3yJqM\nUaNGyTVAbdq0ISEEffvttzrHfPvtt2RpaUkzZ86kPXv2kK2tLbVo0UKeQ6hu3bqFXruktTZffPFF\nkXPDCCEoKCiIvL29SQhB5ubmVKdOHZ397du3L3DN0NBQEkKQpaUlCSHozJkzRCSNuNKOzHqyOUcI\nQfXq1aPc3FwiIlqxYkWxMa1cubLYeyrJfe3bt0+uAQoJCdE5V1sT88MPP5Cvr2+B0VMZGRk0cOBA\nWrJkSYFyHz58SO7u7tS4cWPau3cvff755ySEIF9fXyIi+vrrrwst8+LFi6RSqejHH3+Ut3l4eBQ6\nf5Crq6vcTKZV3Ofk3Llz8n3Xr19f5zmYmprqDJcvT1wDpAfaXvlPZvSsnDRqJHWOHjIE2L9f6kS6\nfDnw/vt6K0IIqRZo0CBg4UJpdBj/atnLRAhR6PZ27drJCztrOyU/PWrpiy++gIWFBWbMmIHs7Gz0\n7t0bv/76K6pUqYLk5OQi/+Xu4uKC5s2bw8rKqtjYvvnmG9y+fRsbNmyASqXCuHHjoFar4e7ujlu3\nbsHV1RXXrl2DjY2NXJMjhICnpydsbGwKHWXVoUMHdOvWDQcPHoSDg4M8Ss3Z2RmWlpZ4++23oVQq\n5WPr1q0LNzc3+W/6+PHjcePGDXz33XdQKBQYOXIk6tati/bt2yM6OloebVaW+3JxcUFERARMTU3R\nvn17nXPbtm0r1xBRIaOkzMzMsHnz5kLLtbS0xF9//YVx48ahd+/eMDExwcKFCzF16lQAQK9eveDj\n44ONGzeiS5cuAICMjAx8+OGHaNOmDT755BMAQG5uLqKiovB+IX9r+/bti6CgIJ1txX1OatWqhSpV\nqqBr166oVq0aOnfujPr166N169ZwcHCQO6iXuwpJs15ye/fuJQDUvXt3Q4fy6sjNJfL2ftxZdMoU\naZue5OU9rmjasEFvl2XM4BYsWEAbN2585nGOjo46swo/LTw8nIKDg3UmrIuOjpY7ub5ozp49S9bW\n1jT7Be7ct2zZMqpevXqh+3x8fGjbtm00a9asUs2fU1Qn7/z8fPLy8iJTU1NasWIF+fn5kZeXF732\n2msUHR0tH5eUlETm5ua0c+fOAteOjIykgQMHFlpuYZ+TFwUnQHqwevVqAkAjRowwdCivHj8/ImNj\nKVPp2ZPo38nD9GHFCumyzZtX6KocjL0QGjVqRJ07d6YMPY9CY2Vz8uRJ6tmzp16veeXKFapVqxYJ\nIcjExIS8vb0pPj5er2W8iATRi7QyWeX07bff4r///S8+++wzzJ8/39DhvHqOHgXeeQe4f19qItu1\nC/h3csqyyMoC6tcHYmOBbduAt9/WQ6yMMcZeCDwRoh7Ex8cDAGrXrm3gSF5RnTsDp08DzZoBkZFA\n27ZS/6AyMjUFvvhCej9rFvDvRK6MMcZeApwA6cHdu3cBAHblMGsrK6F69YBjx4B+/YDUVODNN4Gl\nS6UeQmXwwQdAnTrAhQvAjh36CZUxxpjhcQKkB9rREiWZC4KVI0tLqa3qv/8F8vOBTz6RhnBlZZX6\nklwLxBhjLydOgPRAu3CedtIsZkAKBTBnDrB+vZS9rFoFtG8PXL9e6ktqa4EuXuRaIMYYe1lwAqQH\n2gSoRo0aBo6Eyd59V1pRvn594Px5oFUroJgZT4vDtUCMscouOztbZ6HYiqL9fnwRcQJURhqNRp7O\nmxOgF0yLFsCZM8DAgcDDh9LMhh9/DPy7yvXzeLIWaPv2coiVsQoUGxsLb29vhISEGDoUHZGRkRg/\nfjzCw8MNHYpMH8+qpPel0WgwbNgwKBQKjB07ttTlPSkrKwvx8fHo378/HB0di1364ll27twJlUoF\nJycnxMbGFnkcEeHevXtYs2YNbG1tMWbMmEIncDQ4w47Cr/wePHhAAMjc3NzQobCiaDRES5Y8ni+o\ndWuiq1ef+zI//yyd3qwZzwvEKrcPPvhAXszyWe7du0e9e/fWWTz1hx9+0DmmY8eO9OmnnxZ6/qNH\nj2jSpEl07ty5Z5bl6elJQghaunRpkcfMmzePrKysCl1OwsLCgmJiYujChQt04sSJAudOnTpVXgZk\n9erVZGtrW+h1tCvTEz3fsyrLfRERTZo0iYQQVLVq1QIrtN+8eZMUCgWtWLGiwHnahVa1k1ZmZGSQ\nr68vKRQK+Z46dOhApqamhS5l8aQJEybIi8U+6fDhw2RqairHNnLkyALHBAcHk729vc6z7Nq1K1la\nWtKjR4+IiOjEiRPUrFmzIpcEefqzVZiAgABSKBSFLrj6PDgBKqPk5GQCQNbW1oYOhT3LyZNEjo5S\nFmNhIU2i+Byzk2ZlEdWpI53+xx/lFyZj5c3R0ZEUCgWFhoY+89i+ffuSEIIGDRpE3333HdWtW5cU\nCoXObNK+vr7k6OhY6PnLli0jlUpFR48eLbac1NRUqlq1KikUCoqNjS30GO2aWBYWFrRo0SL69ddf\ndV7a+xk9ejS1aNGiwPlvvPEGjRkzhhISEuQ1wHx9fQtcJygoSD7neZ5Vae+LSEoMhBD0448/UmJi\nIjk5OZGrq6s8g/I///xDQgiysbGhK1eu6JyrTYBu375N2dnZcsJlb29Pc+fOpaCgILp16xadOXOG\nkpKS5PNyc3MpKytL51qTJk0iExMTunv3rrxNo9GQq6sreXh4UHp6On300UckhNBJMletWkVCCDIz\nM6Nhw4bRnj176NixYxQXF6eztpeLiwsJIei9994r8NwDAgJKNGO0tiwhBC1cuPCZxxeFE6Ayunfv\nHgGgGjVqGDoUVhLJyURDhz5eQqNfP6J790p8+i+/SKc1bqzXlTcYq1DampyYmJhijwsICCAhBPXr\n148ePHhARNLfvCZNmpCJiQn9888/REQUFBSks8jok0aNGkVNmjR5ZkzahUCNjIyKPCYnJ4esra1p\n5syZxV5r9OjRJISgxMREedv9+/fJ1NRU/sJs2rRpobUYTyvpsypKSe6LSFqY1cTEhFL+nc1+6dKl\nJISgvXv3EtHjBEgIQb169dI5d9SoUaRUKikhIYHWrl0rJxg5OTnFltm/f39yc3PT2XbhwgVSqVT0\n+eefy9tSUlJICEHLli0jIulZmpiY0ODBg4mIKC0tjaytrcnJyYnCw8OLLXPy5MmlWsrjSVlZWaRW\nq0kIUWTiXRLcB6iM8vPzAfBCqJVGtWrAhg3AunVAlSrAzp3SBIq7dpXo9PffB157DbhyBfj993KO\nlbFyoPm3F79arUadOnWKPI6IMGfOHLzzzjvYsWOHvIipjY0NQkJCUKNGDaxcuRIA0LBhQwCQO9n+\n9ddfyMzMRFZWFnbs2IFevXqVOK7CFjPVMjY2hqmpKa5evYr4+Hj5lZKSonOcdnHSw4cPy9sePnyI\n7OxstGnTBgBgYWGBqKgoxMXFydd5usNuSZ9VWe8LAPbs2QNLS0tYW1sDAPr06YNq1arhzz//BPB4\n0W0ACAwMRHBwsPxzXl4eOnToAFtbW9y4cQNqtRoBAQEwNjaGRqOR70/7faWN69KlS4iJidGJo3nz\n5hg+fDj8/PzkfXv27AEAODg4AJD6u/bq1QsHDx5EZmYmkpOTkZaWhpUrV8LV1RUAkJycjPj4eKSn\np+tc39LSEmlpabhy5YocV0JCwnP1ETIxMUH//v0BADk5OSU+72mcAJWR9pf2PL889gIYNkzq0fzG\nG0BCAvDWW9IK8wkJxZ6mUkkjwQDA17dU/akZw65d0kDFd98FBg8GevaUEuuaNYEuXYCZM4EvvwQ+\n+kj/ZZ8+fRqAlEwUZ8OGDYiLi8OqVasK7LOxscGHH36Ibdu2AQDMzc2hVCoRHh6OXbt2oXfv3pg9\nezZOnjyJBw8ewN7eHoD0d3LNmjW4ePEiYmJiEBMTI8+jpo1LpVIVGRMRIScnB1u3boWzszPUajXs\n7e3Rs2dPnRFOI0eOhIWFBRYsWCBvCwsLg0qlklcaz8nJwbFjx9CoUSOo1Wqo1Wp06tRJp3NvSZ9V\ncUpyX5mZmYiIiICnp6e8rWrVqujYsSMiIyMBALt374aVlRUuXLgANzc3eHt7y8nFtm3b5H+E16tX\nD+np6UhNTQUAbN26Vb4/FxcXjBgxAhqNBtevX8fNmzfRr1+/AvH07dsXycnJGDNmDADg7NmzEEKg\na9euAKSO1QMGDEBSUhISExNRvXp1WFlZITo6Wr6Gg4ODXK6npyfCwsIASM89NTUV7dq1k/c3adIE\n+59z9v6hQ4cCkEaZ3bhx47nOlZWpHopRRkYGASCVSvVCrnbLniEvj+j774nMzaW2rWrViFavLrZv\nUF4eUZMm0uFLllRgrOylsWDB41bY4l4KBVF2tn7LPnjwIAkhyN3dvdjjPD096auvvipy//fff0/G\nxsaUmppKRERt27alwYMHy32G3njjDfrll19IqVTSnTt3iIjI39+/yA7Hq1evJiEEDRgwoMgyjx49\nSkIICggIoKysLNq/fz/duHGj0GOHDRtGQgh5v6enJzk4OBCRtGq9EIJmz55d5Crpz/OsilOS+9I2\nkzVv3pxiYmJo/vz5pFQqSQhBTZs2JSJpNXht09HRo0fJ2NiYunXrRpcvXyYhBHXr1o2IiDIzM6lJ\nkyY0ceJEIpJWew8ODqagoCCaMmUKCSFo+/bt5OfnV6CjtVZKSgqpVCpq2bIlERF5eHiQEIL27dtH\nhw4dIicnJ/n3d+rUKSIimjVrFlWvXl3uUxQeHk5BQUH066+/klqtpn79+hERkbOzM3l6ehIR0enT\np+WO289L+8wUCgUdPny4VNfgBEgPzMzMCAA9fPjQ0KGw0oqKklaT137zdOtG9G//hsJs3y4dZmtL\nlJ5ecWGyl0NEBNH69dJrwwai3bulbXFxRJs3E/33v0S+vkRr1xJlZuq3bB8fHxJCyKOhCpOUlERC\nCLk/ytO0/YDGjBkjb+vatav8pahSqUgIQS4uLtS6dWv5mMzMTNq/fz8FBQXJr4iICCJ63JHX39+/\nyLhCQkLkzr7Pok1etH1Z3Nzc5NFN2i/PZ31xluRZPUtJ7ksbj/ZlZWVFW7dupU6dOpEQghISEnQS\nICKijRs3ykmSEILmzJkj7wsKCiJjY2P666+/dMrJzc0lc3NzWrJkiXxvhT0DbWdzbRKlTYC0r2nT\nptGsWbNICEHz588nIqlfjrOzc6Ejs8aOHUutWrUiIqJ69erRrFmzSv4An/HMOAEyMHt7ewJAt27d\nMnQorCw0GqLffyeqUUPKboyMiKZMISrkS0CjIWrTRjrM19cAsTJWSt7e3iSEoMaNG9O1a9fo7t27\n8is+Pp6ISO5IW5iDBw+Sq6srtW/fXucffYsWLSIhBDVs2JBiY2PJxMSEhBBFDo9/Wr9+/UgIQZ07\nd6bo6GiduBISEoiIaOXKlSSEoLNnz1J2drbOMU+ObiKSaj48PT2pRo0aFB0dTc7OznKH4n379pEQ\ngrZu3Up5eXk617n3xKCIkjwrfdyX9st81KhRFBISIrcmaDuhb9++vUACRET0/vvvy0nA0yPDhgwZ\nQs2bN6eMjAx524QJE0ihUNClS5eKTYCmTp1KDg4OdP/+fSKSEqB69erR7t275fvOzMwkKysr6t+/\nv3ze3r17SaFQ6Iz4+/vvv0mlUtG0adMoJSWFLCws6JNPPiEiKZF+8nnkP8f8IpwAvSDc3NwIAJ09\ne9bQoTB9SEwkGjeOSAgpw7GxIVq5Umr7esKhQ9JuU1OpAomxymDFihVFzsEihKCVK1fKQ7L79u1L\nv/zyC/n5+dHKlSupS5cu8pB47agwrZ9++olMTU0pMDCQiKSaD4VCQefPny9RXF988UWxcQUFBckJ\niTigfsgAACAASURBVLm5OdWpU0dnf/v27QtcMzQ0lIQQZGlpqTNK7fvvv5dHZj3ZnCOEoHr16lHu\nv0M8S/Ksynpf+/btk7/MQ0JCdM7V1sT88MMP5OvrWyABysjIoIEDB9KSQtriHz58SO7u7tS4cWPa\nu3cvff755ySEIN9//8X29ddfF1rmxYsXSaVS0Y8//ihv8/DwKHT+IFdXV7mZTOvbb78lS0tLmjlz\nJu3Zs4dsbW2pRYsWlJaWRufOnZPvu379+jrPwdTUVGe4/LM8WWvGCZABdevWjQDQvn37DB0K06ez\nZ4k6dXrcLNayJdG/7d1a2hH177xjoBgZK4XPPvtM7n8zZswY8vHxocDAQPLz85NrErZs2ULt2rWT\nv2S8vLxo3rx58oR2T8vIyKALFy7IP+fn59PJkydLHJNGo5H77ZiYmNBHH31E3377LQUFBdGKFSso\nOjqafvrpJ7K1taUxY8bQmDFj6P3336d169bRvn37KC0trdDraufEqVu3rrztzz//JCsrKxo5cqR8\nLT8/PwoKCipQs1OSZ1XW+woKCiIzM7MCc/IQEdWqVYt++OGHQmuAnuXevXtyDZSpqSktWrRI3nfm\nzBkSQtCHH34ob3v06BG1b9+eOnbsKG/LycmhunXr0m+//Vbg+jNmzCh0vqUff/yRTE1NSQhBffr0\nkWu64uLiqGrVqtS/f3/5uX/zzTc6TaElFRMTQwqFgiwtLen69evPda4WJ0B6MGjQIAJAGzZsMHQo\nTN80GqmThr3940Ro61Z5d3S0NKciQLRrlwHjZIwV6uzZs2RtbU2zZ882dChFWrZsGVWvXr3QfT4+\nPrRt2zaaNWtWqebPKaqTd35+Pnl5eZGpqSmtWLGC/Pz8yMvLi1577TWKjo6Wj0tKSiJzc3PauXNn\ngWtHRkbSwIEDCy03PDycgoODX+jBQYKIx2+X1cSJE7F8+XL89NNPmDx5sqHDYeXh0SNg4ULpJQQQ\nGQnY2QEAvv8emDpV+jE8XJpqiDHG9OnUqVOYOXMmAgMD9XbNyMhIeHh44N69e1CpVJgwYQK++OIL\n1KpVS29lvMh4HiA9qF69OgDIi6Kyl5CFBeDjA9y4AYweDXz+ubzrk0+Ajh2Bu3cBb2/DhcgYe3n9\nf3v3Hd9U1f8B/HOTpklHuinQXdpSVqHsJRt8GKKIDEUehoqCIgj8HicgCiooUxFkK7gZAipYNiKr\ngFIQoezSFuiA7qYjyff3x/He7lJoIG35vl+v87q3yU1ycinNJ+ece06bNm0sGn4AIDQ0FDdv3oTZ\nbEZOTg4WLlz40IQfgAOQRcirwHMAegh4egKffirC0MmTAAC1GlizBrCzA9atA9avt3IdGWOM3REH\nIAvgAPQQqlcPCA8H/p1aPiRE9I4BwAsvAJcvW7FujDHG7ogDkAVwAHqIqdXK7rhxwJNPAunpwNNP\nA5VYooYxxth9xgHIAjgAMUCMjV61CvD3B44dE2s5McYYq5o4AFkAByAmc3UFvv8esLEB5s0D/l1E\nWTAYgGKrVjPG2IOQm5tbZMHYByU5OfmBv2ZFcQCyAGdnZwBAenq6lWvCqoJ27YAPPhD7I0YAygLJ\ndnbAxYvAa68BV65YrX6MxcfHY+LEidi7d6+1q1JEdHQ0XnzxRZw5c8baVVFY4lxV9H2ZzWYMGzYM\nKpUKL7zwwj2/XmE5OTm4efMmBgwYAH9/f2zYsOGen2vLli2wtbVFYGAg4uPjyzyOiJCYmIgvv/wS\nnp6eGD16NKrkjDvWnYaoZsjMzCQAZGdnZ+2qsCrCZCLq00dMkNiuXbEVvRcuFMt8Dx1K9O/U/Iw9\nSM8//7yyqOWdJCYmUr9+/ZQZoVUqFS1cuLDIMR07dixzza+srCx6+eWX6a+//rrja8mzNn/22Wdl\nHjN79mzS6/WlLivh4OBAcXFxFBUVRUeOHCnx2MmTJysLm65evZo8PT1LfR55hXqiuztXlXlfREQv\nv/wySZJEzs7OJVZqv3z5MqlUKlq2bFmJx8kLrsorq2dnZ9OMGTNIpVIp76lDhw6k0+lKXdKisLFj\nxyqLxha2f/9+0ul0St1GjBhR4pjdu3eTj49PkXPZrVs3cnR0VGYQP3LkCIWFhZW5NEjx363ybN++\nndzc3Ein05W6qv2dcACyAJPJRJIkEQAyFlsvij28kpOJfH1FCJo0qdidQ4cWzCzdrRvRpk1E/64/\nxNj95u/vTyqVig4ePHjHY/v376+s//XJJ5+Qn58fqVQq+v7775VjZsyYQf7+/qU+fsmSJWRra1tk\ngczSpKamkrOzM6lUKoqPjy/1GHltLAcHB5o3bx6tWrWqSJHfz6hRo0pdoqFz5840evRoSkhIUNYC\nmzFjRonniYiIUB5zN+fqXt8XESnrry1atIiSkpIoMDCQGjVqpMykfPbsWZIkiWrVqlViCQ45AMXE\nxFBubq4SuHx8fOijjz6iiIgIunr1Kh0/frzIorH5+fkllt94+eWXSavV0o0bN5TbzGYzNWrUiLp0\n6UKZmZk0fvx4kiSpSMhcuXIlSZJEdnZ2NGzYMPr111/p0KFDdP369SJrfDVs2JAkSaLhw4eXOO/r\n1q2r0MzROTk5yvIi8nb06NF3fFxxHIAsRKVSEQBlET3GiIgOHxaLyhdbQUM0EX30EZGjY0EQ8vMT\ntyUlWa2+7OEgt+TExcWVe5y8GvkTTzyhLH6amJhIjRs3Jq1WS2fPniUiooiIiCKLjRY2cuRIaty4\n8R3rJC9uaWNjU+YxeXl55OLiQtOnTy/3uUaNGkWSJFFSof9LycnJpNPpaO7cuURE1KRJk1JbMYqr\n6LkqS0XeF5FYoFWr1VJKSgoREX322WckSZKygr0cgCRJor59+xZ57MiRI0mtVlNCQgKtXbtWCRh5\neXnlvuaAAQOoWbNmRW6LiooiW1tbevPNN5XbUlJSSJIkWrJkCRGJc6nVamnIkCFERJSWlkYuLi4U\nGBhIZ86cKfc1X3nllXta0qOwmTNnkkajUdbfbN269T0FIB4DZAEmkwlmsxmSJEFd6LJoxtq1Az75\nROyPHi2GAAEAVCoxm3RsLLBwIRAcDFy7Brz1FuDjI2abPn7cWtVmNZjZbAYAeHl5wdvbu8zjiAiz\nZs3CwIEDsXnzZuj1egBArVq1sHfvXri7u2P58uUAgPr16wOAMsh2+/btMBgMyMnJwebNm9G3b98K\n16tt27ZlHqPRaKDT6XD+/HncvHlTKSnFLi7o1KkTAGD//v3KbRkZGcjNzUXr1q0BAA4ODrhy5Qqu\nX7+uPE/xAbsVPVeVfV8A8Ouvv8LR0REuLi4AgMceewyurq7YunUrAMBoNCrH/vbbb9i9e7fys9Fo\nRIcOHeDp6YlLly7By8sL69atg0ajgdlsVt6f6d95y+R6nT59GnFxcUXq0bRpUzz77LNYsWKFct+v\n/17N4evrC0Bc+NO3b1/s2bMHBoMBt2/fRlpaGpYvX45GjRoBAG7fvo2bN28iMzOzyPM7OjoiLS0N\n586dU+qVkJBwV2OERowYgZUrV6JXr15ISUlBVFQUmjZtWuHHKyoVwxgRiT5uAKTVaq1dFVYFmc1i\ntXiAKDycyGAo5SCTiWj7dqK+fYkkqaBVKDyc6LPPiG7ffuD1ZvfP7su76bXfXqPXfnuNJmyfQGN/\nGUujN4+mYRuHUb9v+tGQ9UNo5YmVlJVX+srrlSF3tQQGBpZ73DfffEN6vV5pkSjuvffeU7q9EhIS\nyMbGhpYsWUJbt24lSZLorbfeon379indOkSiK2XNmjUUFRVFsbGxFBsbq6wU/v333ytjRspiNpvJ\nzc2NNBoN2dvbK+N1WrduTZmZmcpx+fn55OjoSG3atFFu27hxI2m1WkpMTCQioubNm5NarSZHR0el\nZSU0NLRIS09Fz1V5KvK+srOzycvLS2lRISK6ffs29e/fX3ncRx99RE5OTnTq1Clq3rw5NW7cmDIy\nMoiIyM7OTjluzZo15OzsrPy7/fjjj8r7CwkJoeHDh5PJZKLo6GiSJImee+65EvXZtGkTSZJEPXv2\nJCKiSZMmkUqlUs6xwWCgr7/+Wul2S0tLIycnJ1q9erXyHPK/j16vp+7du9OJEyeU55IkiZycnJR6\nubu7F+l2vBuDBg2iZs2aKWOM7gYHIAuIiYkhAFSnTh1rV4VVUampREFBItO8+OIdDr5wQQwacnMr\nCEJaLdGwYUR79oiwxKq1j//4mDADdyxuc9zoSsoVi772nj17SJIkatmyZbnH9ejRg6ZNm1bm/QsW\nLCCNRkOpqalERNSmTRsaMmSIMmaoc+fOtHTpUlKr1XTt2jUiEh/OZQ04Xr16NUmSRE899VSZr3ng\nwAGSJInWrVtHOTk5tHPnTrp06VKpx8pjQ+T7e/ToQb6+vkREFBsbS5Ik0cyZM8tcLf1uzlV5KvK+\n5G6ypk2bUlxcHM2ZM4fUajVJkkRNmjQhIrEqvNx1dODAAdJoNNS9e3f6+++/SZIk6t69OxGJcNK4\ncWMaN24cEYkxqrt376aIiAglfPz000+0YsWKEgOtZSkpKWRra0vNmzcnIqIuXbqQJEm0Y8cO2rdv\nHwUGBir/fpGRkUQkArGbm5sypujMmTMUERFBq1atIi8vL3riiSeIiCgoKIh69OhBRETHjh1TBm7f\nC/m8ffvtt/f0eA5AFrBv3z4CQB06dLB2VVgV9uefIscAROvWVeABBgPRd98R9exZEIQAonr1iD74\ngOgexyQw6zsWf4zmH5pP8w/NpwWHF9DnkZ/TyhMrae3JtbTl3BZafnw5tV7empoubVqhQaF34913\n3yVJkpSroUpz69YtkiSpzNYfeRxQ4XEX3bp1Uz4UbW1tSZIkatiwIbVq1Uo5xmAw0M6dOykiIkIp\n//zzDxEVDORds2ZNmfXau3ev0upwJ3J4kceyNGvWTLm6Sf7g3L9/f7nPUZFzdScVeV9yfeSi1+tp\n48aN9Mgjj5AkSZSQkFAkABGJliU5JEmSRLNmzVLui4iIII1GQ9u3by/yOvn5+WRvb0+ffvqp8t5K\nOwfyYHM5RMkBSC5Tpkyh9957jyRJojlz5hCRGJgcFBREzzzzTInne+GFF6hFixZERBQQEEDvvfde\nxU9gGUwmEz322GMUEBBQZGD33eAAZAGrV68mAPTss89auyqsilu2TGQYe3uiO4wVLOrKFaLp04l8\nfAqCkEpF9OijRGvXEv3bFM5qDrPZTLey7+0Pe3kmTpxIkiRRgwYN6MKFC3Tjxg2l3Lx5k4hIGUhb\nmj179lCjRo2offv2ShcMEdG8efNIkiSqX78+xcfHk1arJUmSyrw8vrgnnniCJEmiTp06UWxsbJF6\nyd1ky5cvJ0mS6MSJE5Sbm1vkmOIfgiaTiXr06EHu7u4UGxtLQUFByoDiHTt2kCRJtHHjRjIajUWe\nR+4iq+i5ssT7kgPQyJEjae/evUrolQeh//TTTyUCEBHRc889pwzSLn5l2NChQ6lp06aUnZ2t3DZ2\n7FhSqVR0+vTpcgPQ5MmTydfXl5KTk4lIBKCAgAD65ZdflPdtMBhIr9fTgAEDlMdt27aNVCpVkSv+\n/vjjD7K1taUpU6ZQSkoKOTg40IQJE4hIBOnC58N0F63bx48fV1ql7hUHIAuYOnUqASi3uZgxIjEe\naPhwkV8aNryH3GI0irFCgwYRaTQFYcjenuiZZ4h+/ZXoDld+sIfbsmXLypyDRZIkWr58uTL2pX//\n/rR06VJasWIFLV++nLp27apcEi9fFSZbvHgx6XQ6+u2334hItHyoVCo6efJkher11ltvlVuviIgI\nJZDY29uTt7d3kfvbt29f4jkPHjxIkiQp43zkq9QWLFigXJlVuDtHkiQKCAhQruatyLmq7PvasWOH\nEoD27t1b5LFyS8zChQtpxowZJQJQdnY2DRo0iD799NMSr5uRkUEtW7akBg0a0LZt2+jNN98kSZJo\nxowZRET0/vvvl/qap06dIltbW2XcFpEIQKXNH9SoUSOlm0z24YcfkqOjI02fPp1+/fVX8vT0pPDw\ncEpLS6O//vpLed/BwcFFzoNOpytyuXx50tPTlcvpp0yZQitWrKCVK1cqXa0VxQHIAp599lkCUG4T\nJ2OyjAwRfgCiZ58Voeie3LpF9MUXRI88UrSLrFYtoldfJTp6tBJPzmqy119/XRl/M3r0aHr33Xfp\nt99+oxUrVigtCRs2bKC2bdsqH1C9evWi2bNnlznYNDs7m6KiopSfTSYTHT16tMJ1MpvNyrgdrVZL\n48ePpw8//JAiIiJo2bJlFBsbS4sXLyZPT08aPXo0jR49mp577jn65ptvaMeOHZSWllbq88pz4vj5\n+Sm3bd26lfR6PY0YMUJ5rhUrVlBERESJlp2KnKvKvq+IiAiys7MrMScPEVHt2rVp4cKFpbYA3Uli\nYqLSAqXT6WjevHnKfXILyksvvaTclpWVRe3bt6eOHTsqt+Xl5ZGfnx999dVXJZ7/jTfeKHW+pUWL\nFpFOpyNJkuixxx5TWrquX79Ozs7ONGDAAOW8f/DBB0W6Qivi/PnzpNVqi0z0KEkSTSox4Vr5OABZ\nQPv27QkA7du3z9pVYdXEmTOi0QYQ3WKVdvky0axZRA0aFA1DQUFEb7xBFBnJYYg9lE6cOEEuLi40\nc+ZMa1elTEuWLCE3N7dS73v33Xdp06ZN9N57793T/DllDfI2mUzUq1cv0ul0tGzZMlqxYgX16tWL\n6tWrR7Gxscpxt27dInt7e9qyZUuJ546OjqZBgwaV+rpnzpyh3bt3W3wMmyVJRFVxgY7qpU6dOkhI\nSMC1a9eUeRIYu5Ovvwb++19AqwUOHAD+nZ6kcoiAP/8EvvkG+PZbICGh4D5/f+Cpp0Rp107MRcQY\nqxYiIyMxffp0/PbbbxZ7zujoaHTp0gWJiYmwtbXF2LFj8dZbb6F27doWe42qjANQJWVnZ8PBwQEa\njQYGg4EnQmR35aWXgOXLgdq1gSNHgIAACz65yQT88QewYQOwaRNw/XrBfV5ewIABQL9+QNeugL29\nBV+YMcaqPg5AlXTu3Dk0bNgQQUFBuKhM88tYxeTlAX36AHv2AA0bAgcPAq6u9+GFzGaRsDZuFIHo\n2rWC+3Q6oFs3oG9fUerVuw8VYIyxqoUDUCXt2bMHPXr0QJcuXbBv3z5rV4dVQ6mpwCOPAGfOiMaY\n334T3WL3DZFYZuOXX4Bt20ouuREaCvTuDfToAXTuDDg738fKMMaYdfAggEq6/m+3gpeXl5Vrwqor\nFxeRQ+rWBfbtA154QWSU+0aSxICj994Djh0Dbt4EvvwSGDJEhJ3oaGDRIuDxxwE3N6BtW7FG2c6d\nQHb2fawYY4w9OByAKikxMREA4OnpaeWasOrMzw/49VfAwUEMjv7f/+5zCCqsdm1g5Ejghx+ApCRg\n/35g2jSgY0cxUDoyEpg9G3j0UZHWunQR4enAAdGHxxhj1RAHoEqSV+i1tbW1ck1Ydde8uRieY2MD\nzJsHvP++FSqh0Yhur/ffFwOoU1JEn9zrrwOtWgFGI/D778CMGeI4FxegVy9g1ixxe06OFSrNGGN3\nz8baFajuzGYzAEDFlxQzC+jdW1y9/vTTImPo9cDkyVaskKMj8J//iAKIQPT772LU9p49wN9/A7t2\niQKIwUvt2olWos6dgfbt+QozxliVxAHIQuQgxFhlDR4shtqMGgVMmSLyw9ix1q7Vv1xdgSeeEAUQ\n8wwdOCC6zfbvB06fLtgHRItS69YiDHXpIrrV9Hrr1Z8xxv7FAaiS5LE/8lggxixh5EggMxMYPx4Y\nNw7IzQUmTrR2rUpRuzYwaJAoAHDrlug6k0PQyZPAoUOizJ4txhS1aCHCUJcu4vK3+3LdP2OMlY8v\ng6+k7du3o2/fvujVqxd27Nhh7eqwGuazz4AJE8T+Bx8Ab79t3frctbQ0MbmRHIiOHxcTNBYWGiqu\nNGvTRmybNgV4TB1j7D7jAFRJUVFRCA8PR5MmTXD69GlrV4fVQKtWAWPGiKvCJk8GPv4YqLYTjmdm\nAocPizD0++/A0aMlryTTasWI8DZtRPdZs2ZAgwaiO40xxiyEA1AlJSUlwdPTE25ubrh165a1q8Nq\nqO++A0aMEBdhDRggLpV3cLB2rSwgLw84dUoEochIsY2OLnmcRgM0aiTCULNmopWoWTOgVq0HX2fG\nWI3AAaiSzGYzdDod8vPzYTAYoNPprF0lVkPt3QsMHChmjm7ZEti6VSzpVeOkpoqusqNHxcKuUVHA\npUulH+vpKVqHQkPFVi7+/tW4mYwx9iBwALIAPz8/xMbG4vLlywgMDLR2dVgNdu6cWL/08mXx2f/d\nd0D37tau1QOQmSmuMIuKEuXUKVEyM0s/XqsFQkIKAlFQkFhpNiAA8PERky0xxh5qHIAsoF27djh6\n9Cj++OMPdOzY0drVYTVcUhIwdKhoEVKpxJyFb70l9h8qZjMQHy9SoVyio8U2Pr7sx6nVIgTJgcjf\nXzSl1alTUGrXFovEMsZqLA5AFvD444/j559/xqZNm/Dkk09auzrsIWAyAe++K64MA8Q8hWvWiPXE\nGICMDBGG5EB05Qpw9aoo169XbJ0RZ+eCMFSrFuDhUbAtvC9v7ezu97tijFkQtwNbQK1/B2ImJydb\nuSbsYaFWi9UnOnYEhg8HIiKAJk2ApUvFmqYPPb1eLN3RqlXJ+3JzgdjYgkB09apYEDYhQWzlkpYm\nSmmDsktjb19+QPL0FIu++fkB7u5iUVrGmNVwALIADw8PAOKKMMYepD59xFCY558XIWjoUGDLFjF/\nkJubtWtXRWm1QHCwKGUxm8WyH3IoSk4WJSmp6FbeT0oS03fHxIhyJ3Z2BWHI379gX+6W8/bmcUqM\n3Wf8P8wCeB0wZk3e3sD27cAXXwD/939iLbFdu4CFC8WaYjWxoSE6ORqhHqH37wVUKtFK4+4uLr+/\nEyLR7VZeULpxQ7Q8XbtW0LJUVuuSWg34+opwJIeiwoUHcjNWafw/yAIMBgMAwI7HADArkSSxZEbP\nnsBzz4nVKIYNA776SnSL1aSLE386+xOe+vEpTG4/GTO7zYSdpgr8v5MkwMlJlHr17nx8WpoIQoWL\n3Hokj1OSu+fkddUKKz6Qu3Aw8vQUxcODQxJj5eD/HRbAAYhVFSEh4vNy9Wrgf/8T3WKNGwPvvCNm\nka4Jv6IXb1+EJEmYd3gefj7/M1Y/vhod/arZ1ZfOzkBYmCilyckpOU7p6tWiAUkOTKUFJJmbmwhD\ntWoVhCJn54KwVrwUvs/OrmY2HzL2L74KzAJGjBiBdevW4csvv8TIkSOtXR3GAIjhK5MmibmCANGb\nMncu8NRT1f9z7Vj8MYzaMgr/JP0DCRJebv0yPuzxIZy0Ttau2oNR2kDumBggLk50uyUmioVpzeZ7\nfw21umJBqbzi7CymLK/uv3CsRuIAZAGDBw/Ghg0b8MMPP2AIX4LDqpg9e0QQOnVK/Ny5M7BggViU\nvTrLNeZi1oFZmP3HbBjNRnjpvbCo9yI81fApSPyBK+ZKuH1bhCE5FCUnA+npBSUtrejPhW/PzbVM\nPQp3D95tkNLpRNFqRZH3dbqKde8RifNgNAL5+ZXblnWbjY2ok61tQT0LF51OtKYVLjrdQzhxV9XD\nAcgC+vXrh23btuHnn3/GY489Zu3qMFaCyQSsXAlMnSo+AwFg8GAxiWKDBtatW2WdSjiFF39+EUfj\njwIA/hP0H3za51PUd69v5ZpVc3l5YmB3WSGpvJKWJh6bng5kZd2f+qlURUNE8Y8yosq1gN1vWq2Y\nOqF4OCpcit9vYyNK48bikk9WKRyALKB79+7Yu3cvdu3ahR49eli7OoyVKTVVTJ64eLEYZqJSASNH\nikkV/f2tXbt7ZzKbsOLPFXhr91tIzUmFRqXBpPaTMLXTVOi1emtX7+FmNIolS8prcSotQKWni1ao\nnJyS25yciocblUqEBo2mYtu7PTY/X4TF3NyCUvjnnBzAYBAlO1tsK9u6NnAgsHFj5Z6DcQCyhPbt\n2+PIkSM4ePAgOnToYO3qMHZH8fHAzJnAqlXi80mjERMqvvGGWFe0ukrMSsSbu97EmpNrAAC1HWpj\nVvdZGB0+GmoVL45aoxiNJVt9ind9Fm8lqirM5qLBqHCRQ1JpRe52a9CAW4AsgAOQBYSHhyMqKgp/\n/vknmjdvbu3qMFZhFy8CM2aIgdJms/j8GDhQrC3WsqW1a3fvIuMjMfG3iTgSdwQAEOYZhjk956B3\ncG8eH8QYAwBUwWhc/WRnZwMA7O3trVwTxu5OcDDw9ddiPr4XXxQtQRs3ihUkunYV+0ajtWt599p4\nt8Gh5w7h24Hfwt/ZH6cTT6Pvt33R7atuOBR7yNrVY4xVAdwCZAEBAQGIiYnB5cuXEViTZpxjD53r\n18UVYl98IYZtAGJuvbFjgTFjxFQy1U2OMQdLji3BBwc+wG3DbQBA35C+eK/re2jlVcpaYYyxhwIH\nIAvw9fVFXFwcYmJi4OfnZ+3qMFZpaWnA2rVisPT58+I2jQZ44glg1Cix+nx1m2Q4NScVcw/NxcIj\nC5GVL65Meqz+Y5jWeRraeLexcu0YYw8aByAL8Pb2xvXr1xEXFwdvb29rV4cxizGbgd27RRD6+eeC\nMad16gD//a8IQxVZKqsqScpKwseHPsaSY0uQnS+6r3vV64W3HnkLXQO68hghxh4SHIAsgFuA2MMg\nLg5Ytw748suCViEAaNYMGDJElPIWWK9qErMSseDIAiyOXIzMPNHf19qrNV7v+DqebPAkXzXGWA3H\nAcgCGjZsiHPnzuHMmTNoVN2+DjN2l4iAI0eANWuAH38U3WWyFi1EEBo4UKxLVh2kGFLwWeRn+PTo\np7hluAUACHQJxMS2EzG6+eiHZ3kNxh4yHIAsoE2bNjh27BiOHDmCtm3bWrs6jD0wubnAjh0iCG3Z\nIib/lTVoAPTvL0r79lV/zFB2fjbW/LUGC44swKWUSwAAva0eo8JH4ZXWryDUoxpPkMQYK4EDRT8i\nYgAAIABJREFUkAX07NkTu3fvxo4dO9CrVy9rV4cxq8jJAX77DVi/Hti2Tcw6LXNzA/r0AXr1Anr0\nEFeWVVUmswk/n/8ZC48sxP6YgpXWe9bribEtx+Lx0MehUWusWEPGmCVwALKAJ598Eps3b8bGjRsx\ncOBAa1eHMavLzwcOHgS2bhWDpy9eLHp/aCjQs6cojzwCeHhYp553cirhFD6L/Azfnv5WGTBdx7EO\nRoWPwnPhzyHEvZr08zHGSuAAZAEjRozAunXr8OWXX2LkyJHWrg5jVQoRcO4cEBEhrijbt69gjiFZ\nSIjoJmvXTmybNKlaXWYphhSsO7UOXxz/AmeTzyq3P+L3CEY2G4nBjQbDWedsxRoyxu4WByALGDdu\nHL744gssWbIE48aNs3Z1GKvS8vOByEgRhnbvBo4dE8scFWZvL0JQ06aiNGsGhIUBrq7WqbOMiHAo\n9hBW/LkC6/9Zr7QKadVa9A/tj2eaPIM+wX1gp7GzbkUZY3fEAcgCpkyZgvnz52Pu3LmYMmWKtavD\nWLWSnw+cOiWuLDt8WGwvXSr9WC8v0X1Wv74o8n5AgJio8UHKyM3AxrMb8VXUV9h/dT8I4k+po60j\n+oX0w1MNn0Lv4N68Gj1jVRQHIAuYNm0aZs2ahffeew/Tp0+3dnUYq/Zu3wZOnxbB6NQpICoK+Pvv\nki1FMhsboF69glAUHCxKUBDg63v/u9Pi0uPw/d/f44czP+D49ePK7bZqW3QL6IbH6j+GPsF9EOQW\ndH8rwhirMA5AFvDhhx/inXfewZtvvomPPvrI2tVhrEYymYBr18TCrefPF5ToaHF7WWxsRAtRUFBB\nKJJLvXqAnYV7q66kXMGms5vw07mfcCj2kNIyBABBrkHoFdQLPQN7omtAV7jbu1v2xRljFcYByAJm\nzZqFadOm4Z133sGsWbOsXR3GHjoGg7jSTA5EFy+KbrRLl4D4+PIf6+1dNBTJwahuXbHkR2W61hKz\nErH9wnb8euFX7Lq8Cyk5KUXuD/MMQ5eALujk1wmd/Dqhrr7uvb8YY+yuVKHrLKqv/Px8AIDmQQ9C\nYIwBEK04YWGiFGcwAJcvFwQiuVy8CFy9KgJSfDzw+++lP3etWiIMyYHI1RVwcSkorq6ihcnHR/ys\nUhU81tPBEyPDR2Jk+EgYzUacuH4COy/vxJ4re3Ao9hBOJ57G6cTTWBy5GICYgbqDbwelNPFsAhsV\n/5lm7H7g/1kWkJeXB4ADEGNVkZ0d0LixKMUZjUBsbMlgFBMD3LgBJCQASUminDp159fSaIDu3YEB\nA4AnnhChSWajskFbn7Zo69MWUztPRY4xB8fij2F/zH4cuHYAh2MP40rqFVxJvYJvTn8DAHDQOKCN\ndxu082mHdj7t0Na7LWo71rbQmWHs4cZdYBYwduxYLFu2DJ9//jlefvlla1eHMWYhJhOQmCjC0I0b\nwM2bYu2z1NSyS0YGkJ4u5j9q0QLo1w/o3Rto0wZQl7O+qslswunE0zgUewiHYg/hcNxhXE65XOK4\nAJcAtPNphzZebdDWpy2a12nOl91XApH4dzYaC7ZmswizcincqsdqDg5AFjB48GBs2LABP/zwA4YM\nGWLt6jDGrIxIdL2lpxcEIoMB0GoBJyfRrebmdufnSchMwNH4ozgcdxhH444iMj4SWflZRY6xUdmg\nae2maO3VWhTv1mhUq9FddZ1lZAC3bol6ZmcDeXniPQCAJAG2toBOV1C02oKtjY0ICGq1OPZ+IhIB\nJSdHlPR0EUjlrVwK/5yeLkpmpnifhbcGgwg9d6JWiyAknwd7e9GyKJfiP1fmNvlnDl33HwcgC+ja\ntSv279+PXbt2oUePHtauDmOshjKZTTiTdAZH4o4gMj4SR+OP4kzimSJXmgGAzkaHZrWboUXdFmhe\npznC64SjsWdj2Gvs7+l1b94UXYVxcQVjpgrvx8cDWf/mMrVaFDkUyduybpMkEWwKF6Bg32wWi+7K\nocdsrswZLJ1aLYKcjU1BnYxGEQT/HeL5wGk0BecGKHpuevcGfvnFOvWqSTgAWUBgYCCuXr2K6Oho\n1K9f39rVYYw9RDJyM/DnjT9x7PoxRMZH4sSNE6V2nakkFYLdghHmGYbGno3R0KMhGno0RLBbMBxs\nHSpdj9RUEYri4grCUmysmKLg2jWxn5NT6ZeBSiVaSOTWNGdnUQrvF//ZyQlwdAT0elHkfTu7O7dc\nya1O+fkFQcxgKCjZ2UV/ttRt5endG9i+vfLn8mHHAaiS8vPzodPpQETIycmBra2ttavEGHvIpRhS\n8NfNv3Di+gmcTDiJkzdPIjo5GiYqvb/HW++NYLdgBLkFIdAlEIEugfB38Ye/sz/q6uta7Eo0s1mE\nCbnk54tiMoliNoswIgcSSRKBp3DXW1VaI+5+IRKtT7LC50Mu5Y0nYxXDAaiSLl26hODgYPj5+SEm\nJsba1WGMsVLlGnNxNvksziSewZmkMzibfBZnk87icspl5JvL7udRSSp46b3g4+SjFG+9tyhOBVud\nje4BvhvGKu8hyNL31+XLoqm5Xr16Vq4JY4yVTWujRXidcITXCS9yu9FsRExqDC6lXMKl25dwNe0q\nrqRcQUxaDK6lXcPNzJuIS49DXHpcuc/vbucOHycf+Dr7iq2TL/yc/ZStt5M3bNXcQs6qDg5AlXT+\n/HkAQHBwsJVrwhhjd89GZYMgtyCxTlkpS5XlmfIQnx6P+Ix4xKbFIj4jHnHpcYjPiFduv55xHbcM\nt3DLcAtRCVGlvo4ECV56L/g5+xUpvk6+8HX2ha+TLzzsPSDd70vJGPsXB6BKOnfuHACgQYMGVq4J\nY4xZnq3aFoGugQh0DSzzGDOZkZiViNi0WMSlxyE2PVaUNLGNSY3BjcwbIjRlxONw3OEyX8tb7w0v\nvRe89F6o41gHtR1qw9PBE7UcasHdzh3u9u5w1bnC1c6Vu91YpXAAqqTo6GgAQGhoqJVrwhhj1qGS\nVKjjWAd1HOugtXfrUo/JN+UjPiMe19KuISY1BrHpsbiWdk0JSnHpcUjJSVFmw64IW7Ut9LZ66LV6\n6G31sNfYw8HWAXY2drDT2ClbnY2uoKgL9uVj7DX2sNOIrVwcNA5wtHWEg60DL0dSQ/Eg6Ery8/ND\nbGwszp8/j5CQEGtXhzHGqq2svCxcz7iulISsBNzMvImk7CQkZSWJbrbsW0jJSUGKIaXcwduWZKu2\nFWFI4wAHWwc4aByUoFQ8aGnVWmhttLBV20Kj0kCj0sBGZQMblQ3UKjXUkhoqSaUUSZIgQSqxBVBm\nd6Cfsx86+3d+IO+9JuMAVAlZWVlwdHSERqNBdnY2bB6G6zMZY6yKyDHmIC0nDVn5WcjIzUB2fjay\n8rNgyDcgOz8bOcYcGIwG5Bpzi2zl23OMOTDkG2AwGmDINyArP0s8R14WsvKzlK2Z7sPsi5XwVMOn\nsGHIBmtXo9rjT+xKiIsTV0X4+flx+GGMsQdMZ6ODzvH+jgMiIuSacpGVl4XMvEwlGBmMBmWbnZ9d\nJGTlmnKRZ8pDvjkfRrMR+SaxNZEJJrMJBFK2ZjKDiEAgZSu/bllae5XezcjuDn9qV0JiYiIAoHZt\nXp2ZMcZqIkmSlDFD7vbu1q4OsyBebq0S5ADk6elp5Zowxhhj7G5wAKqEzMxMAIBer7dyTRhjjDF2\nNzgAVYI87sdoNFq5Jowxxhi7GxyAKkGj0QAQC6IyxhhjrPrgAFQJ9vb2AIC0tDQr14Qxxhhjd4MD\nUCXIC6DKC6IyxhhjrHrgiRArIScnB/b29lCpVDAYDEqXGGOMMcaqNm4BqgSdTgcfHx+YTCbExMRY\nuzqMMcYYqyAOQJUUHBwMALh48aKVa8IYY4yxiuIAVEkcgBhjjLHqhwNQJXEAYowxxqofDkCVxAGI\nMcYYq344AFWSHIAuXbpk5ZowxhhjrKL4MvhKSk1NhaurK+zt7ZGZmQlJkqxdJcYYY4zdAbcAVZKz\nszMcHR2RnZ2NlJQUa1eHMcYYu+9iYmKq/TqYHIAqSZIk+Pj4AADi4+OtXBvGGGPs/kpPT0e3bt3Q\nuXPnar0SAgcgC9DpdAB4UVTGGGM1X0xMDPLy8nD48GGEh4fjq6++QnUcTcMByALkcT/V8ReAMcYY\nuxthYWE4deoUBg8ejIyMDIwaNQpDhw6tdsNAbKxdgTu5evVqkZaVwoOMiw84vpf7LHGc3A+alJSE\nGzdu3NfXuh/PIQc3IiqyDwAmkwk5OTnIzc1Fbm6usl/abfn5+cjLy0N+fn6J/fz8fKhUKtjY2ECj\n0RTZarVa6HQ66HQ62NnZldgv7TYbmzv/6hIR8vLyitT5TvtGoxEqlQqSJEGlUiml8M82Njb3VCRJ\nKnGOK7Nviee4077ZbIbJZCpRjEZjqbdb4n75PrkOkiSVWuR/l/KK/F7MZnO5+/d6PxFBo9HA1tZW\nKVqttsjPZd0n75e1lfflNQbL+vcvfl/hc1b897es22oCIoLRaITBYEBOTk6pW5PJpLzfsn6n1Go1\n1Go1bGxslP3SfpaL/DtrNBqRn59f5r78O202m4sU+TZA9CZotVrlb2LxfXt7e+h0OqhU1m+7cHNz\nww8//IB+/fph/PjxWL9+PQ4fPoy1a9eiW7du1q5ehVT5q8CCg4P5EnNWglqtVsKQVqtV/sjIf3Tk\nLWOsfPKH+/0o8pecwgVAkQBZVjGZTEW+oJRWDAZDkYAjB4maTv4yaG9vD3t7e2W/+Lb4bVqttsx/\nl3v5t5RLXFwcJkyYgBMnTkCSJLz66quYNm0aPDw8rH2qylXlA1CPHj1w7do1AEW7mIpX+17us9Rx\nycnJMBqNcHNzK7Ii/P2sk6XeV+FvQ/K28L5KpVJCRnnfTgp/y9VoNMq3Ynlfo9Eo39CKfzPKy8sr\n81tbWd/kKvqHztbWVmk9KlzvwtvC+2q1usQ3/8Lf1Aq3UtxNKd6KWdb5vtt9Sz9fad+Oi3/jLeub\nsCXuL3yfSqW644dk8daY4qV4S1Fp+5W5HxAtwHl5eeWW3NzcEvsV3ebl5d3x37/4fYXPT/FWrOK/\n2zVJ4S9GhbeFW47v9PtUXktlaS2X8u9t8Zbt4vvy73RZLctEVG4re+G/gdWBo6MjMjIyrF2NclX5\nAFQdtGnTBseOHcORI0fQtm1ba1fnoZCfn6/8QcjNzVX+2BQucrcTY6x08od+ad02liyFW2aBsrs1\niwdv+QtKaUWr1ZYIORXpGq/uzGYzcnJykJ2dDYPBUO62+H5eXt49/btV5LisrCxkZWUpobpWrVpI\nTEy08tkqX83/bXkAeBD0gyeHHL1eb+2qMFZtFW7h02q11q4OqwCVSqV0b1UF+/fvx//93//h+PHj\nAIAmTZpg7ty5+M9//mPlmt2Z9UdS1QAcgBhjjD1Mzp49i8cffxxdu3bF8ePH4eXlhVWrVuHkyZPV\nIvwA3AJkERyAGGOMPSxiYmIQFhYGk8kEBwcHvPHGG5g8eTIcHBysXbW7wgHIAjgAMcYYq6mKT7Hg\n7++PoUOHQq/XY8aMGahTp44Va3fvOABZAAcgxhhj1V1ycjKSk5MRFBSkXLlb+OpC+ZiTJ08iJCQE\nXbp0QZ06dYpcUVydcABijDHGHlJyeDlw4AB69+6Njh07YsOGDdBoNMjKysLZs2dx5MgRnDp1CocO\nHUJ0dDRUKhVCQ0PRqVMnACUn160uOABZALcAMcYYq8quXbsGPz8/fPnll/juu+8we/ZsNG/eHEaj\nERqNBp6enmjQoIEy59CKFSvwv//9D/b29ggKCsLp06fh4eGBn376CV5eXggKCoKLi4u131alcACy\nAA5AjDHGqorY2FgcOHAAe/bswdGjR3HmzBl07NgRBw4cAADs3LkTeXl5WLBgAcLDwwEAderUQf36\n9fHHH39Ar9ejf//+aNWqFdzd3dGkSRMMGTIEN27cQP/+/ZXXqa5dXzK+DN4COAAxxhh70NasWYN3\n3nkHgPj8+eWXX+Dj4wN/f3+8/PLL+PPPP9G6dWvMnz8fixcvBgB07NgRAHDgwAEMHToUt2/fBgA4\nOTkhKCgIkiRBp9Ohfv366NKlC5o0aQJALEuVnZ2trMxQeF216opbgCyAAxBjjLEHKS4uDmPGjIGv\nry8GDRqEkJAQvPTSS0hISMDs2bPRvXt3eHp6wsXFBY6OjlCpVDCbzQgJCYGLiwtSU1Nx8eJFPP/8\n81i2bBk8PT0RHR2NuLg4REVFITw8HCaTCYBYYiQkJARbtmzB6dOn4efnZ+V3bxncAmQBHIAYY4xZ\nmtw6k5ubi88++wy//PKLstTEhg0bYDabce3aNcyfPx+Ojo7w8/NDgwYNMH78eLRq1Qp+fn5wcnJS\n1q1TqVRISUmBvb095syZg6+++go7d+7Eo48+irNnzyIoKAgAEB0drdRB/nyrX78+dDodTp069SBP\nwX3FAcgCOAAxxhirjNu3byMiIgKvvfYaWrZsCXt7e7Rs2RLvvvsuiAjbtm3D448/jrlz54KI8Pnn\nn2P+/PmYMGECDh06hPj4eIwYMQK3bt3C999/j82bN2PcuHFo1KgR5syZg7y8PABAeno6nJ2dERsb\ni+HDh2PNmjXIzc1F8+bNMW/ePNStWxdOTk4AoCzYCgD16tWDp6cn/vnnH+W+6o67wCyAAxBjjLGK\nKj54+KOPPlLG8vj7+6NDhw7o378/tm7dipkzZ6JFixYICQlBREQE5s2bh5iYGOTn56NTp06oX78+\ndu3ahaioKLRq1Qomkwljx46Fv78/9Ho9hg0bhpdeegkajQYA4O7ujpYtWyoDogcPHoywsDC89NJL\nOHDgADQajTIwunAd69atC2dnZ1y8eBF5eXmwtbV9UKfrvuEAZAEcgBhjjJXl1q1b+OWXX7B9+3ZE\nRkaiT58+eOONN5SxNFlZWVCpVOjcuTMWLVoEb29vuLm5oVOnTnjsscdw48YN9O7dG9988w2Sk5Ox\ndOlStGzZEi1btsT169fxyCOP4OjRoxg2bJhy/0svvYTk5GR4eHgUqYu9vT2aN2+O9evXK7c1aNAA\nP/74I9544w3Ur18fdevWLfIYs9kMlUqF5ORk2NvbIz09vcTzVkccgCyAAxBjjLHiTCYT1Go1pk+f\njqVLlyIsLAwajQZLly7F8ePHMXPmTDz66KMIDw+H2WxG3759lTW2AODUqVMwGo3o2rUrXF1d0aVL\nF5w+fRrjx4+HWq0GALi5uaF27dr4559/EBAQoNwGoEhIiY+Ph6urK+zt7REaGorc3FxcuHABISEh\nMBqNqF27NpYvX15qy4782fb++++jVq1aNSL8ADwGyCI4ADHG2MNlw4YNWLRoEZKSkgCU/PtvNpuh\nVquxa9cuLF++HJ07d8bGjRtx5MgRfPPNN7hw4QJeeOEFnD9/Hs2aNUNISAgSEhIAiKuutm/fjqlT\np2LgwIFo0KABXF1dERwcDJVKhYkTJ2LcuHEAAJ1Oh4SEBMTExCAvLw8eHh44ePAgzp8/j6VLl+Lx\nxx+Hn58ffH19sXfvXgCAn58f3N3dsX//fgDiM4yIYGtrW+rnmFqtBhGhQ4cOCAkJKfV8GAwGy5zY\nB4gDkAVwAGKMsZpPngPn/PnzmDlzJhYtWoQLFy6UeqxKpYLJZMLvv/8OAJgyZQqCg4Ph6uqKZ555\nBlFRUfDw8MDrr7+uXJo+b948jB49Gj179sSQIUPg5OSEhQsXAgBsbW0RFBSEtLQ0nD9/Hmq1Gkaj\nEQAQEhKClJQUREREoHnz5vj000/RsWNHTJw4EX///TeGDh2Kw4cPo2/fvgBEC1FwcDD27NkDoOiY\npLLm9il8e2pqKrZt24b//e9/aNmyJby9vfHtt99W9vQ+cByALIADEGOM1SxXr17F2rVrMXr0aDRq\n1AgqlQrPPfccACAzMxOnT5+Gg4MDGjVqBKDs4LBx40aMGTMG/fv3L/IZ4evri9dffx1bt25FTEwM\nwsLCoFar8dVXX0GlUmHGjBnYv38/6tatq3SJBQYGws7ODn/++ScAKLeHhoaicePGsLGxgY+PD2xt\nbfHKK6/g888/x4YNG/D++++jbdu2ymt7eXnhhRdeUC6zl7vTCktOTsaCBQuQlZUFAEhLS8OkSZNQ\nt25duLm5YeDAgfjll18QGBiISZMm4fnnn6/U+bYGHgNkARyAGGOs+srJyUF2djbc3Nzw3Xff4bXX\nXkNSUhKcnJwQHByM9u3b45VXXkGvXr0AiEvCAeDy5cvl/t2/du0arl69iu+++w5AQUuLPKhYrVbD\nxcUFmZmZaNGiBVavXo1+/fph7dq1cHR0VK7ckj9jAgMDERISgt27d+Ppp59Wbm/fvj28vb0RGhqK\ndu3a4ZNPPlHGARVWOKSNHDlSCXSlhbdPPvkEn3/+OSRJwmuvvYZvv/0WS5Ysgb+/Pz7++GMEBQWh\ndu3acHd3h6Oj412f86qAA5AFjBs3Dn379kXDhg2tXRXGGGN3ITU1FRMnToS/vz8mT56MiRMnIjk5\nGZ9//jnatGlTZDZlOby4uLjA29sb8fHxyM7Ohqura5FuJHl/9+7daNeuHRwcHJTXK7yExMaNG9Gh\nQwc0btwYcXFx0Gg0aNGihfJ8MnnOnZCQEEycOBGbNm0CAGXAspubmxJ47OzslMfJkyaWNmePjY0N\nkpOTQURISUnBoUOHMGrUKABirqDVq1cjOzsb8+fPx8iRI3Hr1i3k5+dj9uzZGDhwoEXOvbVxF5gF\nPPHEE5gwYYIyAp8xxljVcO7cOaxbtw5ZWVm4du0aIiMji9zv4uKC69ev4+LFi9BqtahVqxY6d+6M\nESNGoGXLlvD19YVery/R0l+/fn0AwMmTJwEUhA2goGvq5s2bytgdIlJafVQqFWbOnIk//vgDEyZM\nACDW2goNDVWev6yWpX79+mHVqlUVeu8qlapE+JGf9/Tp0/D09MScOXOwb98+TJ8+HYMHD8bt27ex\nceNGBAQEYPr06UhPT8fmzZsxfvx4eHt7IzIyEmlpadi+fTsmT56MV199VZk5urr1gnAAYowxViPc\nunULmzdvxquvvormzZtDp9OhUaNGmDp1KgwGA3788Uc8++yzuHTpEgAog4hv3boFo9EIOzs7DBgw\nANeuXcOmTZvw008/4cUXX0Tjxo2xaNEiAGK8TE5ODkJDQwEAV65cASACkBwAbGxE50pubi5iYmIA\nANnZ2UhMTMTBgwfRp08fzJw5E1OmTMGjjz4KQHSr1a5dG8nJyTCZTOXOtFw4bJWHiJQik4PcpUuX\n4Obmhr///hsZGRmws7PDxo0bsXjxYmzevBlhYWGYMGECGjVqhN27d8PFxQVt2rTBxx9/DD8/Pzz1\n1FPYuXMnWrVqBS8vryLPXV1wFxhjjLFqS+5u2rJlC4YPH46srCw0adIEYWFhePbZZxESEoLmzZvD\nw8MDLi4uuHTpEtavX48333xTGfwbFBSEzMxMAEDLli0xf/58vPTSS/Dy8oKLiwtGjRqFUaNGKa9l\nMplgNpuh1WqV7i15vI58lZarqysaNGiAc+fOYcCAAfD19cXp06eh1WoREBCApKQkODs7AyiYaHDM\nmDFwdnZW5g8qS0WXoSgcSG7evIm//voL+/btQ3R0NB599FE888wz2Lx5MwYOHIjw8HBcuHAB77//\nPsxmMzZs2AA3NzeMGTMG06dPByAuda9fvz7OnTuH69evQ6vVwt3d/S7/xaoODkCMMcaqPDl8zJ07\nFwsXLsTmzZvRqlUrSJKEzMxMrFy5EllZWZgzZw6effZZ6HQ6ODk5Ka0xANClSxd06dIFq1evRseO\nHdGpUydl7hy5pcbe3h65ubn48ccfMWjQoFJnU9bpdFCpVMjLy1MCw59//omDBw8iOjoaubm5mD17\nNsaNG4cVK1Zg/fr12L59O0JDQzF8+HD06NGjSNCRA82QIUMscq5yc3MRGRmJo0eP4vDhw4iOjsa5\nc+eg1WpRt25d9OnTB507d8aGDRuQlpaG/Px89O7dG+vXr8fkyZNx6NAh5ZL57t274/nnn1danQIC\nApCTk6O0+siq4/IYHIAYY4xVeXJrhkqlgpeXF44ePQoAaNasGRwdHdGuXTucOHEC3bt3Vz6ciQiJ\niYnw9PQEESEkJARz5sxBu3btMHXqVOzfvx8GgwHHjx9XXicwMBAAUKtWLQBFZ1OOi4uDu7s77Ozs\n0LJlSzg4OCAmJgZjxowBALRr1w4LFy6En58fQkJCoNfr8fzzz2PkyJFFgpisvFaee2U2mzF16lTM\nmzcPDRo0QK9evbB161Y0btwYP/zwA1xdXeHm5gaNRgN/f3/s27cPHh4eCA0NhSRJePrpp/H2229D\np9MBAG7cuAEAuHjxIkJDQ/H777/j5s2b0Ol02LJlCw4ePIht27Zh8ODB+OSTT6rVFWEcgBhjjFVZ\nubm5uHjxIrKysrBjxw6kpqYiMjIShw8fxocffghHR0e8/fbbeOaZZ/D+++/jo48+QnBwMPbt24d/\n/vkH7du3x/bt25VWljZt2mDMmDFYsWIFvv76awwfPhyOjo6oV68eDAYDPDw8oNfr8ccff6BOnTrY\ntWsXtm3bhpMnT+LGjRvYuXMnevTogf/+97+wtbXF/PnzYTabMWHCBPTs2RM+Pj4lgo2NjQ3MZrMy\ntqfwKuuWJg+wfuuttyBJElxdXXH58mXY2NjA39+/yFViDRo0ACCuHPP394enpycOHDiAFi1aKN1y\nHh4eqFu3Lnbs2IHQ0FAsWrQIvXr1QkJCAjIzMxEWFoaxY8fimWeeqVbhB+BB0IwxxqowtVqNBQsW\noF27dpgzZ44SLtq3b48ffvgBly9fRr9+/aDVaqHX6/H777/jwoULaNSoEdauXYuIiIgSY2amTZuG\nZs2aYdKkSbhx4wb8/f3h6uoKOzs73Lx5E4GBgZg2bRrat2+PyZMn4+LFixg1ahQiIyPRo0cPmEwm\n2Nra4r///S/++usvREVF4fnnn4e/v79Sv+JXRKlUKmg0GqjV6vs+WFin08HNzQ2urq5BzYWyAAAD\n40lEQVQAxBVmSUlJuHz5MgDRXQVAaSk7fPgwnJycEBYWhoMHDxapv5ubG5o3b46///5bWQbD29sb\no0ePxtKlS/HNN99g5syZaNy48X19T/cDtwAxxhirsuTZjfV6PTZs2ICePXsq44HkhUVbtGiBlStX\nAgAmTpyId955p8yBxGazGT4+Pvjoo48wduxYdOzYEfHx8crcNk5OTvD398eVK1cwefJk1KlTB61a\ntULDhg2h1WoBlOy6kq8mKxxuqsIVUUajETY2NmjatCl2796NS5cuFQkqLi4uAIC//voLarUajzzy\nCFavXg2gYKC1s7MzZsyYAbVajdDQUGzZsgUdOnSo1oOfZRyAGGOMVWmtWrWCnZ0dXF1dS0w2WLdu\nXYSFheHy5csIDAwssbxDYmIiTp48CQcHB3Ts2FHp2unduzfmzp2Lp59+ushl5XXq1MGSJUvg6Oio\nXKV1J6WN76kK5BDTtGlTmM1mnD9/HkDBuZHH+cjrmY0ePRqrVq1CRkYG9Hq9cmyrVq2U5+zfv7+y\nX95Ei9VB9aw1Y4yxh4Y8QeDixYuV7huZs7MzGjZsiH/++QcNGzbE9u3b8e2332Ly5MkIDg5G3bp1\n8eSTT+LQoUMAioaVJ598EnPmzAEAdOjQAYC4nN3b2xvOzs4gIpjN5grPu1PVyMEkNDQUzs7OuHjx\nIoCCAKTX62FnZ4fExESkpaXB19cXZ8+eVcJPaQrPd1TaRIvVSdWMrYwxxti/AgIC0L9/f/z+++9I\nSkqCt7c3zGYz1Go1srOzcenSJbRv3x7p6ek4d+4cXnjhBbi4uMDLywsrVqxA27ZtSx2jolKpMGnS\nJAwaNAj+/v4l7pckqUp0ZVUGEcHR0RF2dna4evUqUlNTla6vunXrws/PDxcvXsT169fh7OwMe3v7\ncp+vOgee4jgAMcYYq9I0Gg2GDRuG1atXY8eOHRg9ejTUajVMJhM2b96MlJQUvP766/jpp59w9OhR\njB8/Hk2aNEFAQAD8/PwgSVKRtbpkcsApLfzUFHJQVKlUcHNzK9IC5unpibCwMHTo0AHe3t5WrKV1\nSFTdFu9gjDH20ElKSkJAQADatGmDuXPnIiMjA5s2bcLatWuxbNkyDB06FPn5+cqMzEyQg19SUpIy\ntxETOAAxxhirFlq3bo2zZ8/C3d0dsbGx0Gg0WL58OUaOHFnkOHnF9ZrQhcXun5rTmccYY6xGa9Wq\nFbKzs+Ht7Y3o6GhkZ2cr4afwd3m5y4fDDysPjwFijDFWLTz55JM4d+4c3n77bWVSPrmLh8MOu1vc\nBcYYY4yxhw53gTHGGGPsocMBiDHGGGMPnf8HjCsser0d4bQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image at 0x7f70a37df110>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(filename='images/optimizing-what.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fortran"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### F2PY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "F2PY is a program that (almost) automatically wraps fortran code for use in Python: By using the `f2py` program we can compile fortran code into a module that we can import in a Python program.\n",
    "\n",
    "F2PY is a part of NumPy, but you will also need to have a fortran compiler to run the examples below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 0: scalar input, no output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting hellofortran.f\n"
     ]
    }
   ],
   "source": [
    "%%file hellofortran.f\n",
    "C File  hellofortran.f\n",
    "        subroutine hellofortran (n)\n",
    "        integer n\n",
    "       \n",
    "        do 100 i=0, n\n",
    "            print *, \"Fortran says hello\"\n",
    "100     continue\n",
    "        end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a python module using `f2py`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[39mrunning build\u001b[0m\n",
      "\u001b[39mrunning config_cc\u001b[0m\n",
      "\u001b[39munifing config_cc, config, build_clib, build_ext, build commands --compiler options\u001b[0m\n",
      "\u001b[39mrunning config_fc\u001b[0m\n",
      "\u001b[39munifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\u001b[0m\n",
      "\u001b[39mrunning build_src\u001b[0m\n",
      "\u001b[39mbuild_src\u001b[0m\n",
      "\u001b[39mbuilding extension \"hellofortran\" sources\u001b[0m\n",
      "\u001b[39mf2py options: []\u001b[0m\n",
      "\u001b[39mf2py:> /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpz2IPjB/src.linux-x86_64-2.7\u001b[0m\n",
      "Reading fortran codes...\n",
      "\tReading file 'hellofortran.f' (format:fix,strict)\n",
      "Post-processing...\n",
      "\tBlock: hellofortran\n",
      "\t\t\tBlock: hellofortran\n",
      "Post-processing (stage 2)...\n",
      "Building modules...\n",
      "\tBuilding module \"hellofortran\"...\n",
      "\t\tConstructing wrapper function \"hellofortran\"...\n",
      "\t\t  hellofortran(n)\n",
      "\tWrote C/API module \"hellofortran\" to file \"/tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c\"\n",
      "\u001b[39m  adding '/tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c' to sources.\u001b[0m\n",
      "\u001b[39m  adding '/tmp/tmpz2IPjB/src.linux-x86_64-2.7' to include_dirs.\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpz2IPjB/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpz2IPjB/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mbuild_src: building npy-pkg config files\u001b[0m\n",
      "\u001b[39mrunning build_ext\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler using build_ext\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mFound executable /usr/bin/gfortran\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler using build_ext\u001b[0m\n",
      "\u001b[39mbuilding 'hellofortran' extension\u001b[0m\n",
      "\u001b[39mcompiling C sources\u001b[0m\n",
      "\u001b[39mC compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC\n",
      "\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpz2IPjB/tmp\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpz2IPjB/tmp/tmpz2IPjB\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpz2IPjB/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.c:17:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.c:2:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "\u001b[39mcompiling Fortran sources\u001b[0m\n",
      "\u001b[39mFortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpz2IPjB/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mgfortran:f77: hellofortran.f\u001b[0m\n",
      "\u001b[39m/usr/bin/gfortran -Wall -Wall -shared /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7/hellofortranmodule.o /tmp/tmpz2IPjB/tmp/tmpz2IPjB/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpz2IPjB/hellofortran.o -lgfortran -o ./hellofortran.so\u001b[0m\n",
      "Removing build directory /tmp/tmpz2IPjB\n"
     ]
    }
   ],
   "source": [
    "!f2py -c -m hellofortran hellofortran.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example of a python script that use the module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting hello.py\n"
     ]
    }
   ],
   "source": [
    "%%file hello.py\n",
    "import hellofortran\n",
    "\n",
    "hellofortran.hellofortran(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Fortran says hello\r\n",
      " Fortran says hello\r\n",
      " Fortran says hello\r\n",
      " Fortran says hello\r\n",
      " Fortran says hello\r\n",
      " Fortran says hello\r\n"
     ]
    }
   ],
   "source": [
    "# run the script\n",
    "!python hello.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: vector input and scalar output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting dprod.f\n"
     ]
    }
   ],
   "source": [
    "%%file dprod.f\n",
    "\n",
    "       subroutine dprod(x, y, n)\n",
    "    \n",
    "       double precision x(n), y\n",
    "       y = 1.0\n",
    "    \n",
    "       do 100 i=1, n\n",
    "           y = y * x(i)\n",
    "100    continue\n",
    "       end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading fortran codes...\r\n",
      "\tReading file 'dprod.f' (format:fix,strict)\r\n",
      "Post-processing...\r\n",
      "\tBlock: dprod\r\n",
      "{}\r\n",
      "In: :dprod:dprod.f:dprod\r\n",
      "vars2fortran: No typespec for argument \"n\".\r\n",
      "\t\t\tBlock: dprod\r\n",
      "Post-processing (stage 2)...\r\n",
      "Saving signatures to file \"./dprod.pyf\"\r\n"
     ]
    }
   ],
   "source": [
    "!rm -f dprod.pyf\n",
    "!f2py -m dprod -h dprod.pyf dprod.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `f2py` program generated a module declaration file called `dsum.pyf`. Let's look what's in it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "!    -*- f90 -*-\r\n",
      "! Note: the context of this file is case sensitive.\r\n",
      "\r\n",
      "python module dprod ! in \r\n",
      "    interface  ! in :dprod\r\n",
      "        subroutine dprod(x,y,n) ! in :dprod:dprod.f\r\n",
      "            double precision dimension(n) :: x\r\n",
      "            double precision :: y\r\n",
      "            integer, optional,check(len(x)>=n),depend(x) :: n=len(x)\r\n",
      "        end subroutine dprod\r\n",
      "    end interface \r\n",
      "end python module dprod\r\n",
      "\r\n",
      "! This file was auto-generated with f2py (version:2).\r\n",
      "! See http://cens.ioc.ee/projects/f2py2e/\r\n"
     ]
    }
   ],
   "source": [
    "!cat dprod.pyf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The module does not know what Fortran subroutine arguments is input and output, so we need to manually edit the module declaration files and mark output variables with `intent(out)` and input variable with `intent(in)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting dprod.pyf\n"
     ]
    }
   ],
   "source": [
    "%%file dprod.pyf\n",
    "python module dprod ! in \n",
    "    interface  ! in :dprod\n",
    "        subroutine dprod(x,y,n) ! in :dprod:dprod.f\n",
    "            double precision dimension(n), intent(in) :: x\n",
    "            double precision, intent(out) :: y\n",
    "            integer, optional,check(len(x)>=n),depend(x),intent(in) :: n=len(x)\n",
    "        end subroutine dprod\n",
    "    end interface \n",
    "end python module dprod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compile the fortran code into a module that can be included in python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[39mrunning build\u001b[0m\n",
      "\u001b[39mrunning config_cc\u001b[0m\n",
      "\u001b[39munifing config_cc, config, build_clib, build_ext, build commands --compiler options\u001b[0m\n",
      "\u001b[39mrunning config_fc\u001b[0m\n",
      "\u001b[39munifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\u001b[0m\n",
      "\u001b[39mrunning build_src\u001b[0m\n",
      "\u001b[39mbuild_src\u001b[0m\n",
      "\u001b[39mbuilding extension \"dprod\" sources\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpWyCvx1/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mf2py options: []\u001b[0m\n",
      "\u001b[39mf2py: dprod.pyf\u001b[0m\n",
      "Reading fortran codes...\n",
      "\tReading file 'dprod.pyf' (format:free)\n",
      "Post-processing...\n",
      "\tBlock: dprod\n",
      "\t\t\tBlock: dprod\n",
      "Post-processing (stage 2)...\n",
      "Building modules...\n",
      "\tBuilding module \"dprod\"...\n",
      "\t\tConstructing wrapper function \"dprod\"...\n",
      "\t\t  y = dprod(x,[n])\n",
      "\tWrote C/API module \"dprod\" to file \"/tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c\"\n",
      "\u001b[39m  adding '/tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c' to sources.\u001b[0m\n",
      "\u001b[39m  adding '/tmp/tmpWyCvx1/src.linux-x86_64-2.7' to include_dirs.\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpWyCvx1/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpWyCvx1/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mbuild_src: building npy-pkg config files\u001b[0m\n",
      "\u001b[39mrunning build_ext\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler using build_ext\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mFound executable /usr/bin/gfortran\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler using build_ext\u001b[0m\n",
      "\u001b[39mbuilding 'dprod' extension\u001b[0m\n",
      "\u001b[39mcompiling C sources\u001b[0m\n",
      "\u001b[39mC compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC\n",
      "\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpWyCvx1/tmp\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpWyCvx1/tmp/tmpWyCvx1\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpWyCvx1/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c:18:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "/tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]\n",
      " static int f2py_size(PyArrayObject* var, ...)\n",
      "            ^\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.c:2:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "\u001b[39mcompiling Fortran sources\u001b[0m\n",
      "\u001b[39mFortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpWyCvx1/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mgfortran:f77: dprod.f\u001b[0m\n",
      "\u001b[39m/usr/bin/gfortran -Wall -Wall -shared /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7/dprodmodule.o /tmp/tmpWyCvx1/tmp/tmpWyCvx1/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpWyCvx1/dprod.o -lgfortran -o ./dprod.so\u001b[0m\n",
      "Removing build directory /tmp/tmpWyCvx1\n"
     ]
    }
   ],
   "source": [
    "!f2py -c dprod.pyf dprod.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using the module from Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import dprod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on module dprod:\n",
      "\n",
      "NAME\n",
      "    dprod\n",
      "\n",
      "FILE\n",
      "    /home/rob/Desktop/scientific-python-lectures/dprod.so\n",
      "\n",
      "DESCRIPTION\n",
      "    This module 'dprod' is auto-generated with f2py (version:2).\n",
      "    Functions:\n",
      "      y = dprod(x,n=len(x))\n",
      "    .\n",
      "\n",
      "DATA\n",
      "    __version__ = '$Revision: $'\n",
      "    dprod = <fortran object>\n",
      "\n",
      "VERSION\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(dprod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.082818640342675e+62"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dprod.dprod(arange(1,50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.0828186403426752e+62"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compare to numpy\n",
    "prod(arange(1.0,50.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dprod.dprod(arange(1,10), 5) # only the 5 first elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare performance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xvec = rand(500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000000 loops, best of 3: 882 ns per loop\n"
     ]
    }
   ],
   "source": [
    "timeit dprod.dprod(xvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 loops, best of 3: 4.45 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit xvec.prod()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: cummulative sum, vector input and vector output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cummulative sum function for an array of data is a good example of a loop intense algorithm: Loop through a vector and store the cummulative sum in another vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# simple python algorithm: example of a SLOW implementation\n",
    "# Why? Because the loop is implemented in python.\n",
    "def py_dcumsum(a):\n",
    "    b = empty_like(a)\n",
    "    b[0] = a[0]\n",
    "    for n in range(1,len(a)):\n",
    "        b[n] = b[n-1]+a[n]\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fortran subroutine for the same thing: here we have added the `intent(in)` and `intent(out)` as comment lines in the original fortran code, so we do not need to manually edit the fortran module declaration file generated by `f2py`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting dcumsum.f\n"
     ]
    }
   ],
   "source": [
    "%%file dcumsum.f\n",
    "c File dcumsum.f\n",
    "       subroutine dcumsum(a, b, n)\n",
    "       double precision a(n)\n",
    "       double precision b(n)\n",
    "       integer n\n",
    "cf2py  intent(in) :: a\n",
    "cf2py  intent(out) :: b\n",
    "cf2py  intent(hide) :: n\n",
    "\n",
    "       b(1) = a(1)\n",
    "       do 100 i=2, n\n",
    "           b(i) = b(i-1) + a(i)\n",
    "100    continue\n",
    "       end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can directly compile the fortran code to a python module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[39mrunning build\u001b[0m\n",
      "\u001b[39mrunning config_cc\u001b[0m\n",
      "\u001b[39munifing config_cc, config, build_clib, build_ext, build commands --compiler options\u001b[0m\n",
      "\u001b[39mrunning config_fc\u001b[0m\n",
      "\u001b[39munifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\u001b[0m\n",
      "\u001b[39mrunning build_src\u001b[0m\n",
      "\u001b[39mbuild_src\u001b[0m\n",
      "\u001b[39mbuilding extension \"dcumsum\" sources\u001b[0m\n",
      "\u001b[39mf2py options: []\u001b[0m\n",
      "\u001b[39mf2py:> /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpfvrMl6/src.linux-x86_64-2.7\u001b[0m\n",
      "Reading fortran codes...\n",
      "\tReading file 'dcumsum.f' (format:fix,strict)\n",
      "Post-processing...\n",
      "\tBlock: dcumsum\n",
      "\t\t\tBlock: dcumsum\n",
      "Post-processing (stage 2)...\n",
      "Building modules...\n",
      "\tBuilding module \"dcumsum\"...\n",
      "\t\tConstructing wrapper function \"dcumsum\"...\n",
      "\t\t  b = dcumsum(a)\n",
      "\tWrote C/API module \"dcumsum\" to file \"/tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c\"\n",
      "\u001b[39m  adding '/tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c' to sources.\u001b[0m\n",
      "\u001b[39m  adding '/tmp/tmpfvrMl6/src.linux-x86_64-2.7' to include_dirs.\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpfvrMl6/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcopying /usr/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpfvrMl6/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mbuild_src: building npy-pkg config files\u001b[0m\n",
      "\u001b[39mrunning build_ext\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler\u001b[0m\n",
      "\u001b[39mcustomize UnixCCompiler using build_ext\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mFound executable /usr/bin/gfortran\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n",
      "\u001b[39mcustomize Gnu95FCompiler using build_ext\u001b[0m\n",
      "\u001b[39mbuilding 'dcumsum' extension\u001b[0m\n",
      "\u001b[39mcompiling C sources\u001b[0m\n",
      "\u001b[39mC compiler: x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC\n",
      "\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpfvrMl6/tmp\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpfvrMl6/tmp/tmpfvrMl6\u001b[0m\n",
      "\u001b[39mcreating /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpfvrMl6/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c:18:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "/tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]\n",
      " static int f2py_size(PyArrayObject* var, ...)\n",
      "            ^\n",
      "\u001b[39mx86_64-linux-gnu-gcc: /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c\u001b[0m\n",
      "In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
      "                 from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.h:13,\n",
      "                 from /tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.c:2:\n",
      "/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "\u001b[39mcompiling Fortran sources\u001b[0m\n",
      "\u001b[39mFortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n",
      "\u001b[39mcompile options: '-I/tmp/tmpfvrMl6/src.linux-x86_64-2.7 -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'\u001b[0m\n",
      "\u001b[39mgfortran:f77: dcumsum.f\u001b[0m\n",
      "\u001b[39m/usr/bin/gfortran -Wall -Wall -shared /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7/dcumsummodule.o /tmp/tmpfvrMl6/tmp/tmpfvrMl6/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpfvrMl6/dcumsum.o -lgfortran -o ./dcumsum.so\u001b[0m\n",
      "Removing build directory /tmp/tmpfvrMl6\n"
     ]
    }
   ],
   "source": [
    "!f2py -c dcumsum.f -m dcumsum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import dcumsum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py_dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dcumsum.dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cumsum(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benchmark the different implementations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = rand(10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 4.83 ms per loop\n"
     ]
    }
   ],
   "source": [
    "timeit py_dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100000 loops, best of 3: 12.2 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit dcumsum.dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 27.4 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit a.cumsum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. http://www.scipy.org/F2py\n",
    "1. http://dsnra.jpl.nasa.gov/software/Python/F2PY_tutorial.pdf\n",
    "1. http://www.shocksolution.com/2009/09/f2py-binding-fortran-python/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ctypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ctypes is a Python library for calling out to C code. It is not as automatic as `f2py`, and we manually need to load the library and set properties such as the functions return and argument types. On the otherhand we do not need to touch the C code at all. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting functions.c\n"
     ]
    }
   ],
   "source": [
    "%%file functions.c\n",
    "\n",
    "#include <stdio.h>\n",
    "\n",
    "void hello(int n);\n",
    "\n",
    "double dprod(double *x, int n);\n",
    "\n",
    "void dcumsum(double *a, double *b, int n);\n",
    "\n",
    "void\n",
    "hello(int n)\n",
    "{\n",
    "    int i;\n",
    "    \n",
    "    for (i = 0; i < n; i++)\n",
    "    {\n",
    "        printf(\"C says hello\\n\");\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "double \n",
    "dprod(double *x, int n)\n",
    "{\n",
    "    int i;\n",
    "    double y = 1.0;\n",
    "    \n",
    "    for (i = 0; i < n; i++)\n",
    "    {\n",
    "        y *= x[i];\n",
    "    }\n",
    "\n",
    "    return y;\n",
    "}\n",
    "\n",
    "void\n",
    "dcumsum(double *a, double *b, int n)\n",
    "{\n",
    "    int i;\n",
    "    \n",
    "    b[0] = a[0];\n",
    "    for (i = 1; i < n; i++)\n",
    "    {\n",
    "        b[i] = a[i] + b[i-1];\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compile the C file into a shared library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "!gcc -c -Wall -O2 -Wall -ansi -pedantic -fPIC -o functions.o functions.c\n",
    "!gcc -o libfunctions.so -shared functions.o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a compiled shared library `libfunctions.so`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "libfunctions.so: ELF 64-bit LSB  shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=d68173ae6a804f703472af96f413b81a189db4b8, not stripped\r\n"
     ]
    }
   ],
   "source": [
    "!file libfunctions.so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to write wrapper functions to access the C library: To load the library we use the ctypes package, which included in the Python standard library (with extensions from numpy for passing arrays to C). Then we manually set the types of the argument and return values (no automatic code inspection here!). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting functions.py\n"
     ]
    }
   ],
   "source": [
    "%%file functions.py\n",
    "\n",
    "import numpy\n",
    "import ctypes\n",
    "\n",
    "_libfunctions = numpy.ctypeslib.load_library('libfunctions', '.')\n",
    "\n",
    "_libfunctions.hello.argtypes = [ctypes.c_int]\n",
    "_libfunctions.hello.restype  =  ctypes.c_void_p\n",
    "\n",
    "_libfunctions.dprod.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]\n",
    "_libfunctions.dprod.restype  = ctypes.c_double\n",
    "\n",
    "_libfunctions.dcumsum.argtypes = [numpy.ctypeslib.ndpointer(dtype=numpy.float), numpy.ctypeslib.ndpointer(dtype=numpy.float), ctypes.c_int]\n",
    "_libfunctions.dcumsum.restype  = ctypes.c_void_p\n",
    "\n",
    "def hello(n):\n",
    "    return _libfunctions.hello(int(n))\n",
    "\n",
    "def dprod(x, n=None):\n",
    "    if n is None:\n",
    "        n = len(x)\n",
    "    x = numpy.asarray(x, dtype=numpy.float)\n",
    "    return _libfunctions.dprod(x, int(n))\n",
    "\n",
    "def dcumsum(a, n):\n",
    "    a = numpy.asarray(a, dtype=numpy.float)\n",
    "    b = numpy.empty(len(a), dtype=numpy.float)\n",
    "    _libfunctions.dcumsum(a, b, int(n))\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting run_hello_c.py\n"
     ]
    }
   ],
   "source": [
    "%%file run_hello_c.py\n",
    "\n",
    "import functions\n",
    "\n",
    "functions.hello(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C says hello\r\n",
      "C says hello\r\n",
      "C says hello\r\n"
     ]
    }
   ],
   "source": [
    "!python run_hello_c.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Product function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120.0"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "functions.dprod([1,2,3,4,5]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cummulative sum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = rand(100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res_c = functions.dcumsum(a, len(a)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res_fortran = dcumsum.dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0., ...,  0.,  0.,  0.])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_c - res_fortran"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple benchmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 286 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit functions.dcumsum(a, len(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 3: 119 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit dcumsum.dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 261 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit a.cumsum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* http://docs.python.org/2/library/ctypes.html\n",
    "* http://www.scipy.org/Cookbook/Ctypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A hybrid between python and C that can be compiled: Basically Python code with type declarations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting cy_dcumsum.pyx\n"
     ]
    }
   ],
   "source": [
    "%%file cy_dcumsum.pyx\n",
    "\n",
    "cimport numpy\n",
    "\n",
    "def dcumsum(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):\n",
    "    cdef int i, n = len(a)\n",
    "    b[0] = a[0]\n",
    "    for i from 1 <= i < n:\n",
    "        b[i] = b[i-1] + a[i]\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A build file for generating C code and compiling it into a Python module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting setup.py\n"
     ]
    }
   ],
   "source": [
    "%%file setup.py\n",
    "\n",
    "from distutils.core import setup\n",
    "from distutils.extension import Extension\n",
    "from Cython.Distutils import build_ext\n",
    "\n",
    "setup(\n",
    "    cmdclass = {'build_ext': build_ext},\n",
    "    ext_modules = [Extension(\"cy_dcumsum\", [\"cy_dcumsum.pyx\"])]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running build_ext\n",
      "cythoning cy_dcumsum.pyx to cy_dcumsum.c\n",
      "warning: /usr/local/lib/python2.7/dist-packages/Cython/Includes/numpy.pxd:869:17: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.\n",
      "warning: /usr/local/lib/python2.7/dist-packages/Cython/Includes/numpy.pxd:869:24: Non-trivial type declarators in shared declaration (e.g. mix of pointers and values). Each pointer declaration should be on its own line.\n",
      "building 'cy_dcumsum' extension\n",
      "x86_64-linux-gnu-gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/include/python2.7 -c cy_dcumsum.c -o build/temp.linux-x86_64-2.7/cy_dcumsum.o\n",
      "In file included from /usr/include/python2.7/numpy/ndarraytypes.h:1761:0,\n",
      "                 from /usr/include/python2.7/numpy/ndarrayobject.h:17,\n",
      "                 from /usr/include/python2.7/numpy/arrayobject.h:4,\n",
      "                 from cy_dcumsum.c:352:\n",
      "/usr/include/python2.7/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning \"Using deprecated NumPy API, disable it by \" \"#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
      " #warning \"Using deprecated NumPy API, disable it by \" \\\n",
      "  ^\n",
      "In file included from /usr/include/python2.7/numpy/ndarrayobject.h:26:0,\n",
      "                 from /usr/include/python2.7/numpy/arrayobject.h:4,\n",
      "                 from cy_dcumsum.c:352:\n",
      "/usr/include/python2.7/numpy/__multiarray_api.h:1629:1: warning: ‘_import_array’ defined but not used [-Wunused-function]\n",
      " _import_array(void)\n",
      " ^\n",
      "In file included from /usr/include/python2.7/numpy/ufuncobject.h:327:0,\n",
      "                 from cy_dcumsum.c:353:\n",
      "/usr/include/python2.7/numpy/__ufunc_api.h:241:1: warning: ‘_import_umath’ defined but not used [-Wunused-function]\n",
      " _import_umath(void)\n",
      " ^\n",
      "x86_64-linux-gnu-gcc -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -D_FORTIFY_SOURCE=2 -g -fstack-protector --param=ssp-buffer-size=4 -Wformat -Werror=format-security build/temp.linux-x86_64-2.7/cy_dcumsum.o -o /home/rob/Desktop/scientific-python-lectures/cy_dcumsum.so\n"
     ]
    }
   ],
   "source": [
    "!python setup.py build_ext --inplace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cy_dcumsum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = array([1,2,3,4], dtype=float)\n",
    "b = empty_like(a)\n",
    "cy_dcumsum.dcumsum(a,b)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = empty_like(a)\n",
    "cy_dcumsum.dcumsum(a, b)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.,   3.,   6.,  10.,  15.,  21.,  28.,  36.])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py_dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = rand(100000)\n",
    "b = empty_like(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 50.1 ms per loop\n"
     ]
    }
   ],
   "source": [
    "timeit py_dcumsum(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 263 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit cy_dcumsum.dcumsum(a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython in the IPython notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working with the IPython (especially in the notebook), there is a more convenient way of compiling and loading Cython code. Using the `%%cython` IPython magic (command to IPython), we can simply type the Cython code in a code cell and let IPython take care of the conversion to C code, compilation and loading of the function. To be able to use the `%%cython` magic, we first need to load the extension `cythonmagic`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext cythonmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "cimport numpy\n",
    "\n",
    "def cy_dcumsum2(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):\n",
    "    cdef int i, n = len(a)\n",
    "    b[0] = a[0]\n",
    "    for i from 1 <= i < n:\n",
    "        b[i] = b[i-1] + a[i]\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 265 µs per loop\n"
     ]
    }
   ],
   "source": [
    "timeit cy_dcumsum2(a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* http://cython.org\n",
    "* http://docs.cython.org/src/userguide/tutorial.html\n",
    "* http://wiki.cython.org/tutorials/numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "application/json": {
       "Software versions": [
        {
         "module": "Python",
         "version": "2.7.6 (default, Mar 22 2014, 22:59:56) [GCC 4.8.2]"
        },
        {
         "module": "IPython",
         "version": "1.1.0"
        },
        {
         "module": "OS",
         "version": "posix [linux2]"
        },
        {
         "module": "ctypes",
         "version": "1.1.0"
        },
        {
         "module": "Cython",
         "version": "0.20.2"
        }
       ]
      },
      "text/html": [
       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>2.7.6 (default, Mar 22 2014, 22:59:56) [GCC 4.8.2]</td></tr><tr><td>IPython</td><td>1.1.0</td></tr><tr><td>OS</td><td>posix [linux2]</td></tr><tr><td>ctypes</td><td>1.1.0</td></tr><tr><td>Cython</td><td>0.20.2</td></tr><tr><td colspan='2'>Tue Aug 26 23:37:29 2014 JST</td></tr></table>"
      ],
      "text/latex": [
       "\\begin{table}\n",
       "\\begin{tabular}{|l|l|}\\hline\n",
       "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
       "Python & 2.7.6 (default, Mar 22 2014, 22:59:56) [GCC 4.8.2] \\\\ \\hline\n",
       "IPython & 1.1.0 \\\\ \\hline\n",
       "OS & posix [linux2] \\\\ \\hline\n",
       "ctypes & 1.1.0 \\\\ \\hline\n",
       "Cython & 0.20.2 \\\\ \\hline\n",
       "\\hline \\multicolumn{2}{|l|}{Tue Aug 26 23:37:29 2014 JST} \\\\ \\hline\n",
       "\\end{tabular}\n",
       "\\end{table}\n"
      ],
      "text/plain": [
       "Software versions\n",
       "Python 2.7.6 (default, Mar 22 2014, 22:59:56) [GCC 4.8.2]\n",
       "IPython 1.1.0\n",
       "OS posix [linux2]\n",
       "ctypes 1.1.0\n",
       "Cython 0.20.2\n",
       "\n",
       "Tue Aug 26 23:37:29 2014 JST"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%reload_ext version_information\n",
    "\n",
    "%version_information ctypes, Cython"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}