{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using Fortran and C code with Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "this is another notebook from J.R. Johansson\n", "\n", "here is a link to [another F2Py tutorial](https://github.com/thehackerwithin/PyTrieste/wiki/F2Py)" ] }, { "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\n", "\n", "#if you are not on byrd: sudo yum install f2py" ] }, { "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": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAEgCAYAAABchszxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdclXX7B/DP9xw4bMUByhFEwgUOVBy5QkX0lyMtZ5qr1BwVj6OsnifBtByplZapmJjkyp0jUBFHuMUJiQuVIYgsRTbn+v1xd249MkQ4cESv9+t1Xh3u9b3umxPn8jsFEREYY4wxxl4hCkMHwBhjjDFW0TgBYowxxtgrhxMgxhhjjL1yOAFijLEKkp6ejgcPHhg6DMYYOAFijL2AsrOz8ejRI0OHoTcajQZxcXFo0aIFXF1dcefOHQBAcnIy5syZg0uXLhk4QsZePZwAMfYKS0hIgJeXF65du2boUAAAWVlZiI+PR//+/eHo6IgtW7YYNJ4bN26gTp06MDc3x4EDB0p1jXnz5sHIyAj29vaYOHEivLy8UK9ePXzwwQfYs2cPZs6cidGjR+s3cMbYM3ECxFgFe/jwIQYMGIDmzZsjPDxcZ9+6deugUCgKvNRqNQICAp67rJ07d0KhUGDjxo3YvHkzPv30U3kfEWHgwIEIDg6WE434+HiYmZlh4sSJSExMRPXq1QuNx9zcHGfPni2yrKioKCiVSqxcubJATKNHj4ZCocDRo0cBAJmZmZg1axaUSiXMzc2hVqvx4MEDPHr0CLt37y7RfXbq1AmfffZZofsyMjIwefJknD9/HnFxcWjatGmhxz148AC3bt2Sf46Pj0ePHj2QlpaGvLw8eHt7IyMjo0TxPOmzzz7D3LlzsX//fkybNg1Tp04F8Pj516hRA+fPn39hklDGXhVGhg6AsVdJZGQkRo0ahVOnTgEAunbtiv3798PNzQ0A5C9BR0dHTJgwATY2NsjJycG3336LUaNGwdbWFl26dEFKSkqx5VhaWsLS0hIHDx6Uy01JScGSJUugVCoxb948/Pzzzzh16hQUCgW2bNmCL774AsePH0d2djYiIyNhbW2NpUuXIjs7W75ubm4u5syZg9jYWCxevBjr1q2T92nLunr1Klq2bAkiwv/+9z94eHigUaNGBWJ0dHRETk4O+vbti4MHD6JOnTqYPHkyWrVqhUaNGuH+/ftwcnIq0XP18vKCv78/FixYUGDfb7/9hlWrVuHdd9+FtbV1kU1rYWFhmDt3LoKCggAAs2bNQlJSEq5cuYKTJ09i0KBBWLp0KWbMmAEASElJ0Xk2hbGxsYFSqZTPAYDVq1cDAIyNjWFmZobJkyfj66+/xvHjx9GgQYMS3S9jTA+IMVYhkpKSqGXLltSiRQvatGkT+fj4kLW1NVWvXp1Onz5NREQ+Pj4khKCTJ0/qnPvHH3+QEIJ++eUXWr58OQkhSAhBZmZm8vsnX3369CEior1795IQgnx9fek///kPCSGoVq1aFBYWRmq1mhYvXkzz5s0jIQQREWVkZJAQgrp27VroPdy8eZMcHBzI2tqazp49q7NPW9asWbPoypUrciy9evXSOW7UqFGkVCopISGB1q5dS0IIeu+99ygnJ6fUzzYoKIiEEHTmzJkC+0aNGkVNmjQhIqL9+/eTk5NTodcICQmhLl26yD+7ubnR4MGD5Z/d3d3Jzs5O/vn1118nIQQZGRmRsbFxob+HLVu2yMdnZmbK92tsbEynTp0iIiJ/f38SQtBvv/1W6vtnjD0/bgJjrIJ8/vnniI+Px9GjRzF48GD4+vri2rVrqFmzJoYPH478/HwAUq1B27Zt5fNycnKwZMkS2NnZYdiwYfjwww9x7do1hIaGIi0tDY0bN0aDBg2we/duBAYGYv/+/di+fTsAYMeOHfJ16tatCwC4d+8e3N3dUbt2bUyYMAG1a9eWj3ny+KcdOXIETZo0QUpKCvbs2YNWrVrp7H/y3NzcXPl9YGAggoOD5Z/z8vLQoUMH2Nra4saNG3LznrGxMTQaDeLj4xEfHy8/j5Jo2LAhAMi1O3/99RcyMzORlZWFHTt2oFevXgCAv//+Gz169Hjm9WJiYnDx4kU4ODjI20aOHIn4+HicPHkSAHD8+HEcOXIEN27ckGu/xo8fj8DAQAQGBuL8+fMYMGAAAOD3339Ho0aNMGrUKCiVSqxZswZt2rQBAFhZWQEAqlSpUuL7ZYyVHSdAjFWA+Ph4BAQE4Msvv4SlpaW8vWbNmhg3bhyuXbuGbdu2AQASExMxc+ZMxMXF4fTp03jzzTdx4sQJrFmzRv6SdHZ2Rvv27WFsbIxatWqhY8eO6NWrF3r06AFPT08YGUmt29ommqpVq8oJi7e3N0xMTLB8+XKYmZmhY8eOcjza45/+Ml68eDF69+4Ne3t77Nq1Cx06dChwj0+eu3v3blhZWeHChQtwc3ODt7c30tPTAQDbtm2T46tXrx7S09ORmpoKANi6dSvUajXUajVcXFwwYsQIaDSaZz5fc3NzKJVKhIeHY9euXejduzdmz56NkydP4sGDB7C3twcACCFga2srn5ecnIz4+HicO3cOCxYsQExMDACpOQwAPD09AUj9dfr27QtAak7U6tSpE+rWrStff9iwYejRowd69OiB5s2bAwDOnTuHkSNHIjo6Gg4ODggLC8OwYcPkawwYMABbtmxB//79n3mfjDH94QSIsQpw9+5dAMA777xTYN/06dNhZWUFFxcXeducOXNgb2+Pdu3aITQ0FBs2bICXl1epy+/fvz/o32X/pkyZgtDQULkGQqEo+GfgyS/jxYsXY/r06ahTpw5CQ0PRpUuXZ5aVlZWF6tWro1mzZliyZAmuXr2Kfv36ITw8HFlZWRBCAACGDh0Ke3t7fPnllwCkZODAgQMIDAxEnz59sG7dOvz555/PvD9bW1u0atUKhw4dgp+fHwAgNDQU//zzDxQKBd5++20AUiITFBSEuLg4nDt3DrVq1YJarYa7uzsCAwPlGh9tAnTq1Clcv34dnp6ecHZ2BiAls88jLS0NgJR8/f7772jWrFmBYwr7XDDGyhcnQIxVgD179qB27dpQq9WF7lcoFHINCQC4urpi7dq1CAwMRFhYGAYOHFiqcpOSkuTrP+nJ5ivtaKzCjvfz88P06dNRr149BAcHw8bGptiyhBBycqPVqVMnBAQE4PDhw/KXf7du3QAApqamWLx4MVatWoXAwEAoFAp069YNPXr0wIIFC2BmZobo6OgS3auFhQU2b96M3bt3w9jYGEePHsWSJUvQsmVLnaas06dPw97eHu7u7tBoNJgyZQr++usvucbmSbNmzULDhg1hZGSErVu3QgiBQ4cOlSgeLUdHRxgZGYGIcPny5ec6lzFWfngUGGMVwNHRETExMThw4AC6d++us2/VqlVIS0uTawoAwN/fX66hKYujR48WSEiepq0Z0h7/pLlz5wIAli5dijp16iA9PV0nUQOkPiwWFhYFzn3SkCFDsG/fPvj7+0MIoZPQ9ejRA++88w5mzJgBDw8PmJmZAQA+/vhjZGVloWvXriW61z59+uDQoUNo0KABQkJC8Nprr+HKlSuYPn26znE1a9bEjBkz0KxZMwgh4OHhAZVKhfXr1+skW0ZGRvDz80Pbtm3l2rkuXbrgzJkzyMvLk5vxnsXJyUlOsgprOgQezzdkampaomsyxsqOEyDGKsCQIUOwevVqTJo0CZ9++qmclISFhcHf3x/e3t7o2bMnjh8/DrVa/VzJz5MJzJMiIiLw4MEDCCFgbGwMa2vrIq+hrfE5duwYAGmINiB1yL516xaGDh2Kvn37IjAwUO6vozV06FB89dVXOmUV5qeffsKDBw/wxhtvFBgWv2rVKnTp0gWtWrXC4sWLceTIEaxYsQI+Pj5FztvzNBMTE5iYmGDJkiVQq9UYOnQoAgICMHz4cPkYIQQmTZqEadOmFTj/6URRrVZj1KhROtsaNWqEkJAQZGRkPFen5fnz52P+/PmF7tu+fTsGDBiAWrVqYefOnTod4Blj5YcTIMYqgEqlwt69ezFp0iRMnDhR7tjbokUL/Pbbbxg8eLDOsSWVn5+PuLg4vPbaawX23bt3DwBQp04d2NnZwc7ODvXr15drWLRef/11eZLF+/fvQwiBdu3aAZDm0BkxYgTCwsKwadMmAJCTnI8//hheXl5wd3eXl3KoU6dOkc18ZmZm2Lx5c6H7LC0t8ddff2HcuHHo3bs3TExMsHDhQnnSwJJ4//330blzZ7kpS5twaudY0iqqRqxdu3aoWrUqAGkuI20z3ZPeeustLF++vMD2kjbTFaZq1apQKpVQKBRy+Yyx8scJEGMVxMzMDP7+/vD390dOTg7OnTsnJxpaXbp00RlC/ixZWVl48OBBoV/WWk/WyFy9erXA/saNG6Nx48aFntO4cWOcPn26xPFoz3tWs1thbGxssH37doSEhKBWrVpo0qTJc51vZmam049HoVAUqE0xNzeXpwN42sSJE+X3UVFRaN++fYFjevTogTZt2sDc3Fxne1xcHGxsbNCiRYtCr3379m04OjoWuq9bt27P9TtnjOmHoKLqzxljldrx48fRuXNnLFq0CN7e3iU6x9TUFG+99RY2bdr0XEnM02WdOnUKM2fORGBgYGnDf2kcPXoUHh4e+OWXX/Dhhx8aOhzG2L+4BqiUEhIS8M0338De3r7INYgYM6T27dsjLy/vuc7JysrSS1lt27bl5OdfO3bsgEKhQJ06dQwdCmPsCS9FDRARlarKvSyWLFkCb29vGBkZITY2VmdyNcYY03rw4AFu3LiBli1bGjoUxtgTXooESEuj0UChUMijYsorKcrLy0OzZs1w5coVAEDv3r2xa9euCk/CGGOMMVY6lToB+vvvv3H58mXcuHEDAODi4oK+ffsWO1mbPqxatQrjxo2DEEJOtn7++WdMmjSpXMtljDHGmH5U2gTozJkzOiM8zMzMkJmZiZo1a+LLL7+El5dXsaNISttslp6ejgYNGiA+Ph6NGjWS1wUyNTVFWFiYznIGjDHGGHsxVdpO0CqVCj4+PmjRogVUKhXy8/MRExODffv2wcfHB7t27cKyZcvQoEEDnen5g4ODERcXhxEjRpSq3EWLFiE+Ph5t2rRBp06dEBkZCTc3N1y4cAHDhw/HiRMnnmseF8YYY4xVvEpbAwQUXouTnp6OQ4cO4YcffkBKSgoWL14MDw8Peer6QYMGISMjA5s3by4wl8ez3L17F/Xr10dGRgYOHz6M+Ph4DBkyBF5eXrh+/TqioqLw2WefFTnjK2OMMcZeDJV6MVQhBPLy8pCXl4f8/HwA0oyyffr0wYIFC2BlZYWxY8fi+vXr8ro9p0+fhoODg9x35+lhv8UNA/bx8UFGRgb69euHN954Q26CCwsLQ0BAABQKBb777juEhISUx+0yxhhjTE8qdQIESAsWGhkZQalUApBGgmk0GrRq1QqHDh1ClSpVsHjxYuTk5CAvLw937tyBu7s7LCwskJubi++//15OeogI8+fPR2JiYoFyoqOjsXr1aiiVSrmGx9HRETY2NkhKSoKdnR2++uorEBHGjh373POvMMYYY6ziVNoEKCYmBn5+fjhy5Aju378vJzEKhUJe2BEAhg8fjk2bNkGlUuHkyZMQQsjrJoWGhmLu3LnyCsxXr17FrFmzCl3Xx8HBAUeOHMGiRYvkhRyfXDPp1KlT+N///oeRI0diy5YtJV4pmjHGGGMVr9J+S+/evRuTJk2CmZkZWrdujRYtWsDNzQ2NGzeGlZUVHj58iHPnzuG7776TR4MdOnQI9vb28lpAp06dQrVq1XDv3j3Y2toiLCwMVapUgYWFRaFldujQAR06dNDZ5u7ujt27dyMsLAxDhw7Fb7/9Vr43zhhjjLEyq5QJkEajwYQJE/B///d/uHTpErZt24YdO3Zg6dKlMDY2hrm5OdLS0gAAAwYMwJw5cwAAZ8+eRXJyMkxMTAAAR44cgZeXF6pXrw5AqhHKz8/HvXv30KhRoxINldfO7nru3Lnyul3GGGOM6VmlTIAUCgXy8/Oxfv16mJmZ4ddff4VCoUBOTg5Onz6NQ4cOITk5GYMHD0bLli3lYenVqlXDo0ePMGjQILRp0wYHDx7Erl275OaquLg41K5dGzVr1gRQsrmCOAFijDHGKp9KOQz+/PnzeO+99xAREQFTU1N4enpi5cqVsLOzK/a8mzdvYvDgwbC2tsadO3fw4MEDKJVKjBgxAg0bNsSCBQtQv359rF+/HlWqVClRLEQEU1NT5OTkIDMzU+5PxBhjjLEXV6VLgE6fPo0PP/wQubm5WLhwIZKTk/HJJ5/gzTffxNq1a5GTkyN3hH6yM7R2mLxCoYAQAo8ePUJoaCgWLlyI2NhY/PPPPwCALl264ODBg88Vk52dHeLj4xEbGwu1Wq2/m2WMMcZYuah0TWA///wzoqOjsW7dOvTo0QMAkJCQgK+//hoHDhxA9+7dCzRd3bp1CwEBATh//jzu3LkDOzs7dO/eHd26dcO+ffuQlZWF0NBQbNiwQa750U6cWBLVq1dHfHw8kpOTOQFijDHGKoFKNwz++PHjeOedd9CpUyd526BBg6BSqRAcHAyg4CrwN27cwLx587B9+3bY2toiPT0dy5cvR4cOHaBUKtGhQwds3LgRH330ERYvXgwAzzWMXduJOjk5uay3xxhjjLEKUOkSoJs3b0Kj0UClUkGj0YCIkJubi4yMDNSuXVtu6nqSp6cnJkyYIL/fvXs3IiIikJqait9++w1t2rTB8ePHsWnTJjx48OC5Y+IEiDHGGKtcKlUT2Pnz55Gfnw8ikmtoNBoNfHx8YGVlhb59+8ozQj9t0aJFqFq1KhYvXgwzMzNMmDABCoUCw4cPx4ABA5CWlgYhRIk7Pz+pRo0aAICkpKTS3xxjjDHGKkylSoC0tTPh4eEICAhAZGQkbt++je3bt2PIkCHyDM9FmTlzJvLy8jB58mRYWFhgxIgREELAzMwMZmZmpY6La4AYY4yxyqVSJUCtW7fG9OnT8csvv2D8+PHIzs4GAFhYWMDf3x9HjhxB06ZN0bJlS4wcORL16tUDICVOf/75J9atW4fw8HAIIRAREYH8/HwYGRlBo9EAgM6osefBCRBjjDFWuVSqBMjc3BwLFizAggULEBkZifDwcKSlpSE2NhZXrlzBtWvXcOnSJezatQuOjo5yArRu3Tp89dVXcoJib2+PWrVq4erVq7C3ty9Vs9eTOAFijDHGKpdKkwARESIiIpCeno527dqhUaNG8qKkgNQX6MKFC5g3bx7++OMPNGvWTN7n6uqKiRMnIj4+HhqNBsnJyZg7dy6mTZsGS0tLNGrUCJ06dcLXX38NKyur545NmwBxHyDGGGOscqg0CdDFixfx1ltvITs7G6GhoXB2dkZeXh6EEFAqlVAoFIiKisLmzZvxzTffwNjYGAAQFRWFtWvXombNmmjcuDFq1qyJ1q1bw9bWFjdv3sSFCxcQERGBR48elbofkLYTdLnWAIWEAGvWAP7+QCmb6hhjjDEmqTQJ0Pr16xEdHY2dO3fC2dkZQMG5epo2bYrmzZtjzJgx2Lt3L6ysrBAeHg5/f3+oVCrk5eWhbt26cHd3x9ixY1G/fn2MHDmyzMtXlHsTWE4OMHIkEBMDeHpK7xljjDFWapWmKuHOnTuoUaMGHBwcAKDQ+X4aNmyI6dOnIzQ0FNu2bQMgzfuzYcMG+Pn5Yfr06cjLy8OuXbvw5ptvomHDhmjVqhU+/vhj7Nq1C/fv3y9VbOWeAKlUwDffSO+//BLIyCifchhjjLFXRKVJgNq3b4+kpCRs3LgRAKBUKuXRW0969913Ubt2bWiXODMzM8OQIUOwe/duDB06FHv27EFAQAC+/PJL1KlTB1euXMHPP/+Mfv36oX///rh9+/Zzx6ZtAktMTES5La323ntAy5ZAbCzw/fflUwZjjDH2iqg0i6HeunULHTt2xN27d7F582YMGDAAAOR1v7T/vXbtGpo0aYKIiAjUr19fPve1117DiRMn0LZtW/ma+fn52LZtGw4cOICIiAiEhobi4sWLaNq06XPHZ2VlhfT0dKSkpMDa2lo/N/20gwelJjBLS+D6daBWrfIphzHGGHvJVZoaoHr16mHu3LlQqVTw9vaGv78/cnJy5HW/hBBITk7GqFGjYGFhISc/AHDu3DmYmZnBxMQEGo0GOTk50Gg0UCqVGDRoEBYvXow//vgD//zzT6mSH0BaER4A7t69W/abLUq3bkCfPkB6OuDrW37lMMYYYy+5SpMAAcDIkSOxZs0apKWl4YMPPsBXX32FK1euIDw8HIsWLcLrr7+OEydO4NNPPwUAuTnq2LFjcHZ2houLCxQKBVQqFRQKBYgIGo0GFhYWsLOz0xlW/7y0q8DHxcWV/UaLs2ABoFQCK1cCERHlWxZjjDH2kqo0o8AAKaEZOnQo6tevj82bN2Pv3r1YsWIF0tPTodFooFar8eOPP2Ls2LE6512+fBmRkZGYO3cu2rdvj7p168LR0RFmZmYFVo4vrQqpAQIAFxdg/Hjgl1+Azz4Ddu8u3/IYY4yxl1Cl6QNUmBs3biAsLAzGxsawsbGBra0tGjRooHNMTk4OHBwckJWVBbVaDRMTE1hZWcHOzg5OTk5wcnKCp6cn6tevX6ZkaNq0aVi8eDHmz5+Pzz77rKy3Vrx794D69YGHD4EDB6R+QYwxxhgrsUpVA/Q0Z2dneU6goly7dg2ZmZmIi4vD5cuXcfHiRVy6dAnXrl3DgQMHcPfuXWzYsAF79+6FhYVFqWOpsBogALC1Bb74QhoSP2UKEBYGGFXqXyVjjDFWoV7ab03tqLDQ0FCYmJgAAF5//XW8/vrr8jFxcXE4cuQIQkNDy5T8AI8ToHLvA6T1n/9I/YAuXQJWrAAmT66YchljjLGXQKXqBF0aly9fRpcuXWBpaQkiQl5enjx/kFqtxtChQ7F06dIyl6PtBF0hNUAAYGYGLF4svf/qK4DXIWOMMcZK7KVNgLT9ecaNG4eFCxfK242MjKD4dy0tIkJ+fr5eJi+s8BogAOjfH+jeHUhJkZIgxhhjjJVIpe4E/SK5f/8+bGxsUK1atfJdFPVpERFA8+YAEXD2LNCiRcWVzRhjjFVSL20NUEXTzv6cmppa6BId5cbVFfjoI0CjAT75REqEGGOMMVYsToD0xMjICFWqVAERIS0trWIL9/UFatYEjh4F/l0rjTHGGGNF4wRIj5RKJQCU34KoRbG2BubOld5PnQqkplZs+YwxxlglwwmQHuXm5gJ4nAhVqPffB9q3B+Ljgc8/r/jyGWOMsUqEO0HrSVZWFszMzGBkZITs7Gx5pFmFCg8HWrYEcnOBI0eAzp0rPgbGGGOsEuAaID3RDn+3s7MzTPIDAE2aPK79GT8eyM42TByMMcbYC44TID0JDw8HgAJrkVW4L78EGjUCrlwBvv3WsLEwxhhjLyhOgPTk3LlzAICWLVsaNhBTU2mJDEBKgM6cMWw8jDHG2AuIEyA9OXz4MACgbdu2Bo4EwBtvSHMC5eUB770HZGQYOiLGGGPshcKdoPUgMzMT1apVQ3Z2Nu7duwcbGxtDhwRkZgKtW0szRU+aBPz8s6EjYowxxl4YXAOkB8eOHUN2djZatGjxYiQ/gLRY6rp1gLExsGwZsHevoSNijDHGXhicAOnBgQMHAACenp4GjuQpLVoAc+ZI799/H0hMNGw8jDHG2AuCEyA92L9/PwDAy8vLwJEUYto0wMMDSEgAxo3jtcIYY4wxcB+gMktKSoKNjQ2MjY2RkpICc3NzQ4dU0J070orxaWnSCLFx4wwdEWOMMWZQXANURgcPHgQRoVOnTi9m8gMAdetK/YAA4D//Aa5dM2w8jDHGmIFxAlRGx44dAwB07drVwJE8w7BhwLvvSkPihw+XlstgjDHGXlGcAJWRdgLEVq1aGTiSEli2DHBwAE6fBmbPNnQ0jDHGmMFwH6AysrOzQ3x8PG7fvo26desaOpxnO3wY6NoVEAI4ehTo0MHQETHGGGMVjmuAyiA3NxcJCQkQQsDOzs7Q4ZSMhwfw2WeARgN88AEvmMoYY+yVxAlQGaSmpoKIUK1aNRgbGxs6nJLz9QUaNpQWTF2wwNDRMMYYYxWOE6AyyM/PBwAYGRkZOJLnZGoKLF8uvf/mG+DqVcPGwxhjjFUwToDKQKlUAgDy8vIMHEkpdO0KjBolNYFNnMgTJDLGGHulcCfoMsjLy4NKpQIA5OTkVL6aoPv3gcaNgaQkICBAWjmeMcYYewVwDVAZGBkZwdbWFkSEu3fvGjqc51ezJvDdd9L7qVOBlBTDxsMYY4xVEE6AysjFxQUAcPnyZQNHUkqjRwNvvCEtlPrFF4aOhjHGGKsQnACVUbNmzQAAly5dMnAkpSQE8MsvgLExsGIFcPy4oSNijDHGyh0nQGVU6RMgAHB1BT79VHo/fjyQlWXYeBhjjLFyxglQGWkToIsXLxo4kjL673+B+vWBy5eBGTMMHQ1jjDFWrngUWBmlp6fDysoKxsbGePToUeWaEPFpZ84A7dsDeXnAn38CffsaOiLGGGOsXHANUBlZWlrC2dkZubm5iIyMNHQ4ZdO6NTB3rvR+zBggKsqw8TDGGGPlhBMgPWjevDmAl6AZDJCGw/fqJc0N1Lcv8OCBoSNijDHG9I4TID3Q9gOqtEPhn6RQAOvXAy4uQHg48O67wL9LfjDGGGMvC06A9OC1114DANy5c8fAkehJ1arArl1A9erA3r3AtGmGjogxxhjTK06A9ECtVgMAYmNjDRyJHjk7A9u2SfMD/fij9GKMMcZeEpwA6YE2AYqLizNwJHrm4QH4+0vvp0wBtm83bDyMMcaYnnACpAdVqlQBIA2Jf+kMHw7MmSOtFj9sGHDihKEjYowxxsqMEyA9MDExASCtCP9S+vJLYOxYaYbovn2BGzcMHRFjjDFWJpwA6YE2AcrOzjZwJOVECGDZMqBnT+D+feDNN6Vh8owxxlglxQmQHqhUKgAvcQIESJ2hN28G3NyAa9eAfv14zTDGGGOVFidAevBkE9hLvbKIlRWwZw9gbw+EhgIjRwIajaGjYowxxp4bJ0B6oFAoYGRkBAD9okB1AAAgAElEQVTIzc01cDTlrE4daW6gKlWkGqHPPzd0RIwxxthz4wRIT7TNYC9tR+gnNWsGbN0KGBkB330n9Q9ijDHGKhFOgPREuwr8S18DpNW9O+DnJ73/+GPgwAHDxsMYY4w9B06A9EQIAQAvdx+gp40eLQ2R12iAwYOB69cNHRFjjDFWIpwA6YlCIT3KVyoBAoDZs6W5gVJSpJFhDx8aOiLGGGPsmTgB0hNtDZDmVRsVpVAAv/8urR4fEQGMGMEjwxhjjL3wOAHSk1eyCUyrShVg507A2lr6r6+voSNijDHGisUJkJ7k5eUBAJRKpYEjMZAGDYBNm6QaodmzpSHyjDHG2AtK0CtZZaFfRAQjIyNoNBrk5ubKcwK9khYvBqZNA8zNgWPHpJmjGWOMsRcM1wDpQXp6OjQaDSwsLF7t5AcApkyR+gFlZEidou/fN3REjDHGWAGcAOlBQkICAKBGjRoGjuQFIASwciXQpg1w+zbw9ttAZqaho2KMMcZ0cAKkB1evXgUANGjQwMCRvCBMTYHt2wG1Gvj7b2DoUODfPlKMMcbYi4ATID2IjIwEADRs2NDAkbxA6tQB9u0DqlUD/vwTGDcO4O5mjDHGXhCcAOmBtgaoUaNGBo7kBdOkibRwqrk5sGYN8OmnnAQxxhh7IXACpAfaGiBOgArx+uvAtm2AsTGwaBHw1VecBDHGGDM4ToD0gPsAPUPPnsC6dYBSCXzzjbR+GCdBjDHGDIjnASqjrKwsmJmZQalUIisri4fBF2fLFuDdd6UO0dOnAwsWSKPGGGOMsQrGNUBldOvWLQCAo6MjJz/PMnAg8McfgJERsHAhMHYsjw5jjDFmEJwAldHt27cBAPXq1TNsIJXF228DO3YAZmbA6tXAO+9IkyYyxhhjFYgToDJ69OgRAKBq1aoGjqQS6d0bCA4GqlcHdu0CuncHkpMNFs7du9L6rf9W5jHGGHsFcAJURtnZ2QAAlUpl4EgqmfbtpUkSHRyA48eBTp0MloHMmAHMmgUsW2aQ4hljjBkAJ0BlZG5uDuBxTRB7Di4u0oKpTZoA//wDtG0r/VzBPv5Y+q+fH8C/Rsaez31e769Syc7ONsj31Yv4OeEEqIyqVasGAEhKSjJwJJWUvb1UE+TlBSQmAl27Ar//XqEhtGkjTVeUmiqN1meMFY+IcO/ePaxZswa2trYYM2YM9DmgODIyEuPHj0d4eLjerllWsbGx8Pb2RkhISKmvUdL70mg0GDZsGBQKBcaOHVvq8p6UlZWF+Ph49O/fH46OjtiyZUupr7Vz506oVCo4OTkhNja2yOPK+3NSZsTKJCoqigBQ7dq1DR1K5ZabSzRpEpE0QxDRf/9LlJ9fYcWvXy8V26QJkUZTYcUyVuFiY2OpSZMmhe5LS0ujqKgoIiLKyMggtVpN//d//yfvDw4OJnt7exJCyK+uXbuSpaUlPXr0iIiIevbsSSdOnCj0+ufOnStRjJ6eniSEoKVLlxZ5zLx588jKykonFu3LwsKCYmJi6MKFC4XGMnXqVPLx8SEiotWrV5OtrW2h11EqlRQWFkZERB988AEJIWjatGkluofS3hcR0aRJk0gIQVWrViUhBB08eFDed/PmTVIoFLRixYoC540aNYqEEHTkyBEikn6Hvr6+pFAo5Hvq0KEDmZqa0qhRo4qNYcKECTR58uQC2w8fPkympqZybCNHjixwTEk+JydOnKBmzZoV+tyFEPTDDz8UG58+cAJURnl5eWRiYkIA6OHDh4YOp/JbupRIoZCykYEDif79n6W8ZWcT2dlJxQYHV0iRjBlEVFQU1atXr9B9ISEh1KNHDyIiCgsLIyEEOTk5ERHRqlWrSAhBZmZmNGzYMNqzZw8dO3aM4uLi6OLFi/I1unTpQocPHy70+k5OTnT37t1i40tNTaWqVauSQqGg2NjYQo85efKknOgsWrSIfv31V51XaGgoERGNHj2aWrRoUeD8N954g8aMGUMJCQkkhCAjIyPy9fUtcJ2goCD5HEdHR1IoFPK1n1dJ7otISgyEEPTjjz9SYmIiOTk5kaurK2n+/ZfZP//8Q0IIsrGxoStXruicq02Abt++TdnZ2XLCZW9vT3PnzqWgoCC6desWnTlzhpKSkuTzcnNzKSsrS+dakyZNIhMTE53fl0ajIVdXV/Lw8KD09HT66KOPSAihk2SW9HPi4uJCQgh67733Cjz3gIAA+X7LEydAeuDq6koASvyvG/YMgYFEVapI2Yi7O1Exfyz06euvpSL79auQ4hgziP3798tJzdNCQkKoS5cu8s+2trbk5OREDx48IGtra3JycqLw8PAir52dnU3Ozs5yDcTT6tWrR7dv3y42vqioKDkpKUpOTg5ZW1vTzJkzi73W6NGjSQhBiYmJ8rb79++TqakpLVy4kIiImjZtWmgtxtOEEKRQKCgmJuaZxxamJPdFRPT999+TiYkJpaSkEBHR0qVLSQhBe/fuJaLHCZAQgnr16qVz7qhRo0ipVFJCQgKtXbtWTjBycnKKLbN///7k5uams+3ChQukUqno888/l7elpKSQEIKWLVtGRNKzNDExocGDBxORVINYks8JEdHkyZOLTMQrCvcB0gPtEhjXrl0zcCQviZ49pZFhTk7A2bNSJ50zZ8q92PHjAZVKWrw+Kqrci2Ps2TIz9X7Jv//+Gz169HjmccePH0diYiIAIDk5GWlpaVi5ciVcXV3lbfHx8UhPT5fPiYuLQ0pKCtzd3Usdn0ajAQC0a9euyGOMjY1hamqKq1evIj4+Xn6lpKToHNe5c2cAwOHDh+VtDx8+RHZ2Ntq0aQMAsLCwQFRUFOLi4uTrPN1hVxuTWq1GnTp1yu2+AGDPnj2wtLSEtbU1AKBPnz6oVq0a/vzzTwBA3hOTxwYGBiI4OFj+OS8vDx06dICtrS1u3LgBtVqNgIAAGBsbQ6PRyPeXn5+vE9elS5cQExOjE0fz5s0xfPhw+Pn5yfv27NkDAHBwcAAA1KhRA7169cLBgweRmZlZ4s8JAFhaWiItLQ1XrlyR40pISKjQPkKcAOlB/fr1AXACpFeursDJk9Lw+Lg4oHNnYNOmci2yVi1g6FCpE9KSJeVaFHvVbdwodfjv2hXo1k2aC8vDA3BzAxwdAVtbwNISsLPTe9FCCNja2so/a7+gzp07hwULFshfdtopPqpUqYLq1avDysoK0dHR8nkODg5Qq9VQq9Xw9PREWFgYACmh0I6O1Xa8jY2Nxbx58xAdHY3o6GgQEdasWYOLFy8iJiYGMTExuHfvHgDg9OnTAIqfWoSIkJOTg61bt8LZ2RlqtRr29vbo2bOnzginkSNHwsLCAgsWLJC3hYWFQaVSwcXFBQCQk5ODY8eOoVGjRvL9dOrUSadzrzYmY2Pj533cBa5R3H1lZmYiIiICnp6e8raqVauiY8eO8qLbu3fvhpWVFS5cuAA3Nzd4e3vLycW2bdvkFQnq1auH9PR0pKamAgC2bt0q35+LiwtGjBgBjUaD69ev4+bNm+jXr1+BePr27Yvk5GSMGTMGAHD27FkIIdC1a1cA0u93wIABSEpKQmJi4nN9TnJycpCamop27drJ+5s0aYL9+/eX7gGXhkHrn14Sy5cvJwA0ZswYQ4fy8snOJvrgg8edo//3v3LtHH3unFSMhQXRE03kjOnXggWPP9PFvUxMiDIy9Fq0j48PtW3blmJjYyksLIyMjIwKdFYlkprDhBByZ+FZs2ZR9erV5b4i4eHhFBQURL/++iup1Wrq168f3bp1i5RKJR0+fJju3r1LXbt21bm2QqGg27dvk7+/f5EdjlevXk1CCBowYECR93D06FESQlBAQABlZWXR/v376caNG4UeO2zYMBJCyPs9PT3JwcGBiIiio6NJCEGzZ88mjUZDwcHBdPny5QLXOHjwIAkhyN3dvdTPvST3pW0ma968OcXExND8+fNJqVSSEIKaNm1KRNLvT9t0dPToUTI2NqZu3brR5cuXSQhB3bp1IyKizMxMatKkCU2cOJGIiPLz8yk4OJiCgoJoypQpJISg7du3k5+fX4GO1lopKSmkUqmoZcuWRETk4eFBQgjat28fHTp0iJycnOTf36lTp4ioZJ8TIiJnZ2fy9PQkIqLTp08X2WxanjgB0oNdu3YRgALtsUxPNBqiH3543Dn67beJyrHDeY8eUjGzZ5dbEexVd/s20cGD0is4mGjfPqKQEKKwMKKoKKK7d4lSU8tlSKKPj0+BpGTq1KkUGBhIbm5ucgK0ZcsWEkLQrFmziIgoKyuLnJ2d6d133y1wzbFjx1KrVq3o1q1bBRKbZs2aUVBQEE2fPl3uoJuZmUn79++noKAg+RUREUFEjzvy+vv7F3kP2uTsWf2JiB4nL9q+LG5ubvLoJm3CUVSn7aefmTYZLI2S3Jc2Hu3LysqKtm7dSp06dSIhBCUkJOgkQEREGzdulJMkIQTNmTNH3hcUFETGxsb0119/6ZSTm5tL5ubmtGTJEvneCnsG2s7m2iRKmwBpX9OmTaNZs2aREILmz59PRCX7nBBJ/cG0ny1D4dU79UDbVvt0+zPTEyEAb2+gcWNgyBBg+3agY0eps46jo96LmzED2LdPagabNk1atowxvapbV3oZSM2aNTFjxgw0a9YMQgh4eHhApVJh/fr1cvPF0aNHdc4xMTHB0qVL0adPH0yaNAmdOnUCAISGhmLt2rX4WDujKIARI0Zg6NChUCqVaNCgAZycnGBsbIxFixYBAExNTdG9e/dCY9M22axevRrdu3fXWWRaoVDA1tZW7m5w//591K5dG8lPLKWjUqlQvXp1+WcPDw9069YNfn5+mDx5MtLT09G7d28A0LlOfn6+3OcJAJRKJWxsbHRi2rRpE9577z1YWlrKxwkhUKtWrWc88ZLdl9bIkSMxevRoeHh4QAiBjIwMhIaG4lghE8UOGTIE+/btg7+/P4QQGDhwoLyvR48eeOeddzBjxgx4eHjA7N8/Zh9//DGysrLQtWvXYucD2rRpE+zt7TF79mx5m6OjI3766Se0bt0atWrVQlZWFhYuXIjjx48DKNnnJDU1FYmJifL8eYmJiTr9kmxtbaFQVEAPHYOmXy+JkydPEgBq3bq1oUN5+V25QtSggVRFY2ND9Pffei9CoyFq3Voq4uef9X55xgzK19e3yJqMUaNGyTVAbdq0ISEEffvttzrHfPvtt2RpaUkzZ86kPXv2kK2tLbVo0UKeQ6hu3bqFXruktTZffPFFkXPDCCEoKCiIvL29SQhB5ubmVKdOHZ397du3L3DN0NBQEkKQpaUlCSHozJkzRCSNuNKOzHqyOUcIQfXq1aPc3FwiIlqxYkWxMa1cubLYeyrJfe3bt0+uAQoJCdE5V1sT88MPP5Cvr2+B0VMZGRk0cOBAWrJkSYFyHz58SO7u7tS4cWPau3cvff755ySEIF9fXyIi+vrrrwst8+LFi6RSqejHH3+Ut3l4eBQ6f5Crq6vcTKZV3Ofk3Llz8n3Xr19f5zmYmprqDJcvT1wDpAfaXvlPZvSsnDRqJHWOHjIE2L9f6kS6fDnw/vt6K0IIqRZo0CBg4UJpdBj/atnLRAhR6PZ27drJCztrOyU/PWrpiy++gIWFBWbMmIHs7Gz07t0bv/76K6pUqYLk5OQi/+Xu4uKC5s2bw8rKqtjYvvnmG9y+fRsbNmyASqXCuHHjoFar4e7ujlu3bsHV1RXXrl2DjY2NXJMjhICnpydsbGwKHWXVoUMHdOvWDQcPHoSDg4M8Ss3Z2RmWlpZ4++23oVQq5WPr1q0LNzc3+W/6+PHjcePGDXz33XdQKBQYOXIk6tati/bt2yM6OloebVaW+3JxcUFERARMTU3Rvn17nXPbtm0r1xBRIaOkzMzMsHnz5kLLtbS0xF9//YVx48ahd+/eMDExwcKFCzF16lQAQK9eveDj44ONGzeiS5cuAICMjAx8+OGHaNOmDT755BMAQG5uLqKiovB+IX9r+/bti6CgIJ1txX1OatWqhSpVqqBr166oVq0aOnfujPr166N169ZwcHCQO6iXuwpJs15ye/fuJQDUvXt3Q4fy6sjNJfL2ftxZdMoUaZue5OU9rmjasEFvl2XM4BYsWEAbN2585nGOjo46swo/LTw8nIKDg3UmrIuOjpY7ub5ozp49S9bW1jT7Be7ct2zZMqpevXqh+3x8fGjbtm00a9asUs2fU1Qn7/z8fPLy8iJTU1NasWIF+fn5kZeXF7322msUHR0tH5eUlETm5ua0c+fOAteOjIykgQMHFlpuYZ+TFwUnQHqwevVqAkAjRowwdCivHj8/ImNjKVPp2ZPo38nD9GHFCumyzZtX6KocjL0QGjVqRJ07d6YMPY9CY2Vz8uRJ6tmzp16veeXKFapVqxYJIcjExIS8vb0pPj5er2W8iATRi7QyWeX07bff4r///S8+++wzzJ8/39DhvHqOHgXeeQe4f19qItu1C/h3csqyyMoC6tcHYmOBbduAt9/WQ6yMMcZeCDwRoh7Ex8cDAGrXrm3gSF5RnTsDp08DzZoBkZFA27ZS/6AyMjUFvvhCej9rFvDvRK6MMcZeApwA6cHdu3cBAHblMGsrK6F69YBjx4B+/YDUVODNN4GlS6UeQmXwwQdAnTrAhQvAjh36CZUxxpjhcQKkB9rREiWZC4KVI0tLqa3qv/8F8vOBTz6RhnBlZZX6klwLxBhjLydOgPRAu3CedtIsZkAKBTBnDrB+vZS9rFoFtG8PXL9e6ktqa4EuXuRaIMYYe1lwAqQH2gSoRo0aBo6Eyd59V1pRvn594Px5oFUroJgZT4vDtUCMscouOztbZ6HYiqL9fnwRcQJURhqNRp7OmxOgF0yLFsCZM8DAgcDDh9LMhh9/DPy7yvXzeLIWaPv2coiVsQoUGxsLb29vhISEGDoUHZGRkRg/fjzCw8MNHYpMH8+qpPel0WgwbNgwKBQKjB07ttTlPSkrKwvx8fHo378/HB0di1364ll27twJlUoFJycnxMbGFnkcEeHevXtYs2YNbG1tMWbMmEIncDQ4w47Cr/wePHhAAMjc3NzQobCiaDRES5Y8ni+odWuiq1ef+zI//yyd3qwZzwvEKrcPPvhAXszyWe7du0e9e/fWWTz1hx9+0DmmY8eO9OmnnxZ6/qNHj2jSpEl07ty5Z5bl6elJQghaunRpkcfMmzePrKysCl1OwsLCgmJiYujChQt04sSJAudOnTpVXgZk9erVZGtrW+h1tCvTEz3fsyrLfRERTZo0iYQQVLVq1QIrtN+8eZMUCgWtWLGiwHnahVa1k1ZmZGSQr68vKRQK+Z46dOhApqamhS5l8aQJEybIi8U+6fDhw2RqairHNnLkyALHBAcHk729vc6z7Nq1K1laWtKjR4+IiOjEiRPUrFmzIpcEefqzVZiAgABSKBSFLrj6PDgBKqPk5GQCQNbW1oYOhT3LyZNEjo5SFmNhIU2i+Byzk2ZlEdWpI53+xx/lFyZj5c3R0ZEUCgWFhoY+89i+ffuSEIIGDRpE3333HdWtW5cUCoXObNK+vr7k6OhY6PnLli0jlUpFR48eLbac1NRUqlq1KikUCoqNjS30GO2aWBYWFrRo0SL69ddfdV7a+xk9ejS1aNGiwPlvvPEGjRkzhhISEuQ1wHx9fQtcJygoSD7neZ5Vae+LSEoMhBD0448/UmJiIjk5OZGrq6s8g/I///xDQgiysbGhK1eu6JyrTYBu375N2dnZcsJlb29Pc+fOpaCgILp16xadOXOGkpKS5PNyc3MpKytL51qTJk0iExMTunv3rrxNo9GQq6sreXh4UHp6On300UckhNBJMletWkVCCDIzM6Nhw4bRnj176NixYxQXF6eztpeLiwsJIei9994r8NwDAgJKNGO0tiwhBC1cuPCZxxeFE6AyunfvHgGgGjVqGDoUVhLJyURDhz5eQqNfP6J790p8+i+/SKc1bqzXlTcYq1DampyYmJhijwsICCAhBPXr148ePHhARNLfvCZNmpCJiQn9888/REQUFBSks8jok0aNGkVNmjR5ZkzahUCNjIyKPCYnJ4esra1p5syZxV5r9OjRJISgxMREedv9+/fJ1NRU/sJs2rRpobUYTyvpsypKSe6LSFqY1cTEhFL+nc1+6dKlJISgvXv3EtHjBEgIQb169dI5d9SoUaRUKikhIYHWrl0rJxg5OTnFltm/f39yc3PT2XbhwgVSqVT0+eefy9tSUlJICEHLli0jIulZmpiY0ODBg4mIKC0tjaytrcnJyYnCw8OLLXPy5MmlWsrjSVlZWaRWq0kIUWTiXRLcB6iM8vPzAfBCqJVGtWrAhg3AunVAlSrAzp3SBIq7dpXo9PffB157DbhyBfj993KOlbFyoPm3F79arUadOnWKPI6IMGfOHLzzzjvYsWOHvIipjY0NQkJCUKNGDaxcuRIA0LBhQwCQO9n+9ddfyMzMRFZWFnbs2IFevXqVOK7CFjPVMjY2hqmpKa5evYr4+Hj5lZKSonOcdnHSw4cPy9sePnyI7OxstGnTBgBgYWGBqKgoxMXFydd5usNuSZ9VWe8LAPbs2QNLS0tYW1sDAPr06YNq1arhzz//BPB40W0ACAwMRHBwsPxzXl4eOnToAFtbW9y4cQNqtRoBAQEwNjaGRqOR70/7faWN69KlS4iJidGJo3nz5hg+fDj8/PzkfXv27AEAODg4AJD6u/bq1QsHDx5EZmYmkpOTkZaWhpUrV8LV1RUAkJycjPj4eKSnp+tc39LSEmlpabhy5YocV0JCwnP1ETIxMUH//v0BADk5OSU+72mcAJWR9pf2PL889gIYNkzq0fzGG0BCAvDWW9IK8wkJxZ6mUkkjwQDA17dU/akZw65d0kDFd98FBg8GevaUEuuaNYEuXYCZM4EvvwQ++kj/ZZ8+fRqAlEwUZ8OGDYiLi8OqVasK7LOxscGHH36Ibdu2AQDMzc2hVCoRHh6OXbt2oXfv3pg9ezZOnjyJBw8ewN7eHoD0d3LNmjW4ePEiYmJiEBMTI8+jpo1LpVIVGRMRIScnB1u3boWzszPUajXs7e3Rs2dPnRFOI0eOhIWFBRYsWCBvCwsLg0qlklcaz8nJwbFjx9CoUSOo1Wqo1Wp06tRJp3NvSZ9VcUpyX5mZmYiIiICnp6e8rWrVqujYsSMiIyMBALt374aVlRUuXLgANzc3eHt7y8nFtm3b5H+E16tXD+np6UhNTQUAbN26Vb4/FxcXjBgxAhqNBtevX8fNmzfRr1+/AvH07dsXycnJGDNmDADg7NmzEEKga9euAKSO1QMGDEBSUhISExNRvXp1WFlZITo6Wr6Gg4ODXK6npyfCwsIASM89NTUV7dq1k/c3adIE+59z9v6hQ4cCkEaZ3bhx47nOlZWpHopRRkYGASCVSvVCrnbLniEvj+j774nMzaW2rWrViFavLrZvUF4eUZMm0uFLllRgrOylsWDB41bY4l4KBVF2tn7LPnjwIAkhyN3dvdjjPD096auvvipy//fff0/GxsaUmppKRERt27alwYMHy32G3njjDfrll19IqVTSnTt3iIjI39+/yA7Hq1evJiEEDRgwoMgyjx49SkIICggIoKysLNq/fz/duHGj0GOHDRtGQgh5v6enJzk4OBCRtGq9EIJmz55d5Crpz/OsilOS+9I2kzVv3pxiYmJo/vz5pFQqSQhBTZs2JSJpNXht09HRo0fJ2NiYunXrRpcvXyYhBHXr1o2IiDIzM6lJkyY0ceJEIpJWew8ODqagoCCaMmUKCSFo+/bt5OfnV6CjtVZKSgqpVCpq2bIlERF5eHiQEIL27dtHhw4dIicnJ/n3d+rUKSIimjVrFlWvXl3uUxQeHk5BQUH066+/klqtpn79+hERkbOzM3l6ehIR0enTp+WO289L+8wUCgUdPny4VNfgBEgPzMzMCAA9fPjQ0KGw0oqKklaT137zdOtG9G//hsJs3y4dZmtLlJ5ecWGyl0NEBNH69dJrwwai3bulbXFxRJs3E/33v0S+vkRr1xJlZuq3bB8fHxJCyKOhCpOUlERCCLk/ytO0/YDGjBkjb+vatav8pahSqUgIQS4uLtS6dWv5mMzMTNq/fz8FBQXJr4iICCJ63JHX39+/yLhCQkLkzr7Pok1etH1Z3Nzc5NFN2i/PZ31xluRZPUtJ7ksbj/ZlZWVFW7dupU6dOpEQghISEnQSICKijRs3ykmSEILmzJkj7wsKCiJjY2P666+/dMrJzc0lc3NzWrJkiXxvhT0DbWdzbRKlTYC0r2nTptGsWbNICEHz588nIqlfjrOzc6Ejs8aOHUutWrUiIqJ69erRrFmzSv4An/HMOAEyMHt7ewJAt27dMnQorCw0GqLffyeqUUPKboyMiKZMISrkS0CjIWrTRjrM19cAsTJWSt7e3iSEoMaNG9O1a9fo7t278is+Pp6ISO5IW5iDBw+Sq6srtW/fXucffYsWLSIhBDVs2JBiY2PJxMSEhBBFDo9/Wr9+/UgIQZ07d6bo6GiduBISEoiIaOXKlSSEoLNnz1J2drbOMU+ObiKSaj48PT2pRo0aFB0dTc7OznKH4n379pEQgrZu3Up5eXk617n3xKCIkjwrfdyX9st81KhRFBISIrcmaDuhb9++vUACRET0/vvvy0nA0yPDhgwZQs2bN6eMjAx524QJE0ihUNClS5eKTYCmTp1KDg4OdP/+fSKSEqB69erR7t275fvOzMwkKysr6t+/v3ze3r17SaFQ6Iz4+/vvv0mlUtG0adMoJSWFLCws6JNPPiEiKZF+8nnkP8f8IpwAvSDc3NwIAJ09e9bQoTB9SEwkGjeOSAgpw7GxIVq5Umr7esKhQ9JuU1OpAomxymDFihVFzsEihKCVK1fKQ7L79u1Lv/zyC/n5+dHKlSupS5cu8pB47agwrZ9++olMTU0pMDCQiKSaD4VCQefPny9RXF988UWxcQUFBckJiTigfsgAACAASURBVLm5OdWpU0dnf/v27QtcMzQ0lIQQZGlpqTNK7fvvv5dHZj3ZnCOEoHr16lHuv0M8S/Ksynpf+/btk7/MQ0JCdM7V1sT88MMP5OvrWyABysjIoIEDB9KSQtriHz58SO7u7tS4cWPau3cvff755ySEIN9//8X29ddfF1rmxYsXSaVS0Y8//ihv8/DwKHT+IFdXV7mZTOvbb78lS0tLmjlzJu3Zs4dsbW2pRYsWlJaWRufOnZPvu379+jrPwdTUVGe4/LM8WWvGCZABdevWjQDQvn37DB0K06ezZ4k6dXrcLNayJdG/7d1a2hH177xjoBgZK4XPPvtM7n8zZswY8vHxocDAQPLz85NrErZs2ULt2rWTv2S8vLxo3rx58oR2T8vIyKALFy7IP+fn59PJkydLHJNGo5H77ZiYmNBHH31E3377LQUFBdGKFSsoOjqafvrpJ7K1taUxY8bQmDFj6P3336d169bRvn37KC0trdDraufEqVu3rrztzz//JCsrKxo5cqR8LT8/PwoKCipQs1OSZ1XW+woKCiIzM7MCc/IQEdWqVYt++OGHQmuAnuXevXtyDZSpqSktWrRI3nfmzBkSQtCHH34ob3v06BG1b9+eOnbsKG/LycmhunXr0m+//Vbg+jNmzCh0vqUff/yRTE1NSQhBffr0kWu64uLiqGrVqtS/f3/5uX/zzTc6TaElFRMTQwqFgiwtLen69evPda4WJ0B6MGjQIAJAGzZsMHQoTN80GqmThr3940Ro61Z5d3S0NKciQLRrlwHjZIwV6uzZs2RtbU2zZ882dChFWrZsGVWvXr3QfT4+PrRt2zaaNWtWqebPKaqTd35+Pnl5eZGpqSmtWLGC/Pz8yMvLi1577TWKjo6Wj0tKSiJzc3PauXNngWtHRkbSwIEDCy03PDycgoODX+jBQYKIx2+X1cSJE7F8+XL89NNPmDx5sqHDYeXh0SNg4ULpJQQQGQnY2QEAvv8emDpV+jE8XJpqiDHG9OnUqVOYOXMmAgMD9XbNyMhIeHh44N69e1CpVJgwYQK++OIL1KpVS29lvMh4HiA9qF69OgDIi6Kyl5CFBeDjA9y4AYweDXz+ubzrk0+Ajh2Bu3cBb2/DhcgYe3n9f3v3Hd9U1f8B/HOTpklHuinQXdpSVqHsJRt8GKKIDEUehoqCIgj8HicgCiooUxFkK7gZAipYNiKrgFIQoezSFuiA7qYjyff3x/He7lJoIG35vl+v87q3yU1ycinNJ+ece06bNm0sGn4AIDQ0FDdv3oTZbEZOTg4WLlz40IQfgAOQRcirwHMAegh4egKffirC0MmTAAC1GlizBrCzA9atA9avt3IdGWOM3REHIAvgAPQQqlcPCA8H/p1aPiRE9I4BwAsvAJcvW7FujDHG7ogDkAVwAHqIqdXK7rhxwJNPAunpwNNPA5VYooYxxth9xgHIAjgAMUCMjV61CvD3B44dE2s5McYYq5o4AFkAByAmc3UFvv8esLEB5s0D/l1EWTAYgGKrVjPG2IOQm5tbZMHYByU5OfmBv2ZFcQCyAGdnZwBAenq6lWvCqoJ27YAPPhD7I0YAygLJdnbAxYvAa68BV65YrX6MxcfHY+LEidi7d6+1q1JEdHQ0XnzxRZw5c8baVVFY4lxV9H2ZzWYMGzYMKpUKL7zwwj2/XmE5OTm4efMmBgwYAH9/f2zYsOGen2vLli2wtbVFYGAg4uPjyzyOiJCYmIgvv/wSnp6eGD16NKrkjDvWnYaoZsjMzCQAZGdnZ+2qsCrCZCLq00dMkNiuXbEVvRcuFMt8Dx1K9O/U/Iw9SM8//7yyqOWdJCYmUr9+/ZQZoVUqFS1cuLDIMR07dixzza+srCx6+eWX6a+//rrja8mzNn/22WdlHjN79mzS6/WlLivh4OBAcXFxFBUVRUeOHCnx2MmTJysLm65evZo8PT1LfR55hXqiuztXlXlfREQvv/wySZJEzs7OJVZqv3z5MqlUKlq2bFmJx8kLrsorq2dnZ9OMGTNIpVIp76lDhw6k0+lKXdKisLFjxyqLxha2f/9+0ul0St1GjBhR4pjdu3eTj49PkXPZrVs3cnR0VGYQP3LkCIWFhZW5NEjx363ybN++ndzc3Ein05W6qv2dcACyAJPJRJIkEQAyFlsvij28kpOJfH1FCJo0qdidQ4cWzCzdrRvRpk1E/64/xNj95u/vTyqVig4ePHjHY/v376+s//XJJ5+Qn58fqVQq+v7775VjZsyYQf7+/qU+fsmSJWRra1tkgczSpKamkrOzM6lUKoqPjy/1GHltLAcHB5o3bx6tWrWqSJHfz6hRo0pdoqFz5840evRoSkhIUNYCmzFjRonniYiIUB5zN+fqXt8XESnrry1atIiSkpIoMDCQGjVqpMykfPbsWZIkiWrVqlViCQ45AMXExFBubq4SuHx8fOijjz6iiIgIunr1Kh0/frzIorH5+fkllt94+eWXSavV0o0bN5TbzGYzNWrUiLp06UKZmZk0fvx4kiSpSMhcuXIlSZJEdnZ2NGzYMPr111/p0KFDdP369SJrfDVs2JAkSaLhw4eXOO/r1q2r0MzROTk5yvIi8nb06NF3fFxxHIAsRKVSEQBlET3GiIgOHxaLyhdbQUM0EX30EZGjY0EQ8vMTtyUlWa2+7OEgt+TExcWVe5y8GvkTTzyhLH6amJhIjRs3Jq1WS2fPniUiooiIiCKLjRY2cuRIaty48R3rJC9uaWNjU+YxeXl55OLiQtOnTy/3uUaNGkWSJFFSof9LycnJpNPpaO7cuURE1KRJk1JbMYqr6LkqS0XeF5FYoFWr1VJKSgoREX322WckSZKygr0cgCRJor59+xZ57MiRI0mtVlNCQgKtXbtWCRh5eXnlvuaAAQOoWbNmRW6LiooiW1tbevPNN5XbUlJSSJIkWrJkCRGJc6nVamnIkCFERJSWlkYuLi4UGBhIZ86cKfc1X3nllXta0qOwmTNnkkajUdbfbN269T0FIB4DZAEmkwlmsxmSJEFd6LJoxtq1Az75ROyPHi2GAAEAVCoxm3RsLLBwIRAcDFy7Brz1FuDjI2abPn7cWtVmNZjZbAYAeHl5wdvbu8zjiAizZs3CwIEDsXnzZuj1egBArVq1sHfvXri7u2P58uUAgPr16wOAMsh2+/btMBgMyMnJwebNm9G3b98K16tt27ZlHqPRaKDT6XD+/HncvHlTKSnFLi7o1KkTAGD//v3KbRkZGcjNzUXr1q0BAA4ODrhy5QquX7+uPE/xAbsVPVeVfV8A8Ouvv8LR0REuLi4AgMceewyurq7YunUrAMBoNCrH/vbbb9i9e7fys9FoRIcOHeDp6YlLly7By8sL69atg0ajgdlsVt6f6d95y+R6nT59GnFxcUXq0bRpUzz77LNYsWKFct+v/17N4evrC0Bc+NO3b1/s2bMHBoMBt2/fRlpaGpYvX45GjRoBAG7fvo2bN28iMzOzyPM7OjoiLS0N586dU+qVkJBwV2OERowYgZUrV6JXr15ISUlBVFQUmjZtWuHHKyoVwxgRiT5uAKTVaq1dFVYFmc1itXiAKDycyGAo5SCTiWj7dqK+fYkkqaBVKDyc6LPPiG7ffuD1ZvfP7su76bXfXqPXfnuNJmyfQGN/GUujN4+mYRuHUb9v+tGQ9UNo5YmVlJVX+srrlSF3tQQGBpZ73DfffEN6vV5pkSjuvffeU7q9EhISyMbGhpYsWUJbt24lSZLorbfeon379indOkSiK2XNmjUUFRVFsbGxFBsbq6wU/v333ytjRspiNpvJzc2NNBoN2dvbK+N1WrduTZmZmcpx+fn55OjoSG3atFFu27hxI2m1WkpMTCQioubNm5NarSZHR0elZSU0NLRIS09Fz1V5KvK+srOzycvLS2lRISK6ffs29e/fX3ncRx99RE5OTnTq1Clq3rw5NW7cmDIyMoiIyM7OTjluzZo15OzsrPy7/fjjj8r7CwkJoeHDh5PJZKLo6GiSJImee+65EvXZtGkTSZJEPXv2JCKiSZMmkUqlUs6xwWCgr7/+Wul2S0tLIycnJ1q9erXyHPK/j16vp+7du9OJEyeU55IkiZycnJR6ubu7F+l2vBuDBg2iZs2aKWOM7gYHIAuIiYkhAFSnTh1rV4VVUampREFBItO8+OIdDr5wQQwacnMrCEJaLdGwYUR79oiwxKq1j//4mDADdyxuc9zoSsoVi772nj17SJIkatmyZbnH9ejRg6ZNm1bm/QsWLCCNRkOpqalERNSmTRsaMmSIMmaoc+fOtHTpUlKr1XTt2jUiEh/OZQ04Xr16NUmSRE899VSZr3ngwAGSJInWrVtHOTk5tHPnTrp06VKpx8pjQ+T7e/ToQb6+vkREFBsbS5Ik0cyZM8tcLf1uzlV5KvK+5G6ypk2bUlxcHM2ZM4fUajVJkkRNmjQhIrEqvNx1dODAAdJoNNS9e3f6+++/SZIk6t69OxGJcNK4cWMaN24cEYkxqrt376aIiAglfPz000+0YsWKEgOtZSkpKWRra0vNmzcnIqIuXbqQJEm0Y8cO2rdvHwUGBir/fpGRkUQkArGbm5sypujMmTMUERFBq1atIi8vL3riiSeIiCgoKIh69OhBRETHjh1TBm7fC/m8ffvtt/f0eA5AFrBv3z4CQB06dLB2VVgV9uefIscAROvWVeABBgPRd98R9exZEIQAonr1iD74gOgexyQw6zsWf4zmH5pP8w/NpwWHF9DnkZ/TyhMrae3JtbTl3BZafnw5tV7empoubVqhQaF349133yVJkpSroUpz69YtkiSpzNYfeRxQ4XEX3bp1Uz4UbW1tSZIkatiwIbVq1Uo5xmAw0M6dOykiIkIp//zzDxEVDORds2ZNmfXau3ev0upwJ3J4kceyNGvWTLm6Sf7g3L9/f7nPUZFzdScVeV9yfeSi1+tp48aN9Mgjj5AkSZSQkFAkABGJliU5JEmSRLNmzVLui4iIII1GQ9u3by/yOvn5+WRvb0+ffvqp8t5KOwfyYHM5RMkBSC5Tpkyh9957jyRJojlz5hCRGJgcFBREzzzzTInne+GFF6hFixZERBQQEEDvvfdexU9gGUwmEz322GMUEBBQZGD33eAAZAGrV68mAPTss89auyqsilu2TGQYe3uiO4wVLOrKFaLp04l8fAqCkEpF9OijRGvXEv3bFM5qDrPZTLey7+0Pe3kmTpxIkiRRgwYN6MKFC3Tjxg2l3Lx5k4hIGUhbmj179lCjRo2offv2ShcMEdG8efNIkiSqX78+xcfHk1arJUmSyrw8vrgnnniCJEmiTp06UWxsbJF6yd1ky5cvJ0mS6MSJE5Sbm1vkmOIfgiaTiXr06EHu7u4UGxtLQUFByoDiHTt2kCRJtHHjRjIajUWeR+4iq+i5ssT7kgPQyJEjae/evUrolQeh//TTTyUCEBHRc889pwzSLn5l2NChQ6lp06aUnZ2t3DZ27FhSqVR0+vTpcgPQ5MmTydfXl5KTk4lIBKCAgAD65ZdflPdtMBhIr9fTgAEDlMdt27aNVCpVkSv+/vjjD7K1taUpU6ZQSkoKOTg40IQJE4hIBOnC58N0F63bx48fV1ql7hUHIAuYOnUqASi3uZgxIjEeaPhwkV8aNryH3GI0irFCgwYRaTQFYcjenuiZZ4h+/ZXoDld+sIfbsmXLypyDRZIkWr58uTL2pX///rR06VJasWIFLV++nLp27apcEi9fFSZbvHgx6XQ6+u2334hItHyoVCo6efJkher11ltvlVuviIgIJZDY29uTt7d3kfvbt29f4jkPHjxIkiQp43zkq9QWLFigXJlVuDtHkiQKCAhQruatyLmq7PvasWOHEoD27t1b5LFyS8zChQtpxowZJQJQdnY2DRo0iD799NMSr5uRkUEtW7akBg0a0LZt2+jNN98kSZJoxowZRET0/vvvl/qap06dIltbW2XcFpEIQKXNH9SoUSOlm0z24YcfkqOjI02fPp1+/fVX8vT0pPDwcEpLS6O//vpLed/BwcFFzoNOpytyuXx50tPTlcvpp0yZQitWrKCVK1cqXa0VxQHIAp599lkCUG4TJ2OyjAwRfgCiZ58Voeie3LpF9MUXRI88UrSLrFYtoldfJTp6tBJPzmqy119/XRl/M3r0aHr33Xfpt99+oxUrVigtCRs2bKC2bdsqH1C9evWi2bNnlznYNDs7m6KiopSfTSYTHT16tMJ1MpvNyrgdrVZL48ePpw8//JAiIiJo2bJlFBsbS4sXLyZPT08aPXo0jR49mp577jn65ptvaMeOHZSWllbq88pz4vj5+Sm3bd26lfR6PY0YMUJ5rhUrVlBERESJlp2KnKvKvq+IiAiys7MrMScPEVHt2rVp4cKFpbYA3UliYqLSAqXT6WjevHnKfXILyksvvaTclpWVRe3bt6eOHTsqt+Xl5ZGfnx999dVXJZ7/jTfeKHW+pUWLFpFOpyNJkuixxx5TWrquX79Ozs7ONGDAAOW8f/DBB0W6Qivi/PnzpNVqi0z0KEkSTSox4Vr5OABZQPv27QkA7du3z9pVYdXEmTOi0QYQ3WKVdvky0axZRA0aFA1DQUFEb7xBFBnJYYg9lE6cOEEuLi40c+ZMa1elTEuWLCE3N7dS73v33Xdp06ZN9N57793T/DllDfI2mUzUq1cv0ul0tGzZMlqxYgX16tWL6tWrR7Gxscpxt27dInt7e9qyZUuJ546OjqZBgwaV+rpnzpyh3bt3W3wMmyVJRFVxgY7qpU6dOkhISMC1a9eUeRIYu5Ovvwb++19AqwUOHAD+nZ6kcoiAP/8EvvkG+PZbICGh4D5/f+Cpp0Rp107MRcQYqxYiIyMxffp0/PbbbxZ7zujoaHTp0gWJiYmwtbXF2LFj8dZbb6F27doWe42qjANQJWVnZ8PBwQEajQYGg4EnQmR35aWXgOXLgdq1gSNHgIAACz65yQT88QewYQOwaRNw/XrBfV5ewIABQL9+QNeugL29BV+YMcaqPg5AlXTu3Dk0bNgQQUFBuKhM88tYxeTlAX36AHv2AA0bAgcPAq6u9+GFzGaRsDZuFIHo2rWC+3Q6oFs3oG9fUerVuw8VYIyxqoUDUCXt2bMHPXr0QJcuXbBv3z5rV4dVQ6mpwCOPAGfOiMaY334T3WL3DZFYZuOXX4Bt20ouuREaCvTuDfToAXTuDDg738fKMMaYdfAggEq6/m+3gpeXl5VrwqorFxeRQ+rWBfbtA154QWSU+0aSxICj994Djh0Dbt4EvvwSGDJEhJ3oaGDRIuDxxwE3N6BtW7FG2c6dQHb2fawYY4w9OByAKikxMREA4OnpaeWasOrMzw/49VfAwUEMjv7f/+5zCCqsdm1g5Ejghx+ApCRg/35g2jSgY0cxUDoyEpg9G3j0UZHWunQR4enAAdGHxxhj1RAHoEqSV+i1tbW1ck1Ydde8uRieY2MDzJsHvP++FSqh0Yhur/ffFwOoU1JEn9zrrwOtWgFGI/D778CMGeI4FxegVy9g1ixxe06OFSrNGGN3z8baFajuzGYzAEDFlxQzC+jdW1y9/vTTImPo9cDkyVaskKMj8J//iAKIQPT772LU9p49wN9/A7t2iQKIwUvt2olWos6dgfbt+QozxliVxAHIQuQgxFhlDR4shtqMGgVMmSLyw9ix1q7Vv1xdgSeeEAUQ8wwdOCC6zfbvB06fLtgHRItS69YiDHXpIrrV9Hrr1Z8xxv7FAaiS5LE/8lggxixh5EggMxMYPx4YNw7IzQUmTrR2rUpRuzYwaJAoAHDrlug6k0PQyZPAoUOizJ4txhS1aCHCUJcu4vK3+3LdP2OMlY8vg6+k7du3o2/fvujVqxd27Nhh7eqwGuazz4AJE8T+Bx8Ab79t3frctbQ0MbmRHIiOHxcTNBYWGiquNGvTRmybNgV4TB1j7D7jAFRJUVFRCA8PR5MmTXD69GlrV4fVQKtWAWPGiKvCJk8GPv4YqLYTjmdmAocPizD0++/A0aMlryTTasWI8DZtRPdZs2ZAgwaiO40xxiyEA1AlJSUlwdPTE25ubrh165a1q8NqqO++A0aMEBdhDRggLpV3cLB2rSwgLw84dUoEochIsY2OLnmcRgM0aiTCULNmopWoWTOgVq0HX2fGWI3AAaiSzGYzdDod8vPzYTAYoNPprF0lVkPt3QsMHChmjm7ZEti6VSzpVeOkpoqusqNHxcKuUVHApUulH+vpKVqHQkPFVi7+/tW4mYwx9iBwALIAPz8/xMbG4vLlywgMDLR2dVgNdu6cWL/08mXx2f/dd0D37tau1QOQmSmuMIuKEuXUKVEyM0s/XqsFQkIKAlFQkFhpNiAA8PERky0xxh5qHIAsoF27djh69Cj++OMPdOzY0drVYTVcUhIwdKhoEVKpxJyFb70l9h8qZjMQHy9SoVyio8U2Pr7sx6nVIgTJgcjfXzSl1alTUGrXFovEMsZqLA5AFvD444/j559/xqZNm/Dkk09auzrsIWAyAe++K64MA8Q8hWvWiPXEGICMDBGG5EB05Qpw9aoo169XbJ0RZ+eCMFSrFuDhUbAtvC9v7ezu97tijFkQtwNbQK1/B2ImJydbuSbsYaFWi9UnOnYEhg8HIiKAJk2ApUvFmqYPPb1eLN3RqlXJ+3JzgdjYgkB09apYEDYhQWzlkpYmSmmDsktjb19+QPL0FIu++fkB7u5iUVrGmNVwALIADw8PAOKKMMYepD59xFCY558XIWjoUGDLFjF/kJubtWtXRWm1QHCwKGUxm8WyH3IoSk4WJSmp6FbeT0oS03fHxIhyJ3Z2BWHI379gX+6W8/bmcUqM3Wf8P8wCeB0wZk3e3sD27cAXXwD/939iLbFdu4CFC8WaYjWxoSE6ORqhHqH37wVUKtFK4+4uLr+/EyLR7VZeULpxQ7Q8XbtW0LJUVuuSWg34+opwJIeiwoUHcjNWafw/yAIMBgMAwI7HADArkSSxZEbPnsBzz4nVKIYNA776SnSL1aSLE386+xOe+vEpTG4/GTO7zYSdpgr8v5MkwMlJlHr17nx8WpoIQoWL3Hokj1OSu+fkddUKKz6Qu3Aw8vQUxcODQxJj5eD/HRbAAYhVFSEh4vNy9Wrgf/8T3WKNGwPvvCNmka4Jv6IXb1+EJEmYd3gefj7/M1Y/vhod/arZ1ZfOzkBYmCilyckpOU7p6tWiAUkOTKUFJJmbmwhDtWoVhCJn54KwVrwUvs/OrmY2HzL2L74KzAJGjBiBdevW4csvv8TIkSOtXR3GAIjhK5MmibmCANGbMncu8NRT1f9z7Vj8MYzaMgr/JP0DCRJebv0yPuzxIZy0Ttau2oNR2kDumBggLk50uyUmioVpzeZ7fw21umJBqbzi7CymLK/uv3CsRuIAZAGDBw/Ghg0b8MMPP2AIX4LDqpg9e0QQOnVK/Ny5M7BggViUvTrLNeZi1oFZmP3HbBjNRnjpvbCo9yI81fApSPyBK+ZKuH1bhCE5FCUnA+npBSUtrejPhW/PzbVMPQp3D95tkNLpRNFqRZH3dbqKde8RifNgNAL5+ZXblnWbjY2ok61tQT0LF51OtKYVLjrdQzhxV9XDAcgC+vXrh23btuHnn3/GY489Zu3qMFaCyQSsXAlMnSo+AwFg8GAxiWKDBtatW2WdSjiFF39+EUfjjwIA/hP0H3za51PUd69v5ZpVc3l5YmB3WSGpvJKWJh6bng5kZd2f+qlURUNE8Y8yosq1gN1vWq2YOqF4OCpcit9vYyNK48bikk9WKRyALKB79+7Yu3cvdu3ahR49eli7OoyVKTVVTJ64eLEYZqJSASNHikkV/f2tXbt7ZzKbsOLPFXhr91tIzUmFRqXBpPaTMLXTVOi1emtX7+FmNIolS8prcSotQKWni1aonJyS25yciocblUqEBo2mYtu7PTY/X4TF3NyCUvjnnBzAYBAlO1tsK9u6NnAgsHFj5Z6DcQCyhPbt2+PIkSM4ePAgOnToYO3qMHZH8fHAzJnAqlXi80mjERMqvvGGWFe0ukrMSsSbu97EmpNrAAC1HWpjVvdZGB0+GmoVL45aoxiNJVt9ind9Fm8lqirM5qLBqHCRQ1JpRe52a9CAW4AsgAOQBYSHhyMqKgp//vknmjdvbu3qMFZhFy8CM2aIgdJms/j8GDhQrC3WsqW1a3fvIuMjMfG3iTgSdwQAEOYZhjk956B3cG8eH8QYAwBUwWhc/WRnZwMA7O3trVwTxu5OcDDw9ddiPr4XXxQtQRs3ihUkunYV+0ajtWt599p4t8Gh5w7h24Hfwt/ZH6cTT6Pvt33R7atuOBR7yNrVY4xVAdwCZAEBAQGIiYnB5cuXEViTZpxjD53r18UVYl98IYZtAGJuvbFjgTFjxFQy1U2OMQdLji3BBwc+wG3DbQBA35C+eK/re2jlVcpaYYyxhwIHIAvw9fVFXFwcYmJi4OfnZ+3qMFZpaWnA2rVisPT58+I2jQZ44glg1Cix+nx1m2Q4NScVcw/NxcIjC5GVL65Meqz+Y5jWeRraeLexcu0YYw8aByAL8Pb2xvXr1xEXFwdvb29rV4cxizGbgd27RRD6+eeCMad16gD//a8IQxVZKqsqScpKwseHPsaSY0uQnS+6r3vV64W3HnkLXQO68hghxh4SHIAsgFuA2MMgLg5Ytw748suCViEAaNYMGDJElPIWWK9qErMSseDIAiyOXIzMPNHf19qrNV7v+DqebPAkXzXGWA3HAcgCGjZsiHPnzuHMmTNoVN2+DjN2l4iAI0eANWuAH38U3WWyFi1EEBo4UKxLVh2kGFLwWeRn+PTop7hluAUACHQJxMS2EzG6+eiHZ3kNxh4yHIAsoE2bNjh27BiOHDmCtm3bWrs6jD0wubnAjh0iCG3ZIib/lTVoAPTvL0r79lV/zFB2fjbW/LUGC44swKWUSwAAva0eo8JH4ZXWryDUoxpPkMQYK4EDRT8iYgAAIABJREFUkAX07NkTu3fvxo4dO9CrVy9rV4cxq8jJAX77DVi/Hti2Tcw6LXNzA/r0AXr1Anr0EFeWVVUmswk/n/8ZC48sxP6YgpXWe9bribEtx+Lx0MehUWusWEPGmCVwALKAJ598Eps3b8bGjRsxcOBAa1eHMavLzwcOHgS2bhWDpy9eLHp/aCjQs6cojzwCeHhYp553cirhFD6L/Azfnv5WGTBdx7EORoWPwnPhzyHEvZr08zHGSuAAZAEjRozAunXr8OWXX2LkyJHWrg5jVQoRcO4cEBEhrijbt69gjiFZSIjoJmvXTmybNKlaXWYphhSsO7UOXxz/AmeTzyq3P+L3CEY2G4nBjQbDWedsxRoyxu4WByALGDduHL744gssWbIE48aNs3Z1GKvS8vOByEgRhnbvBo4dE8scFWZvL0JQ06aiNGsGhIUBrq7WqbOMiHAo9hBW/LkC6/9Zr7QKadVa9A/tj2eaPIM+wX1gp7GzbkUZY3fEAcgCpkyZgvnz52Pu3LmYMmWKtavDWLWSnw+cOiWuLDt8WGwvXSr9WC8v0X1Wv74o8n5AgJio8UHKyM3AxrMb8VXUV9h/dT8I4k+po60j+oX0w1MNn0Lv4N68Gj1jVRQHIAuYNm0aZs2ahffeew/Tp0+3dnUYq/Zu3wZOnxbB6NQpICoK+Pvvki1FMhsboF69glAUHCxKUBDg63v/u9Pi0uPw/d/f44czP+D49ePK7bZqW3QL6IbH6j+GPsF9EOQWdH8rwhirMA5AFvDhhx/inXfewZtvvomPPvrI2tVhrEYymYBr18TCrefPF5ToaHF7WWxsRAtRUFBBKJJLvXqAnYV7q66kXMGms5vw07mfcCj2kNIyBABBrkHoFdQLPQN7omtAV7jbu1v2xRljFcYByAJmzZqFadOm4Z133sGsWbOsXR3GHjoGg7jSTA5EFy+KbrRLl4D4+PIf6+1dNBTJwahuXbHkR2W61hKzErH9wnb8euFX7Lq8Cyk5KUXuD/MMQ5eALujk1wmd/Dqhrr7uvb8YY+yuVKHrLKqv/Px8AIDmQQ9CYIwBEK04YWGiFGcwAJcvFwQiuVy8CFy9KgJSfDzw+++lP3etWiIMyYHI1RVwcSkorq6ihcnHR/ysUhU81tPBEyPDR2Jk+EgYzUacuH4COy/vxJ4re3Ao9hBOJ57G6cTTWBy5GICYgbqDbwelNPFsAhsV/5lm7H7g/1kWkJeXB4ADEGNVkZ0d0LixKMUZjUBsbMlgFBMD3LgBJCQASUminDp159fSaIDu3YEBA4AnnhChSWajskFbn7Zo69MWUztPRY4xB8fij2F/zH4cuHYAh2MP40rqFVxJvYJvTn8DAHDQOKCNdxu082mHdj7t0Na7LWo71rbQmWHs4cZdYBYwduxYLFu2DJ9//jlefvlla1eHMWYhJhOQmCjC0I0bwM2bYu2z1NSyS0YGkJ4u5j9q0QLo1w/o3Rto0wZQl7O+qslswunE0zgUewiHYg/hcNxhXE65XOK4AJcAtPNphzZebdDWpy2a12nOl91XApH4dzYaC7ZmswizcincqsdqDg5AFjB48GBs2LABP/zwA4YMGWLt6jDGrIxIdL2lpxcEIoMB0GoBJyfRrebmdufnSchMwNH4ozgcdxhH444iMj4SWflZRY6xUdmgae2maO3VWhTv1mhUq9FddZ1lZAC3bol6ZmcDeXniPQCAJAG2toBOV1C02oKtjY0ICGq1OPZ+IhIBJSdHlPR0EUjlrVwK/5yeLkpmpnifhbcGgwg9d6JWiyAknwd7e9GyKJfiP1fmNvlnDl33HwcgC+jatSv279+PXbt2oUePHtauDmOshjKZTTiTdAZH4o4gMj4SR+OP4kzimSJXmgGAzkaHZrWboUXdFmhepznC64SjsWdj2Gvs7+l1b94UXYVxcQVjpgrvx8cDWf/mMrVaFDkUyduybpMkEWwKF6Bg32wWi+7KocdsrswZLJ1aLYKcjU1BnYxGEQT/HeL5wGk0BecGKHpuevcGfvnFOvWqSTgAWUBgYCCuXr2K6Oho1K9f39rVYYw9RDJyM/DnjT9x7PoxRMZH4sSNE6V2nakkFYLdghHmGYbGno3R0KMhGno0RLBbMBxsHSpdj9RUEYri4grCUmysmKLg2jWxn5NT6ZeBSiVaSOTWNGdnUQrvF//ZyQlwdAT0elHkfTu7O7dcya1O+fkFQcxgKCjZ2UV/ttRt5endG9i+vfLn8mHHAaiS8vPzodPpQETIycmBra2ttavEGHvIpRhS8NfNv3Di+gmcTDiJkzdPIjo5GiYqvb/HW++NYLdgBLkFIdAlEIEugfB38Ye/sz/q6uta7Eo0s1mECbnk54tiMoliNoswIgcSSRKBp3DXW1VaI+5+IRKtT7LC50Mu5Y0nYxXDAaiSLl26hODgYPj5+SEmJsba1WGMsVLlGnNxNvksziSewZmkMzibfBZnk87icspl5JvL7udRSSp46b3g4+SjFG+9tyhOBVudje4BvhvGKu8hyNL31+XLoqm5Xr16Vq4JY4yVTWujRXidcITXCS9yu9FsRExqDC6lXMKl25dwNe0qrqRcQUxaDK6lXcPNzJuIS49DXHpcuc/vbucOHycf+Dr7iq2TL/yc/ZStt5M3bNXcQs6qDg5AlXT+/HkAQHBwsJVrwhhjd89GZYMgtyCxTlkpS5XlmfIQnx6P+Ix4xKbFIj4jHnHpcYjPiFduv55xHbcMt3DLcAtRCVGlvo4ECV56L/g5+xUpvk6+8HX2ha+TLzzsPSDd70vJGPsXB6BKOnfuHACgQYMGVq4JY4xZnq3aFoGugQh0DSzzGDOZkZiViNi0WMSlxyE2PVaUNLGNSY3BjcwbIjRlxONw3OEyX8tb7w0vvRe89F6o41gHtR1qw9PBE7UcasHdzh3u9u5w1bnC1c6Vu91YpXAAqqTo6GgAQGhoqJVrwhhj1qGSVKjjWAd1HOugtXfrUo/JN+UjPiMe19KuISY1BrHpsbiWdk0JSnHpcUjJSVFmw64IW7Ut9LZ66LV66G31sNfYw8HWAXY2drDT2ClbnY2uoKgL9uVj7DX2sNOIrVwcNA5wtHWEg60DL0dSQ/Eg6Ery8/NDbGwszp8/j5CQEGtXhzHGqq2svCxcz7iulISsBNzMvImk7CQkZSWJbrbsW0jJSUGKIaXcwduWZKu2FWFI4wAHWwc4aByUoFQ8aGnVWmhttLBV20Kj0kCj0sBGZQMblQ3UKjXUkhoqSaUUSZIgQSqxBVBmd6Cfsx86+3d+IO+9JuMAVAlZWVlwdHSERqNBdnY2bB6G6zMZY6yKyDHmIC0nDVn5WcjIzUB2fjay8rNgyDcgOz8bOcYcGIwG5Bpzi2zl23OMOTDkG2AwGmDINyArP0s8R14WsvKzlK2Z7sPsi5XwVMOnsGHIBmtXo9rjT+xKiIsTV0X4+flx+GGMsQdMZ6ODzvH+jgMiIuSacpGVl4XMvEwlGBmMBmWbnZ9dJGTlmnKRZ8pDvjkfRrMR+SaxNZEJJrMJBFK2ZjKDiEAgZSu/bllae5XezcjuDn9qV0JiYiIAoHZtXp2ZMcZqIkmSlDFD7vbu1q4OsyBebq0S5ADk6elp5Zowxhhj7G5wAKqEzMxMAIBer7dyTRhjjDF2NzgAVYI87sdoNFq5Jowxxhi7GxyAKkGj0QAQC6IyxhhjrPrgAFQJ9vb2AIC0tDQr14Qxxhhjd4MDUCXIC6DKC6IyxhhjrHrgiRArIScnB/b29lCpVDAYDEqXGGOMMcaqNm4BqgSdTgcfHx+YTCbExMRYuzqMMcYYqyAOQJUUHBwMALh48aKVa8IYY4yxiuIAVEkcgBhjjLHqhwNQJXEAYowxxqofDkCVxAGIMcYYq344AFWSHIAuXbpk5ZowxhhjrKL4MvhKSk1NhaurK+zt7ZGZmQlJkqxdJcYYY4zdAbcAVZKzszMcHR2RnZ2NlJQUa1eHMcYYu+9iYmKq/TqYHIAqSZIk+Pj4AADi4+OtXBvGGGPs/kpPT0e3bt3QuXPnar0SAgcgC9DpdAB4UVTGGGM1X0xMDPLy8nD48GGEh4fjq6++QnUcTcMByALkcT/V8ReAMcYYuxthYWE4deoUBg8ejIyMDIwaNQpDhw6tdsNAbKxdgTu5evVqkZaVwoOMiw84vpf7LHGc3A+alJSEGzdu3NfXuh/PIQc3IiqyDwAmkwk5OTnIzc1Fbm6usl/abfn5+cjLy0N+fn6J/fz8fKhUKtjY2ECj0RTZarVa6HQ66HQ62NnZldgv7TYbmzv/6hIR8vLyitT5TvtGoxEqlQqSJEGlUiml8M82Njb3VCRJKnGOK7Nviee4077ZbIbJZCpRjEZjqbdb4n75PrkOkiSVWuR/l/KK/F7MZnO5+/d6PxFBo9HA1tZWKVqttsjPZd0n75e1lfflNQbL+vcvfl/hc1b897es22oCIoLRaITBYEBOTk6pW5PJpLzfsn6n1Go11Go1bGxslP3SfpaL/DtrNBqRn59f5r78O202m4sU+TZA9CZotVrlb2LxfXt7e+h0OqhU1m+7cHNzww8//IB+/fph/PjxWL9+PQ4fPoy1a9eiW7du1q5ehVT5q8CCg4P5EnNWglqtVsKQVqtV/sjIf3TkLWOsfPKH+/0o8pecwgVAkQBZVjGZTEW+oJRWDAZDkYAjB4maTv4yaG9vD3t7e2W/+Lb4bVqttsx/l3v5t5RLXFwcJkyYgBMnTkCSJLz66quYNm0aPDw8rH2qylXlA1CPHj1w7do1AEW7mIpX+17us9RxycnJMBqNcHNzK7Ii/P2sk6XeV+FvQ/K28L5KpVJCRnnfTgp/y9VoNMq3Ynlfo9Eo39CKfzPKy8sr81tbWd/kKvqHztbWVmk9KlzvwtvC+2q1usQ3/8Lf1Aq3UtxNKd6KWdb5vtt9Sz9fad+Oi3/jLeubsCXuL3yfSqW644dk8daY4qV4S1Fp+5W5HxAtwHl5eeWW3NzcEvsV3ebl5d3x37/4fYXPT/FWrOK/2zVJ4S9GhbeFW47v9PtUXktlaS2X8u9t8Zbt4vvy73RZLctEVG4re+G/gdWBo6MjMjIyrF2NclX5AFQdtGnTBseOHcORI0fQtm1ba1fnoZCfn6/8QcjNzVX+2BQucrcTY6x08od+ad02liyFW2aBsrs1iwdv+QtKaUWr1ZYIORXpGq/uzGYzcnJykJ2dDYPBUO62+H5eXt49/btV5LisrCxkZWUpobpWrVpITEy08tkqX83/bXkAeBD0gyeHHL1eb+2qMFZtFW7h02q11q4OqwCVSqV0b1UF+/fvx//93//h+PHjAIAmTZpg7ty5+M9//mPlmt2Z9UdS1QAcgBhjjD1Mzp49i8cffxxdu3bF8ePH4eXlhVWrVuHkyZPVIvwA3AJkERyAGGOMPSxiYmIQFhYGk8kEBwcHvPHGG5g8eTIcHBysXbW7wgHIAjgAMcYYq6mKT7Hg7++PoUOHQq/XY8aMGahTp44Va3fvOABZAAcgxhhj1V1ycjKSk5MRFBSkXLlb+OpC+ZiTJ08iJCQEXbp0QZ06dYpcUVydcABijDHGHlJyeDlw4AB69+6Njh07YsOGDdBoNMjKysLZs2dx5MgRnDp1CocOHUJ0dDRUKhVCQ0PRqVMnACUn160uOABZALcAMcYYq8quXbsGPz8/fPnll/juu+8we/ZsNG/eHEajERqNBp6enmjQoIEy59CKFSvwv//9D/b29ggKCsLp06fh4eGBn376CV5eXggKCoKLi4u131alcACyAA5AjDHGqorY2FgcOHAAe/bswdGjR3HmzBl07NgRBw4cAADs3LkTeXl5WLBgAcLDwwEAderUQf369fHHH39Ar9ejf//+aNWqFdzd3dGkSRMMGTIEN27cQP/+/ZXXqa5dXzK+DN4COAAxxhh70NasWYN33nkHgPj8+eWXX+Dj4wN/f3+8/PLL+PPPP9G6dWvMnz8fixcvBgB07NgRAHDgwAEMHToUt2/fBgA4OTkhKCgIkiRBp9Ohfv366NKlC5o0aQJALEuVnZ2trMxQeF216opbgCyAAxBjjLEHKS4uDmPGjIGvry8GDRqEkJAQvPTSS0hISMDs2bPRvXt3eHp6wsXFBY6OjlCpVDCbzQgJCYGLiwtSU1Nx8eJFPP/881i2bBk8PT0RHR2NuLg4REVFITw8HCaTCYBYYiQkJARbtmzB6dOn4efnZ+V3bxncAmQBHIAYY4xZmtw6k5ubi88++wy//PKLstTEhg0bYDabce3aNcyfPx+Ojo7w8/NDgwYNMH78eLRq1Qp+fn5wcnJS1q1TqVRISUmBvb095syZg6+++go7d+7Eo48+irNnzyIoKAgAEB0drdRB/nyrX78+dDodTp069SBPwX3FAcgCOAAxxhirjNu3byMiIgKvvfYaWrZsCXt7e7Rs2RLvvvsuiAjbtm3D448/jrlz54KI8Pnnn2P+/PmYMGECDh06hPj4eIwYMQK3bt3C999/j82bN2PcuHFo1KgR5syZg7y8PABAeno6nJ2dERsbi+HDh2PNmjXIzc1F8+bNMW/ePNStWxdOTk4AoCzYCgD16tWDp6cn/vnnH+W+6o67wCyAAxBjjLGKKj54+KOPPlLG8vj7+6NDhw7o378/tm7dipkzZ6JFixYICQlBREQE5s2bh5iYGOTn56NTp06oX78+du3ahaioKLRq1Qomkwljx46Fv78/9Ho9hg0bhpdeegkajQYA4O7ujpYtWyoDogcPHoywsDC89NJLOHDgADQajTIwunAd69atC2dnZ1y8eBF5eXmwtbV9UKfrvuEAZAEcgBhjjJXl1q1b+OWXX7B9+3ZERkaiT58+eOONN5SxNFlZWVCpVOjcuTMWLVoEb29vuLm5oVOnTnjsscdw48YN9O7dG9988w2Sk5OxdOlStGzZEi1btsT169fxyCOP4OjRoxg2bJhy/0svvYTk5GR4eHgUqYu9vT2aN2+O9evXK7c1aNAAP/74I9544w3Ur18fdevWLfIYs9kMlUqF5ORk2NvbIz09vcTzVkccgCyAAxBjjLHiTCYT1Go1pk+fjqVLlyIsLAwajQZLly7F8ePHMXPmTDz66KMIDw+H2WxG3759lTW2AODUqVMwGo3o2rUrXF1d0aVLF5w+fRrjx4+HWq0GALi5uaF27dr4559/EBAQoNwGoEhIiY+Ph6urK+zt7REaGorc3FxcuHABISEhMBqNqF27NpYvX15qy4782fb++++jVq1aNSL8ADwGyCI4ADHG2MNlw4YNWLRoEZKSkgCU/PtvNpuhVquxa9cuLF++HJ07d8bGjRtx5MgRfPPNN7hw4QJeeOEFnD9/Hs2aNUNISAgSEhIAiKuutm/fjqlTp2LgwIFo0KABXF1dERwcDJVKhYkTJ2LcuHEAAJ1Oh4SEBMTExCAvLw8eHh44ePAgzp8/j6VLl+Lxxx+Hn58ffH19sXfvXgCAn58f3N3dsX//fgDiM4yIYGtrW+rnmFqtBhGhQ4cOCAkJKfV8GAwGy5zYB4gDkAVwAGKMsZpPngPn/PnzmDlzJhYtWoQLFy6UeqxKpYLJZMLvv/8OAJgyZQqCg4Ph6uqKZ555BlFRUfDw8MDrr7+uXJo+b948jB49Gj179sSQIUPg5OSEhQsXAgBsbW0RFBSEtLQ0nD9/Hmq1GkajEQAQEhKClJQUREREoHnz5vj000/RsWNHTJw4EX///TeGDh2Kw4cPo2/fvgBEC1FwcDD27NkDoOiYpLLm9il8e2pqKrZt24b//e9/aNmyJby9vfHtt99W9vQ+cByALIADEGOM1SxXr17F2rVrMXr0aDRq1AgqlQrPPfccACAzMxOnT5+Gg4MDGjVqBKDs4LBx40aMGTMG/fv3L/IZ4evri9dffx1bt25FTEwMwsLCoFar8dVXX0GlUmHGjBnYv38/6tatq3SJBQYGws7ODn/++ScAKLeHhoaicePGsLGxgY+PD2xtbfHKK6/g888/x4YNG/D++++jbdu2ymt7eXnhhRdeUC6zl7vTCktOTsaCBQuQlZUFAEhLS8OkSZNQt25duLm5YeDAgfjll18QGBiISZMm4fnnn6/U+bYGHgNkARyAGGOs+srJyUF2djbc3Nzw3Xff4bXXXkNSUhKcnJwQHByM9u3b45VXXkGvXr0AiEvCAeDy5cvl/t2/du0arl69iu+++w5AQUuLPKhYrVbDxcUFmZmZaNGiBVavXo1+/fph7dq1cHR0VK7ckj9jAgMDERISgt27d+Ppp59Wbm/fvj28vb0RGhqKdu3a4ZNPPlHGARVWOKSNHDlSCXSlhbdPPvkEn3/+OSRJwmuvvYZvv/0WS5Ysgb+/Pz7++GMEBQWhdu3acHd3h6Oj412f86qAA5AFjBs3Dn379kXDhg2tXRXGGGN3ITU1FRMnToS/vz8mT56MiRMnIjk5GZ9//jnatGlTZDZlOby4uLjA29sb8fHxyM7Ohqura5FuJHl/9+7daNeuHRwcHJTXK7yExMaNG9GhQwc0btwYcXFx0Gg0aNGihfJ8MnnOnZCQEEycOBGbNm0CAGXAspubmxJ47OzslMfJkyaWNmePjY0NkpOTQURISUnBoUOHMGrUKABirqDVq1cjOzsb8+fPx8iRI3Hr1i3k5+dj9uzZGDhwoEXOvbVxF5gFPPHEE5gwYYIyAp8xxljVcO7cOaxbtw5ZWVm4du0aIiMji9zv4uKC69ev4+LFi9BqtahVqxY6d+6MESNGoGXLlvD19YVery/R0l+/fn0AwMmTJwEUhA2goGvq5s2bytgdIlJafVQqFWbOnIk//vgDEyZMACDW2goNDVWev6yWpX79+mHVqlUVeu8qlapE+JGf9/Tp0/D09MScOXOwb98+TJ8+HYMHD8bt27exceNGBAQEYPr06UhPT8fmzZsxfvx4eHt7IzIyEmlpadi+fTsmT56MV199VZk5urr1gnAAYowxViPcunULmzdvxquvvormzZtDp9OhUaNGmDp1KgwGA3788Uc8++yzuHTpEgAog4hv3boFo9EIOzs7DBgwANeuXcOmTZvw008/4cUXX0Tjxo2xaNEiAGK8TE5ODkJDQwEAV65cASACkBwAbGxE50pubi5iYmIAANnZ2UhMTMTBgwfRp08fzJw5E1OmTMGjjz4KQHSr1a5dG8nJyTCZTOXOtFw4bJWHiJQik4PcpUuX4Obmhr///hsZGRmws7PDxo0bsXjxYmzevBlhYWGYMGECGjVqhN27d8PFxQVt2rTBxx9/DD8/Pzz11FPYuXMnWrVqBS8vryLPXV1wFxhjjLFqS+5u2rJlC4YPH46srCw0adIEYWFhePbZZxESEoLmzZvDw8MDLi4uuHTpEtavX48333xTGfwbFBSEzMxMAEDLli0xf/58vPTSS/Dy8oKLiwtGjRqFUaNGKa9lMplgNpuh1WqV7i15vI58lZarqysaNGiAc+fOYcCAAfD19cXp06eh1WoREBCApKQkODs7AyiYaHDMmDFwdnZW5g8qS0WXoSgcSG7evIm//voL+/btQ3R0NB599FE888wz2Lx5MwYOHIjw8HBcuHAB77//PsxmMzZs2AA3NzeMGTMG06dPByAuda9fvz7OnTuH69evQ6vVwt3d/S7/xaoODkCMMcaqPDl8zJ07FwsXLsTmzZvRqlUrSJKEzMxMrFy5EllZWZgzZw6effZZ6HQ6ODk5Ka0xANClSxd06dIFq1evRseOHdGpUydl7hy5pcbe3h65ubn48ccfMWjQoFJnU9bpdFCpVMjLy1MCw59//omDBw8iOjoaubm5mD17NsaNG4cVK1Zg/fr12L59O0JDQzF8+HD06NGjSNCRA82QIUMscq5yc3MRGRmJo0eP4vDhw4iOjsa5c+eg1WpRt25d9OnTB507d8aGDRuQlpaG/Px89O7dG+vXr8fkyZNx6NAh5ZL57t274/nnn1danQICApCTk6O0+siq4/IYHIAYY4xVeXJrhkqlgpeXF44ePQoAaNasGRwdHdGuXTucOHEC3bt3Vz6ciQiJiYnw9PQEESEkJARz5sxBu3btMHXqVOzfvx8GgwHHjx9XXicwMBAAUKtWLQBFZ1OOi4uDu7s77Ozs0LJlSzg4OCAmJgZjxowBALRr1w4LFy6En58fQkJCoNfr8fzzz2PkyJFFgpisvFaee2U2mzF16lTMmzcPDRo0QK9evbB161Y0btwYP/zwA1xdXeHm5gaNRgN/f3/s27cPHh4eCA0NhSRJePrpp/H2229Dp9MBAG7cuAEAuHjxIkJDQ/H777/j5s2b0Ol02LJlCw4ePIht27Zh8ODB+OSTT6rVFWEcgBhjjFVZubm5uHjxIrKysrBjxw6kpqYiMjIShw8fxocffghHR0e8/fbbeOaZZ/D+++/jo48+QnBwMPbt24d//vkH7du3x/bt25VWljZt2mDMmDFYsWIFvv76awwfPhyOjo6oV68eDAYDPDw8oNfr8ccff6BOnTrYtWsXtm3bhpMnT+LGjRvYuXMnevTogf/+97+wtbXF/PnzYTabMWHCBPTs2RM+Pj4lgo2NjQ3MZrMytqfwKuuWJg+wfuuttyBJElxdXXH58mXY2NjA39+/yFViDRo0ACCuHPP394enpycOHDiAFi1aKN1yHh4eqFu3Lnbs2IHQ0FAsWrQIvXr1QkJCAjIzMxEWFoaxY8fimWeeqVbhB+BB0IwxxqowtVqNBQsWoF27dpgzZ44SLtq3b48ffvgBly9fRr9+/aDVaqHX6/H777/jwoULaNSoEdauXYuIiIgSY2amTZuGZs2aYdKkSbhx4wb8/f3h6uoKOzs73Lx5E4GBgZg2bRrat2+PyZMn4+LFixg1ahQiIyPRo0cPmEwm2Nra4r///S/++usvREVF4fnnn4e/v79Sv+JXRKlUKmg0GqjV6vs+WFin08HNzQ2urq5BzYWyAAAD40lEQVQAxBVmSUlJuHz5MgDRXQVAaSk7fPgwnJycEBYWhoMHDxapv5ubG5o3b46///5bWQbD29sbo0ePxtKlS/HNN99g5syZaNy48X19T/cDtwAxxhirsuTZjfV6PTZs2ICePXsq44HkhUVbtGiBlStXAgAmTpyId955p8yBxGazGT4+Pvjoo48wduxYdOzYEfHx8crcNk5OTvD398eVK1cwefJk1KlTB61atULDhg2h1WoBlOy6kq8mKxxuqsIVUUajETY2NmjatCl2796NS5cuFQkqLi4uAIC//voLarUajzzyCFavXg2gYKC1s7MzZsyYAbVajdDQUGzZsgUdOnSo1oOfZRyAGGOMVWmtWrWCnZ0dXF1dS0w2WLduXYSFheHy5csIDAwssbxDYmIiTp48CQcHB3Ts2FHp2unduzfmzp2Lp59+ushl5XXq1MGSJUvg6OioXKV1J6WN76kK5BDTtGlTmM1mnD9/HkDBuZHH+cjrmY0ePRqrVq1CRkYG9Hq9cmyrVq2U5+zfv7+yX95Ei9VB9aw1Y4yxh4Y8QeDixYuV7huZs7MzGjZsiH/++QcNGzbE9u3b8e2332Ly5MkIDg5G3bp18eSTT+LQoUMAioaVJ598EnPmzAEAdOjQAYC4nN3b2xvOzs4gIpjN5grPu1PVyMEkNDQUzs7OuHjxIoCCAKTX62FnZ4fExESkpaXB19cXZ8+eVcJPaQrPd1TaRIvVSdWMrYwxxti/AgIC0L9/f/z+++9ISkqCt7c3zGYz1Go1srOzcenSJbRv3x7p6ek4d+4cXnjhBbi4uMDLywsrVqxA27ZtSx2jolKpMGnSJAwaNAj+/v4l7pckqUp0ZVUGEcHR0RF2dna4evUqUlNTla6vunXrws/PDxcvXsT169fh7OwMe3v7cp+vOgee4jgAMcYYq9I0Gg2GDRuG1atXY8eOHRg9ejTUajVMJhM2b96MlJQUvP766/jpp59w9OhRjB8/Hk2aNEFAQAD8/PwgSVKRtbpkcsApLfzUFHJQVKlUcHNzK9IC5unpibCwMHTo0AHe3t5WrKV1SFTdFu9gjDH20ElKSkJAQADatGmDuXPnIiMjA5s2bcLatWuxbNkyDB06FPn5+cqMzEyQg19SUpIytxETOAAxxhirFlq3bo2zZ8/C3d0dsbGx0Gg0WL58OUaOHFnkOHnF9ZrQhcXun5rTmccYY6xGa9WqFbKzs+Ht7Y3o6GhkZ2cr4afwd3m5y4fDDysPjwFijDFWLTz55JM4d+4c3n77bWVSPrmLh8MOu1vcBcYYY4yxhw53gTHGGGPsocMBiDHGGGMPnf8HjCsser0d4bQAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "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=1, n\n", " print *,i, \"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/tmp32w1xbxr/src.linux-x86_64-3.6/hellofortranmodule.c\u001b[0m\n", "\u001b[39mcreating /tmp/tmp32w1xbxr/src.linux-x86_64-3.6\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/tmp32w1xbxr/src.linux-x86_64-3.6/hellofortranmodule.c\"\n", "\u001b[39m adding '/tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.c' to sources.\u001b[0m\n", "\u001b[39m adding '/tmp/tmp32w1xbxr/src.linux-x86_64-3.6' to include_dirs.\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmp32w1xbxr/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmp32w1xbxr/src.linux-x86_64-3.6\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[39mget_default_fcompiler: matching types: '['gnu95', 'intel', 'lahey', 'pg', 'absoft', 'nag', 'vast', 'compaq', 'intele', 'intelem', 'gnu', 'g95', 'pathf95', 'nagfor']'\u001b[0m\n", "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n", "\u001b[39mFound executable /home2/nhn2/miniconda3/envs/pangeo/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: gcc -pthread -B /home2/nhn2/miniconda3/envs/pangeo/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC\n", "\u001b[0m\n", "\u001b[39mcreating /tmp/tmp32w1xbxr/tmp\u001b[0m\n", "\u001b[39mcreating /tmp/tmp32w1xbxr/tmp/tmp32w1xbxr\u001b[0m\n", "\u001b[39mcreating /tmp/tmp32w1xbxr/tmp/tmp32w1xbxr/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmp32w1xbxr/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgcc: /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/hellofortranmodule.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/hellofortranmodule.c:15:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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[39mgcc: /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.c:2:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.c: In function 'format_def':\n", "/tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.c:138:18: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]\n", " if (size < sizeof(notalloc)) {\n", " ^\n", "\u001b[39mcompiling Fortran sources\u001b[0m\n", "\u001b[39mFortran f77 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran f90 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran fix compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmp32w1xbxr/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgfortran:f77: hellofortran.f\u001b[0m\n", "\u001b[39m/home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmp32w1xbxr/tmp/tmp32w1xbxr/src.linux-x86_64-3.6/hellofortranmodule.o /tmp/tmp32w1xbxr/tmp/tmp32w1xbxr/src.linux-x86_64-3.6/fortranobject.o /tmp/tmp32w1xbxr/hellofortran.o -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -lgfortran -o ./hellofortran.cpython-36m-x86_64-linux-gnu.so\u001b[0m\n", "Removing build directory /tmp/tmp32w1xbxr\n" ] } ], "source": [ "!source deactivate pangeo; 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(4)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 1 Fortran says hello\n", " 2 Fortran says hello\n", " 3 Fortran says hello\n", " 4 Fortran says hello\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", " integer 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...\n", "\tReading file 'dprod.f' (format:fix,strict)\n", "Post-processing...\n", "\tBlock: dprod\n", "\t\t\tBlock: dprod\n", "Post-processing (stage 2)...\n", "Saving signatures to file \"./dprod.pyf\"\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 `dprod.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 -*-\n", "! Note: the context of this file is case sensitive.\n", "\n", "python module dprod ! in \n", " interface ! in :dprod\n", " subroutine dprod(x,y,n) ! in :dprod:dprod.f\n", " double precision dimension(n) :: x\n", " double precision :: y\n", " integer, optional,check(len(x)>=n),depend(x) :: n=len(x)\n", " end subroutine dprod\n", " end interface \n", "end python module dprod\n", "\n", "! This file was auto-generated with f2py (version:2).\n", "! See http://cens.ioc.ee/projects/f2py2e/\n" ] } ], "source": [ "!cat dprod.pyf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The module does not know what Fortran subroutine arguments are 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": {}, "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/tmpsds5fq9n/src.linux-x86_64-3.6\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/tmpsds5fq9n/src.linux-x86_64-3.6/dprodmodule.c\"\n", "\u001b[39m adding '/tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.c' to sources.\u001b[0m\n", "\u001b[39m adding '/tmp/tmpsds5fq9n/src.linux-x86_64-3.6' to include_dirs.\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpsds5fq9n/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpsds5fq9n/src.linux-x86_64-3.6\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[39mget_default_fcompiler: matching types: '['gnu95', 'intel', 'lahey', 'pg', 'absoft', 'nag', 'vast', 'compaq', 'intele', 'intelem', 'gnu', 'g95', 'pathf95', 'nagfor']'\u001b[0m\n", "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n", "\u001b[39mFound executable /home2/nhn2/miniconda3/envs/pangeo/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: gcc -pthread -B /home2/nhn2/miniconda3/envs/pangeo/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC\n", "\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsds5fq9n/tmp\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsds5fq9n/tmp/tmpsds5fq9n\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsds5fq9n/tmp/tmpsds5fq9n/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmpsds5fq9n/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgcc: /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/dprodmodule.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/dprodmodule.c:16:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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/tmpsds5fq9n/src.linux-x86_64-3.6/dprodmodule.c:109:12: warning: 'f2py_size' defined but not used [-Wunused-function]\n", " static int f2py_size(PyArrayObject* var, ...)\n", " ^~~~~~~~~\n", "\u001b[39mgcc: /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.c:2:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.c: In function 'format_def':\n", "/tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.c:138:18: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]\n", " if (size < sizeof(notalloc)) {\n", " ^\n", "\u001b[39mcompiling Fortran sources\u001b[0m\n", "\u001b[39mFortran f77 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran f90 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran fix compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmpsds5fq9n/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgfortran:f77: dprod.f\u001b[0m\n", "\u001b[39m/home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmpsds5fq9n/tmp/tmpsds5fq9n/src.linux-x86_64-3.6/dprodmodule.o /tmp/tmpsds5fq9n/tmp/tmpsds5fq9n/src.linux-x86_64-3.6/fortranobject.o /tmp/tmpsds5fq9n/dprod.o -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -lgfortran -o ./dprod.cpython-36m-x86_64-linux-gnu.so\u001b[0m\n", "Removing build directory /tmp/tmpsds5fq9n\n" ] } ], "source": [ "!source deactivate pangeo; f2py -c dprod.pyf dprod.f; source activate pangeo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Using the module from Python" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import dprod" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on module dprod:\n", "\n", "NAME\n", " dprod\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", " dprod = \n", "\n", "VERSION\n", " b'$Revision: $'\n", "\n", "FILE\n", " /home/home2/nhn2/notebooks/dprod.cpython-36m-x86_64-linux-gnu.so\n", "\n", "\n" ] } ], "source": [ "help(dprod)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6.082818640342675e+62" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dprod.dprod(arange(1.0,50.0))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6.082818640342675e+62" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compare to numpy\n", "prod(arange(1.0,50.0))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "120.0" ] }, "execution_count": 17, "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": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xvec = rand(500)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "933 ns ± 22.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "timeit dprod.dprod(xvec)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.16 µs ± 7.07 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\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": 21, "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": 22, "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": 23, "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/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.c\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsj8xct9a/src.linux-x86_64-3.6\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/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.c\"\n", "\u001b[39m adding '/tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.c' to sources.\u001b[0m\n", "\u001b[39m adding '/tmp/tmpsj8xct9a/src.linux-x86_64-3.6' to include_dirs.\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpsj8xct9a/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcopying /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpsj8xct9a/src.linux-x86_64-3.6\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[39mget_default_fcompiler: matching types: '['gnu95', 'intel', 'lahey', 'pg', 'absoft', 'nag', 'vast', 'compaq', 'intele', 'intelem', 'gnu', 'g95', 'pathf95', 'nagfor']'\u001b[0m\n", "\u001b[39mcustomize Gnu95FCompiler\u001b[0m\n", "\u001b[39mFound executable /home2/nhn2/miniconda3/envs/pangeo/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: gcc -pthread -B /home2/nhn2/miniconda3/envs/pangeo/compiler_compat -Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC\n", "\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsj8xct9a/tmp\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsj8xct9a/tmp/tmpsj8xct9a\u001b[0m\n", "\u001b[39mcreating /tmp/tmpsj8xct9a/tmp/tmpsj8xct9a/src.linux-x86_64-3.6\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmpsj8xct9a/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgcc: /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.c:16:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.c:109:12: warning: 'f2py_size' defined but not used [-Wunused-function]\n", " static int f2py_size(PyArrayObject* var, ...)\n", " ^~~~~~~~~\n", "\u001b[39mgcc: /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.c\u001b[0m\n", "In file included from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1818:0,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,\n", " from /home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n", " from /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.h:13,\n", " from /tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.c:2:\n", "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-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/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.c: In function 'format_def':\n", "/tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.c:138:18: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]\n", " if (size < sizeof(notalloc)) {\n", " ^\n", "\u001b[39mcompiling Fortran sources\u001b[0m\n", "\u001b[39mFortran f77 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran f90 compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n", "Fortran fix compiler: /home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\u001b[0m\n", "\u001b[39mcompile options: '-I/tmp/tmpsj8xct9a/src.linux-x86_64-3.6 -I/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/numpy/core/include -I/home2/nhn2/miniconda3/envs/pangeo/include/python3.6m -c'\u001b[0m\n", "\u001b[39mgfortran:f77: dcumsum.f\u001b[0m\n", "\u001b[39m/home2/nhn2/miniconda3/envs/pangeo/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmpsj8xct9a/tmp/tmpsj8xct9a/src.linux-x86_64-3.6/dcumsummodule.o /tmp/tmpsj8xct9a/tmp/tmpsj8xct9a/src.linux-x86_64-3.6/fortranobject.o /tmp/tmpsj8xct9a/dcumsum.o -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -L/usr/lib/gcc/x86_64-redhat-linux/6.4.1 -lgfortran -o ./dcumsum.cpython-36m-x86_64-linux-gnu.so\u001b[0m\n", "Removing build directory /tmp/tmpsj8xct9a\n" ] } ], "source": [ "!source deactivate pangeo; f2py -c dcumsum.f -m dcumsum; source activate pangeo" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import dcumsum" ] }, { "cell_type": "code", "execution_count": 25, "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": 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": [ "py_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": [ "dcumsum.dcumsum(a)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 1., 3., 6., 10., 15., 21., 28., 36.])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cumsum(a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Benchmark the different implementations:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = rand(10000)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.77 ms ± 220 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "timeit py_dcumsum(a)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11.6 µs ± 425 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "timeit dcumsum.dcumsum(a)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "28.6 µs ± 437 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" ] } ], "source": [ "timeit a.cumsum()" ] }, { "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%file functions.c\n", "\n", "#include \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": null, "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": null, "metadata": { "collapsed": false }, "outputs": [], "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": null, "metadata": { "collapsed": false }, "outputs": [], "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%file run_hello_c.py\n", "\n", "import functions\n", "\n", "functions.hello(3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "!python run_hello_c.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Product function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "functions.dprod([1,2,3,4,5]) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cumulative sum:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = rand(100000)\n", "b = np.asarray(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import xarray as xr\n", "da = xr.DataArray(a, dims=['x'])\n", "da" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res_c = functions.dcumsum(a, len(a)) " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res_fortran = dcumsum.dcumsum(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "res_c - res_fortran" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simple benchmarking" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timeit dcumsum.dcumsum(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timeit functions.dcumsum(a, len(a))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timeit a.cumsum()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "timeit b.cumsum()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "timeit da.cumsum('x')" ] }, { "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%file cy_dcumsum.pyx\n", "# ln -s ~/miniconda3/lib/python3.6/site-packages/numpy/core/include/numpy ~/miniconda3/envs/pangeo/include/python3.6m/\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": null, "metadata": { "collapsed": false }, "outputs": [], "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "!python setup.py build_ext --inplace" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import cy_dcumsum" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "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": null, "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "b = empty_like(a)\n", "cy_dcumsum.dcumsum(a, b)\n", "b" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "py_dcumsum(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "a = rand(100000)\n", "b = empty_like(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timeit py_dcumsum(a)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "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 `Cython`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%load_ext Cython" ] }, { "cell_type": "code", "execution_count": null, "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": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "timeit cy_dcumsum2(a,b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Further reading" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* http://cython.org" ] } ], "metadata": { "gist_info": { "gist_id": null, "gist_url": null }, "kernelspec": { "display_name": "Py3 pangeo", "language": "python", "name": "pangeo" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }