{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystencils.session import *\n",
    "import timeit\n",
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo: Benchmark numpy, Cython, pystencils\n",
    "\n",
    "In this notebook we compare and benchmark different ways of implementing a simple stencil kernel in Python.\n",
    "Our simple example computes the average of the four neighbors in 2D and stores it in a second array. To prevent out-of-bounds accesses, we skip the cells at the border and compute values only in the range `[1:-1, 1:-1]`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementations\n",
    "\n",
    "The first implementation is a pure Python implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_pure_python(src, dst):       \n",
    "    for x in range(1, src.shape[0] - 1):\n",
    "        for y in range(1, src.shape[1] - 1):\n",
    "            dst[x, y] = (src[x + 1, y] + src[x - 1, y] +\n",
    "                         src[x, y + 1] + src[x, y - 1]) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously, this will be a rather slow version, since the loops are written directly in Python. \n",
    "\n",
    "Next, we use *numpy* functions to delegate the looping to numpy. The first version uses the `roll` function to shift the array by one element in each direction. This version has to allocate a new array for each accessed neighbor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_numpy_roll(src, dst):\n",
    "    neighbors = [np.roll(src, axis=a, shift=s) for a in (0, 1) for s in (-1, 1)]\n",
    "    np.divide(sum(neighbors), 4, out=dst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using views, we can get rid of the additional copies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def avg_numpy_slice(src, dst):\n",
    "    dst[1:-1, 1:-1] = src[2:, 1:-1] + src[:-2, 1:-1] + \\\n",
    "                      src[1:-1, 2:] + src[1:-1, :-2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To further optimize the kernel we switch to Cython, to get a compiled C version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import cython\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def avg_cython(object[double, ndim=2] src, object[double, ndim=2] dst):\n",
    "    cdef int xs, ys, x, y\n",
    "    xs, ys = src.shape\n",
    "    for x in range(1, xs - 1):\n",
    "        for y in range(1, ys - 1):\n",
    "            dst[x, y] = (src[x + 1, y] + src[x - 1, y] +\n",
    "                         src[x, y + 1] + src[x, y - 1]) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If available, we also try the numba just-in-time compiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    from numba import jit\n",
    "\n",
    "    @jit(nopython=True)\n",
    "    def avg_numba(src, dst):\n",
    "        dst[1:-1, 1:-1] = src[2:, 1:-1] + src[:-2, 1:-1] + \\\n",
    "                          src[1:-1, 2:] + src[1:-1, :-2]\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we also create a *pystencils* version of the same stencil code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "src, dst = ps.fields(\"src, dst: [2D]\")\n",
    "\n",
    "update = ps.Assignment(dst[0,0], \n",
    "                       (src[1, 0] + src[-1, 0] + src[0, 1] + src[0, -1]) / 4)\n",
    "avg_pystencils = ps.create_kernel(update).compile()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_implementations = {\n",
    "    'pure Python': avg_pure_python,\n",
    "    'numpy roll': avg_numpy_roll,\n",
    "    'numpy slice': avg_numpy_slice,\n",
    "    'pystencils': avg_pystencils,\n",
    "}\n",
    "if 'avg_cython' in globals():\n",
    "    all_implementations['Cython'] = avg_cython\n",
    "if 'avg_numba' in globals():\n",
    "    all_implementations['numba'] = avg_numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark functions\n",
    "\n",
    "We implement a short function to get in- and output arrays of a given shape and to measure the runtime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_arrays(shape):\n",
    "    in_arr = np.random.rand(*shape)\n",
    "    out_arr = np.empty_like(in_arr)\n",
    "    return in_arr, out_arr\n",
    "\n",
    "def do_benchmark(func, shape):\n",
    "    in_arr, out_arr = get_arrays(shape)\n",
    "    func(src=in_arr, dst=out_arr) # warmup\n",
    "    timer = timeit.Timer('f(src=src, dst=dst)', globals={'f': func, 'src': in_arr, 'dst': out_arr})\n",
    "    calls, time_taken = timer.autorange()\n",
    "    return time_taken / calls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHkCAYAAAAuH2ukAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde5zWc/7/8cezopNOmqxxyKBy2BBGrJSi5LSlXZbypQhrEVk57Wad95vl57TLNxWyCruWrIjF6sAqTNFJtlghViJyqph6/f64PpNrrrnmmpmaZqY877fbdfO53p/34fUZ1216zfv9/nwuRQRmZmZmll292g7AzMzMrC5zsmRmZmaWg5MlMzMzsxycLJmZmZnl4GTJzMzMLAcnS2ZmZmY5NKjtAKxuyMvLi4KCgtoOw8zMrEbMnDnzk4hoU5m6TpYMgIKCAoqKimo7DDMzsxoh6d3K1vUynJmZmVkOTpbMzMzMcnCyZGZmZpaDkyUzMzOzHLzB2wCY+8EKCi57srbDMDMzy2rxiGNqbWzPLJmZmZnl4GTJzMzMLAcnS2ZmZmY5OFkyMzMzy8HJkpmZmVkOm0SyJGmNpNclzZP0sKQmOeoWSBqQ9n6QpD9tpLiOk/S75PiWJMbXJS2U9HlavYGSFiWvgWnlUyT9O63dNlnG6Jx2frakflWIb0tJ0yT5rkczM7P1VGv/iEpqEBHFlay+MiI6Je3GA2cDN5dTtwAYADywwUFW7BKgD0BEXFhSKGkIsG9yvDVwJVAIBDBT0uMR8VlS/eSIyPWlbPOAwogolpQPzJY0sTI/u4j4VtI/gROB8VW/PDMzM1vvmaVkBudNSfdJmiPpbyUzPpIWS8pLjgslTUmOr5I0StIzwJ8l1Zd0o6RXkz5+WYmhXwDaSbpW0gVp8Vwv6XxgBNA1mYkpSWC2k/R0MrPzh7Q2/SXNTWasbkgr/yrpb7akGZJ+lOX6OwCrI+KTLDH2Bx5MjnsDz0bE8iRBehY4shLXCUBEfJOWGDUilXBlxrJTcm15kupJekHSEcnpx4CTKzuemZmZlbahy3C7AaMiYm/gC+CcSrTZH+gbEQOAwcCKiDgAOAA4U9LO5TVMlpOOAuYCdwMDk/J6wEmkZk8uA16IiE4RcUvStBOp2ZW9gBMl7ShpO+AG4LDk/AGSjkvqNwVmRMQ+wDTgzCzhdAFmZYlxJ2Bn4PmkaHvg/bQqS5KyEvcmid0VklTOdR8oaX5y3WdnzipFxLvJtYwELgLeiIhnktPzSP1szczMbD1saLL0fkT8KzkeBxxSiTaPR8TK5PgI4FRJrwMvA62B9lnaNE7qFAHvAXdHxGLgU0n7Jv28FhGfljPmPyNiRUSsAt4AdiKVQEyJiGVJ8jEe6JbU/xZ4IjmeSWppL1M+sCxL+UnA3yJiTfI+WwJUMjt0ckTsBXRNXqdkCz4iXo6IHycxXy6pUZY6Y4BmpJYoh6WVrwG+ldQss42ksyQVSSpa882KbEObmZn94G3onqXMJaGS98V8n4hl/sP+ddqxgCER8Y8Kxlm3ZynDGGAQsC1wT472q9OO15C67qyzOInvIiIy6peJCWiRpfwk4Ny090uA7mnvdwCmAETEB8l/v5T0ANAZ+HN5QUXEAklfAx1JJY7rJEugOyRvtwK+TDvdEFiVpb9RwCiAhvntyyzvmZmZ2YbPLLWV9JPkuD/wYnK8mNRyG8DPc7T/B/ArSVtAah+QpKZVGH8Cqf0/ByR9QSpJKDOLksXLwKHJPp/6SfxTqzD2AqBdeoGk3YBWwPS04n8AR0hqJakVqVmwf0hqkLavawvgWFJLZmT0uXPJ3WzJEt9upH6+mW4gNTv2O2B0WvvWwLKI+K4K12ZmZmaJDU2WFgADJc0Btgb+Lym/GrhN0gukZmbKM4bUstgsSfOAu6jCbFdEfAtMBv6atuw1ByhONmdfmKPtf4HLk/azgVkR8ffKjk1qL9O+GfuM+gMPpc1KERHLgWuBV5PXNUlZQ1JJ0xzgdeADkiRHUh9J1yRdHELqDrjXSSWH52RuKpd0KKmE8YaIGE9q2e205HQPYFIVrsvMzMzSKO3f9ao1lAqAJyKiY3UGVMUY6pHaZH1CRCyqhfFvAyZGxHM1PXZlSXoUuDwi/p2rXsP89pE/8NYaisrMzKxqFo84plr7kzQzIgorU3eTeChlNpL2BN4itXm7xhOlxO+Bch+QWdskbQk8VlGiZGZmZuVb7w3eyd1otTarFBFvALvU1vhJDEuBx2szhlySZcpyN4ybmZlZxTbZmSUzMzOzmuBkyczMzCwHf8GqAbDX9i0oqubNc2ZmZpsDzyyZmZmZ5eBkyczMzCwHJ0tmZmZmOThZMjMzM8vBG7wNgLkfrKDgsidrOwwz28iq+ynIZj8EnlkyMzMzy8HJkpmZmVkOTpbMzMzMcnCyZGZmZpaDkyUzMzOzHJwsVYGk7pKeSI4HSfpTOfWOk/S7jLLjJYWkwuR9a0mTJX1VXj9JvX0kTZc0V9JESc2rEO9eksZWtr6ZmZmV5WQpg6TqeJzCJcCdaX02A84HXk6rswq4AhhWQV9jgMsiYi9gAnBxZYOIiLnADpLaVraNmZmZlVajyZKkAkkLJI2WNF/SM5IaJ+empM265ElanBwPkvRYMqvyjqTzJP1a0muSZkjaOq39rZJekjRPUmdJ9SQtktQmqVNP0luS8jLiukrSKEnPAH+W1EjSvclszmuSelThGjsAqyPik7Tia4E/kEqQAIiIryPixfSycuwGTEuOnwV+nmXMfpKeU0q+pIWStk1OTwROqmz8ZmZmVlptzCy1B+6IiB8Dn5PlH/8sOgIDgM7A9cA3EbEvMB04Na1e04g4GDgHuCci1gLjgJOT8z2B2RmJTIn9gb4RMQA4FyCZzekP3CepUSWvrwswq+SNpH2BHSPiiUq2zzQP6JMcnwDsmFkhIiYAHyVxjwaujIiPktNFQNf1HNvMzOwHrzaSpXci4vXkeCZQUIk2kyPiy4hYBqwgNVsCMDej/YMAETENaC6pJXAP3ydUpwP3ljPG4xGxMjk+BLg/6etN4F2gQyXiBMgHlkFqJgu4Bbiokm2zOR04V9JMoBnwbTn1hgCXk5rVejCt/GNgu2wNJJ0lqUhS0ZpvVmxAiGZmZpuv2kiWVqcdr+H7r1wp5vt4Mmdx0tusTXu/ltJf2RIZ7SIi3geWSjoMOBB4qpy4vk47VrnRV2wl38ffjNSs2JRkWfEg4PGS5cbKiIg3I+KIiNifVDL4djlVtyf18/hRkqSVaJTElK3vURFRGBGF9Zu0qGxIZmZmPyh1aYP3YlJLYQDHr2cfJwJIOgRYEREl0yVjSC3H/TUi1lSin2kkS3fJHqS2wL8rGcMCoB1ARKyIiLyIKIiIAmAG0CciiirZF5K2Sf5bDxgOjMxSpwGpGbMByfi/TjvdgdRSnpmZma2HupQs3QT8StJLQF5FlcvxWdJ+JDA4rfxxYCvKX4LLdCdQX9Jc4C/AoIhYXUGbEtOAfSVVODuVzDbdDAyStETSnkn5mLTZp/6SFgJvAh+Wcw2/AV6IiBdIJUpnSNojOdcD8DfkmpmZrSdFZK5cbZokTQGGZZu1SRKPWyKiRjY6S7oNmBgRz9XEeDniaAhMBQ6JiOJcdRvmt4/8gbfWTGBmVmsWjzimtkMwqxMkzYyISm2LqUszSxuFpMuAR0htfq4pvwea1OB45WlL6hlNORMlMzMzK191PICxToiI7uWUjwBG1HAsS0kt/dWqiFgELKrtOMzMzDZlm/3MkpmZmdmGcLJkZmZmloOTJTMzM7McNps9S7Zh9tq+BUW+S8bMzKwMzyyZmZmZ5eBkyczMzCwHJ0tmZmZmOXjPkgEw94MVFFzmb0UxP+HZzCyTZ5bMzMzMcnCyZGZmZpaDkyUz22Q99thjdOvWjW222YbGjRuz0047cdxxx/H000+vqzN27FgksXjx4hqNbciQIfz0pz9d9/7dd9+lb9++7LTTTjRu3Ji8vDy6d+/OU089VapdUVERZ511FrvvvjtNmjShbdu2nHzyybzzzjvrHctll13G3nvvTcuWLWnSpAm777471157Ld988826Ol988QXXXHMNBx98MK1bt6Zly5YcfPDBPPbYY2X669u3L+eee+56x2O2qXGyZGabpNtvv51+/frRvn177r77bp588kmGDx8OwPPPP7+u3jHHHMP06dPJz8+vsdjefvtt7rrrLq688sp1ZV999RV5eXlcd911TJo0ibvvvputttqKo48+mkcffXRdvYceeoj58+dz/vnn89RTTzFixAhmzZpFYWEh77///nrF88UXX3DaaafxwAMPMHHiRE4++WSuv/56+vfvv67Oe++9x5133smhhx7KuHHj+Mtf/kKHDh3o168fd9xxR6n+rrrqKkaPHs3ChQvXKx6zTY0iorZjsDqgYX77yB94a22HYXXAprLBu23btuy///5MmDChzLm1a9dSr17t/S04ZMgQZsyYwauvvpqzXnFxMTvvvDOdOnVi4sSJACxbtow2bdqUqvfuu++y8847M3z4cK655ppqifHyyy9nxIgRLFu2jLy8PL7++msk0aRJk1L1Dj/8cBYtWsR7771Xqrxz584UFhZy5513Vks8ZjVN0syIKKxMXc8slUNSgaR5yXGhpNurqd9bJXVLjiXpekkLJS2QdH45bQZKWpS8BlZxvOcktaqO2M3qkuXLl7PttttmPZeeKGUuww0aNAhJWV9TpkxZ12727Nn06dOHVq1a0bhxY7p06cILL7xQYVyrV69m3LhxDBgwoMK6DRo0oEWLFmyxxRbryjITJYCddtqJNm3a8MEHH1TYZ2W1bt0aYN3YTZs2LZMoARQWFvLhhx+WKT/ppJMYP348K1eurLaYzOoqJ0uVEBFFEZE1kakKSVsDB0XEtKRoELAjsHtE7AE8VE6bK4EDgc7AlVVMfu4HztmQuM3qos6dO3Pfffdx4403Vmk56IorrmD69OmlXl26dFm3Pwhg1qxZHHzwwSxfvpzRo0fzyCOP0Lp1a3r27MnMmTNz9j9jxgw+//xzunbtmvX82rVrKS4u5qOPPuLaa69l4cKFFe7/WbBgAR9//DF77LFHpa8zm+LiYr766iuee+45br75Zk4//XRatGiRs820adPYfffdy5R369aNL774gunTp29QTGabghpNlpLZmgWSRkuaL+kZSY2Tc1MkFSbHeZIWJ8eDJD0maaKkdySdJ+nXkl6TNCNJJkra3yrpJUnzJHWWVC+ZjWmT1Kkn6S1JeRlxHSrp9eT1mqRmGee7S3oiOd5K0r2S5kqaI+nnSfkRkqZLmiXpYUlbZfkRHA88nfb+V8A1EbEWICI+ztKmN/BsRCyPiM+AZ4EjM+JrIenfknZL3j8o6czk9ONAf8w2MyNHjqRdu3Zccskl7LbbbuTl5dG/f3+eeeaZnO123XVXDjrooHWvF198kenTpzN+/Hh22WUXAC6++GLatm3L888/z/HHH8/RRx/NhAkT2GWXXbj22mtz9j9jxgwksffee2c9f8kll7DFFluQn5/PH/7wBx566CEOP/zwcvsrLi7m7LPPpk2bNgwePLiCn0r55s2bxxZbbEGzZs3o1asXvXr1YtSoUTnbjBo1ihkzZnD55ZeXObfPPvtQr149ZsyYsd4xmW0qamNmqT1wR0T8GPgc+Hkl2nQEBpCaWbke+CYi9gWmA6em1WsaEQeTmkm5J0lCxgEnJ+d7ArMj4pOM/ocB50ZEJ6ArkGte+QpgRUTsFRF7A88nyddwoGdE7AcUAb/O0rYLkP5n6a7AiZKKJD0lqX2WNtsD6bs6lyRl60TECuA8YKykk4BWETE6OfcZ0FBS6xzXZLbJ6dChA6+99hpTp07lt7/9LZ06dWLChAn07t2b6667rlJ9TJw4kUsvvZQbbriB4447DoCVK1cydepUTjjhBOrVq0dxcTHFxcVEBD179mTatGk5+/zwww9p3rw5W265ZdbzQ4cO5dVXX2XixIkcddRRDBgwgCeeeKLc/s477zxeeuklxo0bR6tW67+i3q5dO1599VWmTJnC73//eyZMmMCpp55abv0pU6Zw/vnnc8opp3DyySeXOb/FFlvQokWLrEt0Zpub2niC9zsR8XpyPBMoqESbyRHxJfClpBXAxKR8LpD+59uDABExTVJzSS2Be4C/A7cCpwP3Zun/X8DNksYDj0bEEknlxdITOKnkTUR8JulYYE/gX0m7LUklcpnygWVp7xsCqyKiUNLPklgz5+6zBVJmV35EPCvpBOAOYJ+M0x8D2wGflupYOgs4C6B+87L7JMzquvr169OtWze6desGpBKVI488kquvvppzzz03Z3Ixe/ZsBgwYwODBgxk2bNi68uXLl7NmzRquvfbacmeRcm0gX7VqFQ0bNix33B122IEddtgBgGOPPZbu3bszbNgwjj322DJ1L7/8ckaNGsV9993HEUccUW6fldGoUSMKC1N7WQ899FDy8/M57bTTGDJkCAcddFCpuq+++ip9+vThsMMO4+677y63z8aNG3vPkv0g1MbM0uq04zV8n7AV8308jXK0WZv2fi2lE77MJCIi4n1gqaTDSO37eYqylUYAZwCNgRmSyi7Qf09ZxhGppbJOyWvPiMg2X76S0te2BHgkOZ5A6cQvvc6Oae93AMr8KSepHrBHMsbWGacbkWW2LCJGRURhRBTWb5J734LZpmC77bbjjDPOoLi4mEWLFpVbb+nSpfTp04eDDjqozN1cLVu2pF69egwZMoRXX3016yvXnXatW7fms88+q3TMhYWFvPXWW2XKr7/+ekaMGMFtt93GKaecUun+qjIuUGbsuXPn0rt3bzp16sQjjzxSavN5puXLl5OXl1fuebPNRV3a4L0Y2D85Pn49+zgRQNIhpJbKViTlY0gtx/01ItZkNpK0a0TMjYgbSC2h5UqWniG15FXSthUwA+giqV1S1kRShyxtFwDt0t4/BhyWHB8KZNul+g/gCEmtkrGOSMoyXZj03x+4R9IWSSwCtiX18zXbbJT3zKE333wToNw75VatWkXfvn1p2rQpDz/8MA0alJ5gb9q0KV27dmX27Nnst99+FBYWlnnlsvvuu/Pdd9+xZMmSCq9h7dq1vPjii+y6666lym+//XaGDx/O9ddfz5AhQyrsZ31MnToVoNTYixYtolevXuyyyy488cQTNG7cuNz2H330EatWrWK33XbbKPGZ1SV16Yt0bwL+KukU4PmKKpfjM0kvAc1JLbmVeJzU8lu2JTiAoZJ6kJrpeoPU7FN5T7C7DrgjeazAGuDqiHhU0iDgQUkl8+/DKZv8PAn8klTyBjACGC/pQuArUrNbJBvdz46IMyJiuaRrgZIHtlwTEcvTO00SszOAzhHxpaRpyfhXkkpAZ0REcTnXY7ZJ6tixIz169KBfv37svPPOfPHFF0yaNImRI0fyi1/8Yt2dbZmGDh3KrFmzGDt27LrEqsSee+5J8+bNufnmm+nWrRu9e/dm8ODB5Ofn88knnzBr1izWrFnDiBEjyo2rZEnwlVdeWbfcBqkHOS5fvpwuXbqw7bbb8tFHH3H33Xfzyiuv8MADD6yr99BDDzF06FCOPPJIDjvssFIbqJs3b86ee+657n337t1ZvHhxzqeTz5kzh2HDhnHCCSewyy67sHr1aqZNm8Ztt93GUUcdxU9+8hMAPv74Y3r16sW3337L1VdfzRtvvFGqn3333bfU8uLLL79c6nrNNmebzUMpJU0BhkVEUZZzhcAtEZH9Xt4aJOlF4NiI+LyGxrsNeDwi/pmrnh9KaSU2lYdSjhw5kkmTJjF79myWLl1K/fr16dChA/3792fo0KHrNliPHTuW0047jXfeeYeCggK6d+++blYl0+TJk+nevTuQul3/6quv5vnnn2fFihW0adOG/fbbj7PPPpujjz46Z2wHHngge+65J/fe+/3fZ48//ji33nor8+bNY8WKFWy77bbss88+XHrppXTp0mVdvUGDBnHfffdl7ffQQw8t9SyoAw44gPr16+e8I23p0qVceOGFTJ8+nY8++ogmTZqwyy67MGjQIM4444x1CdCUKVPo0aNHuf2U/PxKnHnmmbz22msUFZX5lWu2SVAVHkq52SdLki4jdYv+yRHxYm3ElhHPgcDKiJhTQ+OdWXJnXC5OlqzEppIs1WVjx47lggsu4L///W/WBz1Wh6+//ppWrVoxbtw4fvGLX2yUMcqzatUq8vPzuemmmzbocQZmtakqyVJd2rO0QSKie7ZZpYgYERE71YVECSAiXq6pRCkZr8JEycyq1ymnnML222+/Ub8K5KWXXmLXXXfl+OPXd4vn+rvrrrvYZpttGDiwSl8oYLbJ2mySJTOzuqJ+/frcc889G21WCaBXr14sWLCgVr4Dr2HDhowdO7bM5nizzdVmswxnG8bLcFbCy3Bm9kPwg1yGMzMzM9sYPIdqAOy1fQuKPKNgZmZWhmeWzMzMzHJwsmRmZmaWg5MlMzMzsxycLJmZmZnl4A3eBsDcD1ZQcNmTtR3GJsu325uZbb48s2RmZmaWg5MlMzMzsxycLJmZmZnl4GTJzMzMLIfNLlmStK2khyS9LekNSZMkdchRv5Oko9PeXyVp2EaKbaikU5PjGyW9KWmOpAmSWmapv5uk19NeX0gaWoXxbpJ0WHVeg5mZ2Q/NZpUsSRIwAZgSEbtGxJ7Ab4Af5WjWCTg6x/nqiq0BcDrwQFL0LNAxIvYGFgKXZ7aJiH9HRKeI6ATsD3xD6voq64/AZRsUuJmZ2Q/cZpUsAT2A7yJiZElBRLweES9Iul9S35JySeMl9QGuAU5MZm5OTE7vKWmKpP9IOj+tza8lzUteQ5OyAkkLJI2WNF/SM5IaZ4ntMGBWRBQncT1TcgzMAHao4NoOB96OiHczT0j6e9qM1S8ljU/GeBdoLWnbCvo2MzOzcmxuyVJHYGY558YApwFIagEcDEwCfgf8JZnB+UtSd3egN9AZuFLSFpL2T9ofCBwEnClp36R+e+COiPgx8Dnw8yzjd8kR2+nAUxVc20nAg+WcOwv4naSuwEXAkLRzs5KxzczMbD1sbslSuSJiKtBO0jZAf+CRtJmdTE9GxOqI+AT4mNQy3iHAhIj4OiK+Ah4Fuib134mI15PjmUBBlj7zgWWZhZJ+CxQD48uLXdKWQB/g4XKubSmppG8ycFFELE87/TGwXTn9niWpSFLRmm9WlDe8mZnZD9rmlizNJ7W3pzz3AyeTmiG6N0e91WnHa0g96VxVrJ9pJdAovUDSQOBY4OSIiBz9H0VqCW9pjjp7AZ9SNjFqlIxdRkSMiojCiCis36RFjq7NzMx+uDa3ZOl5oKGkM0sKJB0g6dDk7VhgKEBEzE/KvgSaVaLvacBxkppIagr0A16oQmwLgHZpcR0JXAr0iYhvKmjbn/KX4JDUmVRCtS8wTNLOaac7APOqEKeZmZml2aySpWR2ph/QK3l0wHzgKuDD5PxSUklL+qzSZFIbutM3eGfrexapZOsV4GVgTES8VoXwngK6pb3/E6kk7dlk7JEAkraTNKmkkqQmQC9Sy35lSGoIjAZOj4gPSe1ZukcpW5BK0IqqEKeZmZmlUe7Vn81LknjMBfaLiBrfpCNpAnBJRCyqofH6kbrWKyqq2zC/feQPvLUGoto8+Yt0zcw2LZJmRkRhZepuVjNLuUjqCbwJ/LE2EqXEZaQ2eteUBsD/q8HxzMzMNjvZNiJvliLiOaBtLcfwb+DfNThe1rvnzMzMrPJ+MDNLZmZmZuvDyZKZmZlZDk6WzMzMzHL4wexZstz22r4FRb6jy8zMrAzPLJmZmZnl4GTJzMzMLAcnS2ZmZmY5OFkyMzMzy8EbvA2AuR+soOCyJ6utP3/9h5mZbS48s2RmZmaWg5MlMzMzsxycLJmZmZnl4GTJzMzMLAcnS2ZmZmY5bNbJkqRBkrbbCP2OkbRncrxYUl4V2t4qqVtyfJ6ktyRFrj4kDZS0KHkNrGKsz0lqVZU2ZmZm9r3NOlkCBgHVnixFxBkR8UZV20naGjgoIqYlRf8CegLvVtDmSuBAoDNwZRWTn/uBc6oaq5mZmaXUyWRJUoGkNyXdJ2mOpL9JaiLpcEkT0ur1kvSopPqSxkqaJ2mupAslHQ8UAuMlvS6psaT9JU2VNFPSPyTlJ/1MkXSDpFckLZTUNSmvL+mmpM85koak1S/MiLmppCclzU7iODHLpR0PPF3yJiJei4jFFfw4egPPRsTyiPgMeBY4MmPsFpL+LWm35P2Dks5MTj8O9K9gDDMzMytHnUyWErsBoyJib+ALUrMjzwN7SGqT1DkNuBfoBGwfER0jYi/g3oj4G1AEnBwRnYBi4I/A8RGxP3APcH3aeA0iojMwlNRMDsBZwM7Avkkc43PEeyTwYUTsExEdSUuK0nQBZlbppwDbA++nvV+SlK0TESuA84Cxkk4CWkXE6OTcZ0BDSa2rOK6ZmZlRt5Ol9yPiX8nxOOCQiAhSy0r/I6kl8BPgKeA/wC6S/ijpSFLJVabdgI7As5JeB4YDO6SdfzT570ygIDnuCYyMiGKAiFieI965QM9khqprksBkygeW5broLJSlLMoURDybxHAHcEbG6Y/Jshwp6SxJRZKK1nyTLVwzMzOry8lSZkJQ8v5e4H9ILS09HBHFyezJPsAU4FxgTJb+BMyPiE7Ja6+IOCLt/Orkv2v4/mtglCWO7MFGLAT2J5Ww/K+k32WpthJoVJn+0iwBdkx7vwPwYWYlSfWAPZIxts443Sgpz4x5VEQURkRh/SYtqhiWmZnZD0NdTpbaSvpJctwfeBEgIj4klSwMB8YCJHeS1YuIR4ArgP2Sdl8CzZLjfwNtSvqUtIWkH1cQwzPA2ZIaJG0yk5B1krvuvomIccBNaTGkWwC0q2DMTP8AjpDUKtnYfURSlunCpP/+wD2StkjiErAtsLiK45qZmRl1O1laAAyUNIfUTMn/pZ0bT2qZruSOtO2BKcny2ljg8qR8LDAyKa9PaoP1DZJmA68DB1cQwxjgPWBO0mZAjrp7Aa8kY/0WuC5LnSeB7iVvJJ0vaQmp2aI5ksYk5YUlx8nS37XAq8nrmszlQEkdSC29XY4DB6oAACAASURBVBQRLwDTSCWTkJrtmlGylGhmZmZVo9Q2oLpFUgHwRLJROtv5PwGvRcTdNRlXdZD0InBsRHxeQ+PdBjweEf/MVa9hfvvIH3hrtY27eMQx1daXmZlZdZM0MyIKK65Zt2eWspI0E9ib1KbvTdFFQNsaHG9eRYmSmZmZla9BxVVqXvLsoayzSslt/5usiHi5hscbXZPjmZmZbW42uZklMzMzs5rkZMnMzMwshzq5DGc1b6/tW1DkTdlmZmZleGbJzMzMLAcnS2ZmZmY5OFkyMzMzy8HJkpmZmVkOTpbMzMzMcqiTX3diNU/SMuDdtKIWwIosVfOAT2okqA1X3jXUxTHWt5+qtKts3Yrq5Tqf69ym8tmpic9NdY2zuXxucp3fVD434N8561N3Q+tsyO+cnSKiTQVjp0SEX36VeQGjyikvqu3YNvQa6uIY69tPVdpVtm5F9XKdr+DcJvHZqYnPTXWNs7l8bnKd31Q+N9X1/7Smxqgrn50NrVNTv3O8DGflmVjbAVSDmriG6hpjffupSrvK1q2oXq7z/tzU7Diby+emKuPUZf6dU/W6G1qnRj43XoazKpFUFJX8lmazdP7s2Prw58bWV3V+djyzZFU1qrYDsE2WPzu2Pvy5sfVVbZ8dzyyZmZmZ5eCZJTMzM7McnCyZmZmZ5eBkyczMzCwHJ0tWbSTVk3S9pD9KGljb8dimQVJ3SS9IGimpe23HY5sWSU0lzZR0bG3HYpsGSXskv2/+JulXlWnjZMkAkHSPpI8lzcsoP1LSvyW9JemyCrrpC2wPfAcs2VixWt1RTZ+bAL4CGuHPzQ9GNX12AC4F/rpxorS6pjo+NxGxICLOBn4BVOrRAr4bzgCQ1I3UP1h/joiOSVl9YCHQi9Q/Yq8C/YH6wP9mdHF68vosIu6S9LeIOL6m4rfaUU2fm08iYq2kHwE3R8TJNRW/1Z5q+uzsTeorLRqR+hw9UTPRW22pjs9NRHwsqQ9wGfCniHigonEbVN8l2KYsIqZJKsgo7gy8FRH/AZD0ENA3Iv4XKDPlLWkJ8G3yds3Gi9bqiur43KT5DGi4MeK0uqeafuf0AJoCewIrJU2KiLUbNXCrVdX1OyciHgcel/Qk4GTJNsj2wPtp75cAB+ao/yjwR0ldgWkbMzCr06r0uZH0M6A30BL408YNzeq4Kn12IuK3AJIGkcxQbtTorK6q6u+c7sDPSP1xNqkyAzhZslyUpazcdduI+AYYvPHCsU1EVT83j5JKtM2q9NlZVyFibPWHYpuQqv7OmQJMqcoA3uBtuSwBdkx7vwPwYS3FYpsOf25sffmzY+tjo39unCxZLq8C7SXtLGlL4CTg8VqOyeo+f25sffmzY+tjo39unCwZAJIeBKYDu0laImlwRBQD5wH/ABYAf42I+bUZp9Ut/tzY+vJnx9ZHbX1u/OgAMzMzsxw8s2RmZmaWg5MlMzMzsxycLJmZmZnl4GTJzMzMLAcnS2ZmZmY5OFkyMzMzy8HJkpnVCklrJL0uaZ6kiZJabkBf3SUdnPb+bEmnVk+klRq/q6T5yfU0rkT9MZL2XM+xCiTNq2Kbl9LaDlifcXP0/ZtsY5ltTvycJQMgLy8vCgoKajsMMzOzGjFz5sxPIqJNZer6i3QNgIKCAoqKimo7DDMzsxoh6d3K1vUynJmZmVkOTpbMzMzMcnCyZGZmZpaD9ywZAHM/WEHBZU/WdhhmZmZZLR5xTK2N7ZklMzMzsxycLJmZmZnl4GTJzMzMLAcnS2ZmZmY5OFkyMzMzy6HOJEsZ3xP1sKQmOeqW+n4jSYMk/WkjxXWcpN8lx90kzZJULOn4tDo9kthLXqskHZecO0/SW5JCUl6W/g9Irv34LOeaZfT7iaRbqxD7eZJOW78rNzMzM9jIyZKkqjyaYGVEdIqIjsC3wNk56hYA1fplkDlcAtyZHL8HDAIeSK8QEZOT2DsBhwHfAM8kp/8F9ATKPFZdUn3gBuAf2QaOiC9L+k36fhd4tAqx3wOcX4X6ZmZmliFnspTM4Lwp6T5JcyT9rWTGR9LikpkSSYWSpiTHV0kaJekZ4M+S6ku6UdKrSR+/rERcLwDtJF0r6YK0eK6XdD4wAuiazLZcmJzeTtLTkhZJ+kNam/6S5iYzVjeklX+V9Ddb0gxJP8py/R2A1RHxCUBELI6IOcDaHLEfDzwVEd8kbV6LiMXl1B0CPAJ8XNEPRFJ7YJvkZ5N57va02a/ekqZJqpfEsFhS54r6NzMzs+wqM7O0GzAqIvYGvgDOqUSb/YG+ETEAGAysiIgDgAOAMyXtXF7DZDbqKGAucDcwMCmvB5wEjAcuA15IZlxuSZp2Ak4E9gJOlLSjpO1Izdwclpw/oGR5DGgKzIiIfYBpwJlZwukCzKrE9aY7CXiwokqStgf6ASMr2W9/4C8REVnOXUbqmnsAtwOnRURJQlcEdK3kGGZmZpahMsnS+xHxr+R4HHBIJdo8HhErk+MjgFMlvQ68DLQG2mdp0zipU0RquevuZEbmU0n7Jv28FhGfljPmPyNiRUSsAt4AdiKVnE2JiGURUUwq0eqW1P8WeCI5nklqaS9TPrCsEtcLgKR8Usla1mW1DLcCl0bEmkp2X24SlswgnQk8C/wpIt5OO/0xsF058Z4lqUhS0ZpvVlQyDDMzsx+WyuwpypzJKHlfzPfJVqOMOl+nHQsYEhEVJRArk305mcaQ2ie0Lak9OOVZnXa8htS1KUf979JmaUrql4kJaJGjj0y/ACZExHeVqFsIPCQJIA84WlJxRDyWWVHSPkCDiJiZo7+9gE8pmxg1InUdZUTEKGAUQMP89tlmrMzMzH7wKjOz1FbST5Lj/sCLyfFiUsttAD/P0f4fwK8kbQGpfUCSmlYhxgnAkaRmiUoSri+BZpVo+zJwqKS8ZDN1f2BqFcZeALSrQv3+VGIJDiAido6IgogoAP4GnJMtUapMv5J2Ai4C9gWOknRg2ukOwLzKxGRmZmZlVSZZWgAMlDQH2Br4v6T8auA2SS+QmpkpzxhSy2KzJM0D7qIKX+AbEd8Ck4G/pi1ZzQGKk83ZF+Zo+1/g8qT9bGBWRPy9smOT2su0r5Lpn+Q2/yXACcBdkuaXVJRUAOxIRjIm6fykzQ7AHEljKho0WY5M9wvKSZaS2O4GhkXEh6T2iI2RVDLb1wV4rqIxzczMLDtl3y+cnEwlAE8kt/PXimRj9yzghIhYVAvj3wZMjIhNLuFI9nr9OiJOqahuw/z2kT+w0o9wMjMzq1GLRxxTrf1JmhkRhZWpW2ceSpmNpD2Bt0ht3q7xRCnxe6DcB2TWcXnAFbUdhJmZ2aYs53JYcjdarc0qRcQbwC61NX4Sw1Lg8dqMYX1FxLO1HYOZmdmmrk7PLJmZmZnVNidLZmZmZjk4WTIzMzPLoSpfdGubsb22b0FRNd9pYGZmtjnwzJKZmZlZDk6WzMzMzHJwsmRmZmaWg5MlMzMzsxy8wdsAmPvBCgoue7K2w7DNXHV/XYGZWU3wzJKZmZlZDk6WzMzMzHJwsmRmZmaWg5MlMzMzsxycLFWBpO6SnkiOB0n6Uzn1jpP0u+S4m6RZkoolHZ9R72lJn5f0WU5fbSVNlvSapDmSjq5CvHtJGlvZ+mZmZlaWk6UMkqrjDsFLgDuT4/eAQcADWerdCJxSQV/Dgb9GxL7ASWn9Vigi5gI7SGpb2TZmZmZWWo0mS5IKJC2QNFrSfEnPSGqcnJsiqTA5zpO0ODkeJOkxSRMlvSPpPEm/TmZaZkjaOq39rZJekjRPUmdJ9SQtktQmqVNP0luS8jLiukrSKEnPAH+W1EjSvZLmJuP0qMI1dgBWR8QnABGxOCLmAGsz60bEP4EvK+gygObJcQvgwyxj9pP0nFLyJS2UtG1yeiKpJMvMzMzWQ23MLLUH7oiIHwOfAz+vRJuOwACgM3A98E0y0zIdODWtXtOIOBg4B7gnItYC44CTk/M9gdkliUyG/YG+ETEAOBcgIvYC+gP3SWpUyevrAsyqZN3KuAr4H0lLgEnAkMwKETEB+IhU3KOBKyPio+R0EdC1GuMxMzP7QamNZOmdiHg9OZ4JFFSizeSI+DIilgErSM2WAMzNaP8gQERMA5pLagncw/cJ1enAveWM8XhErEyODwHuT/p6E3gX6FCJOAHygWWVrFsZ/YGxEbEDcDRwv6Rs/9+GAJeTmtV6MK38Y2C7bB1LOktSkaSiNd+sqMaQzczMNh+1kSytTjtew/dPES/m+3gyZ3HS26xNe7+W0k8hj4x2ERHvA0slHQYcCDxVTlxfpx2r3OgrtpKy8W+IwcBfASJietJ3XpZ625P6efwoI5lqlMRURkSMiojCiCis36RFNYZsZma2+ahLG7wXk1oKAzg+R71cTgSQdAiwIiJKpkvGkFqO+2tErKlEP9NIlu6SPUhtgX9XMoYFQLuqBF2B94DDk1j2IJX8lJq5Sjal30tqqXIB8Ou00x2AedUYj5mZ2Q9KXUqWbgJ+Jeklss+cVMZnSfuRpGZkSjwObEX5S3CZ7gTqS5oL/AUYFBGrK2hTYhqwryQBSDog2W90AnCXpPklFSW9ADwMHC5piaTeSfk1kvok1S4CzpQ0m9Qy46CIyJxB+w3wQkS8QCpROiNJrAB6AP7SNzMzs/Wksv/ubpokTQGGRURRlnOFwC0RUSMbnSXdBkyMiOdqYrwccTQEpgKHRERxrroN89tH/sBbayYw+8HyF+maWV0haWZEFFambl2aWdooJF0GPEJq83NN+T3QpAbHK09b4LKKEiUzMzMrX3U8gLFOiIju5ZSPAEbUcCxLSS391aqIWAQsqu04zMzMNmWb/cySmZmZ2YZwsmRmZmaWg5MlMzMzsxw2mz1LtmH22r4FRb5TyczMrAzPLJmZmZnl4GTJzMzMLAcnS2ZmZmY5OFkyMzMzy8EbvA2AuR+soOAyf4VcXeevCzEzq3meWTIzMzPLwcmSmZmZWQ5Olsysznnsscfo1q0b22yzDY0bN2annXbiuOOO4+mnn15XZ+zYsUhi8eLFNRrbkCFD+OlPf1qq7De/+Q1HHHEErVu3RhJjx44t0+6///0vl19+OYWFhbRo0YI2bdpw+OGHM23atDJ116xZwy233ELHjh1p2rQp+fn59OvXjzlz5qxXzF9++SXDhg2je/fuNG/eHElMmTKlTL2FCxdywQUXsPfee7PVVluRn59Pnz59mD17dpm633zzDVdeeSUdOnSgcePG7Ljjjpx66qll/n/07duXc889d73iNqsrnCyZWZ1y++23069fP9q3b8/dd9/Nk08+yfDhwwF4/vnn19U75phjmD59Ovn5+TUW29tvv81dd93FlVdeWar8j3/8IytXruTYY48tt+3MmTP5y1/+Qt++ffnb3/7G2LFjadSoEd27d+eJJ54oVfeKK65g2LBhHHfccUycOJHbbruNt99+mx49erBkyZIqx/3pp59yzz330KBBA3r16lVuvWeeeYbJkyczcOBAJk6cyJ133smyZcs48MADmTlzZqm6Z5xxBjfeeCNnnnkmkyZN4rrrrmPatGkcfvjhfPXVV+vqXXXVVYwePZqFCxdWOW6zukIRUdsxWB3QML995A+8tbbDsAr8EDZ4t23blv33358JEyaUObd27Vrq1au9v/GGDBnCjBkzePXVV0uVl8T11ltv0b59e+69914GDRpUqs7nn3/OVlttRYMG399XU1xczI9//GN+9KMflZph2m677ejevTsPPPDAurI333yTPfbYg5EjR/LLX/6ySnFHBJIAeO655+jVqxeTJ0+me/fupep98skn62bHSqxYsYKCggJ++tOf8uc//xmAlStX0qxZMy655BJ+//vfr6v79NNPc9RRR/H000/Tu3fvdeWdO3emsLCQO++8s0pxm21MkmZGRGFl6npmqRySCiTNS44LJd1eTf3eKqlbRtkfJX1VTv0tJd0raa6k2ZK6V3G85yS12oCQzWrU8uXL2XbbbbOeS0+UMpfhBg0ahKSsr/Qlp9mzZ9OnTx9atWpF48aN6dKlCy+88EKFca1evZpx48YxYMCAnHGVp2XLlqUSJYAGDRrQqVMnPvjgg1Ll3377Lc2bNy/THlKJWVWlJz+55OXllanbokULOnToUCrG4uJi1qxZU+kYTzrpJMaPH8/KlSurHLtZXeBkqRIioigizt/QfiRtDRwUEdPSygqBljmanZnEsBfQC/h/kqry/+1+4Jz1CNesVnTu3Jn77ruPG2+8sUpLN1dccQXTp08v9erSpQtNmjShbdu2AMyaNYuDDz6Y5cuXM3r0aB555BFat25Nz549yywzZZoxYwaff/45Xbt23aDrS/ftt98yffp09thjj1Ll55xzDuPGjePvf/87X3zxBf/5z38455xz2GGHHTjxxBOrbfzKWL58OfPmzSsVY7NmzTjllFO4/fbbmTx5Ml999RXz58/n4osvZp999uHwww8v1Ue3bt344osvmD59eo3GblZdajRZSmZrFkgaLWm+pGckNU7OTUkSByTlSVqcHA+S9JikiZLekXSepF9Lek3SjCQBKWl/q6SXJM2T1FlSPUmLJLVJ6tST9JakvIy4DpX0evJ6TVKzjPPdJT2RHG+VNtMzR9LPk/IjJE2XNEvSw5K2yvIjOB54Oq3f+sCNwCU5fmx7Av8EiIiPgc+BUtOGklpI+rek3ZL3D0o6Mzn9ONA/R/9mdcrIkSNp164dl1xyCbvttht5eXn079+fZ555Jme7XXfdlYMOOmjd68UXX2T69OmMHz+eXXbZBYCLL76Ytm3b8vzzz3P88cdz9NFHM2HCBHbZZReuvfbanP3PmDEDSey9997Vdq1XXXUVS5Ys4dJLLy1Vfs0113D55Zfzs5/9jBYtWrDrrrsyf/58pkyZwtZbb11t41fGkCFDiAiGDh1aqvzee++lX79+HHbYYTRr1oyOHTvy3Xff8eyzz7LllluWqrvPPvtQr149ZsyYUZOhm1Wb2phZag/cERE/JvUP/88r0aYjMADoDFwPfBMR+wLTgVPT6jWNiINJzaTcExFrgXHAycn5nsDsiPgko/9hwLkR0QnoCuSaK74CWBERe0XE3sDzSfI1HOgZEfsBRcCvs7TtAqT/+Xoe8HhE/DfHeLOBvpIaSNoZ2B/YMb1CRKxI+hor6SSgVUSMTs59BjSU1DrHGGZ1RocOHXjttdeYOnUqv/3tb+nUqRMTJkygd+/eXHfddZXqY+LEiVx66aXccMMNHHfccUBqn83UqVM54YQTqFevHsXFxRQXFxMR9OzZM+tdaek+/PBDmjdvXiYRWF8PPPAAI0aM4IorrigzW/V///d/XHfddQwfPpzJkyfz8MMP06xZM4444gg+/PDDahm/Mv73f/+XBx54gD/96U+0a9eu1Lnhw4czbtw4brrpJqZOncr999/Pp59+ylFHHcXXX39dqu4WW2xBixYtajR2s+pUG0/wficiXk+OZwIFlWgzOSK+BL6UtAKYmJTPBdL/zHsQICKmSWouqSVwD/B34FbgdODeLP3/C7hZ0njg0YhYkmONvydwUsmbiPhM0rGkZoD+lbTbklQilykfWAYgaTvgBKB77kvnHmAPUgnYu8BLQHFmpYh4VtIJwB3APhmnPwa2Az5NL5R0FnAWQP3mbSoIw6zm1K9fn27dutGtW2p734cffsiRRx7J1VdfzbnnnkurVuVvw5s9ezYDBgxg8ODBDBs2bF358uXLWbNmDddee225s0i5NpCvWrWKhg0bbsBVfW/ixIkMGjSIwYMHc/XVV5c6t3z5ci688EIuvvjiUucOO+wwCgoKuPHGG7nllluqJY5cRo4cyW9+8xuuu+46Tj/99FLn5s+fz4gRIxgzZgyDBw9eV37ggQfSoUMHxowZwwUXXFCqTePGjb1nyTZZtZEsrU47XgM0To6L+X6mq1GONmvT3q+l9DVk3toXEfG+pKWSDgMO5PtZpvRKIyQ9CRwNzJDUE1hVTvzKMo6AZyOiouWulXx/bfsC7YC3kgSriaS3IqLUn28RUQxcuG4g6SVgUZmgUvuY9kjG2BpIv7+4EVlmyyJiFDAKUnfDVRC7Wa3ZbrvtOOOMM7jgggtYtGgRnTt3zlpv6dKl9OnTh4MOOqjMnVctW7akXr16nHvuuZx66qlZ2+faqN26dWs+++yz9b+IxD//+U9OOOEE+vXrx1133VXm/MKFC1m9ejUHHHBAqfKtt96aXXfdlQULFmxwDBW5//77Oeecc7jooov47W9/W+b83LlzAcrE2L59e1q2bJk1xuXLl5OXl1em3GxTUJc2eC8mtcQEqb096+NEAEmHkFoqW5GUjyG1HPfXiFiT2UjSrhExNyJuIDWDs3uOMZ4hteRV0rYVMAPoIqldUtZEUocsbReQSpCIiCcjYtuIKIiIAlJLi+0yGyR9NU2OewHFEfFGlr4vTPrvD9wjaYukjYBtSf18zeq8999/P2v5m2++CVDunXKrVq2ib9++NG3alIcffrjMnWdNmzala9euzJ49m/3224/CwsIyr1x23313vvvuu/V6zlGJ6dOn07dvXw4//HDGjRuXNTkrub5XXnmlVPny5ct566232H777dd7/MqYMGECp512GmeccQY33XRT1jrlxbhw4UI+//zzMjF+9NFHrFq1it12223jBG22kdWlL9K9CfirpFOA5yuqXI7PkpmX5qSW3Eo8Tmr5LdsSHMBQST1IzXS9ATxFasksm+uAO5LHCqwBro6IRyUNAh6UVDJPPxzIvJXnSeCXpJK3cknqAxRGxO+AbYB/SFoLfACckqV+B+AMoHNEfClpWjL+laQS0BnJDJVZndexY0d69OhBv3792Hnnnfniiy+YNGkSI0eO5Be/+MW6O9syDR06lFmzZjF27Nh1iVWJPffck+bNm3PzzTfTrVs3evfuzeDBg8nPz+eTTz5h1qxZrFmzhhEjRpQbV8mS4CuvvMIOO+xQ6tzUqVNZtmwZH330EQBFRUVstVXqHo/jj0/97ffmm29yzDHHkJeXx8UXX1zm7ruDDjoIgIKCAo499lhuvPFG6tWrx6GHHsqnn37KH/7wB1avXs2vfvWrdW2mTJlCjx49sj7XKdNTTz3F119/vW5WaOrUqXzyySc0bdqUo446CoBp06bRv39/9t57bwYNGlRqQ3bDhg3Zd999AejatSv77LMPF110EZ999hmFhYW89957XHfddbRo0YKBAweWGvvll18u9TM029RsNg+llDQFGBYRRVnOFQK3RET13fO7niS9CBwbEZ/X0Hi3kdpE/s9c9fxQyk3DD+GhlCNHjmTSpEnMnj2bpUuXUr9+fTp06ED//v0ZOnToug3WY8eO5bTTTuOdd96hoKCA7t27M3Xq1Kx9pj+AccGCBVx99dU8//zzrFixgjZt2rDffvtx9tln8//bu+84qar7/+OvN6CARBEFkxXLWgBFFAREowERiRLjz5IYFRtYYxdjCZYEGwrRRGM0EiGAX8ESjRoUexQwCiq9iDWiIioqghoQAT+/P+buOjs7e3cXtrG8n4/HPLhz7imfu3cWPpxzZuaQQw5JjW3vvfemffv2jBpV8v9daWMX/R1bFG9Zsv8uXr58OX/84x+59957ee+999hss83o3LkzgwYNKrEEOX78eA499FCeeOIJ+vTpkxp7YWEh7733Xqny7bffvvizqq666qpSe6jy1YPMp4Jff/31jBs3joULF9KyZUv23XdfrrnmmlIzSKeffjozZsxg6tRSfz2b1RpV4kMp632yJGkgcBZwfET8pzZiy4lnb2BFRKzdlzxVfrzTi94Zl8bJ0vphQ0iW6rLRo0dzwQUX8NFHH7HJJpvUdjhcfvnljBs3jjlz5lT4gydr2jfffENBQQE33XRTic3gZrWtMslSXdqztE4iome+WaWIGBIR29eFRAkgIl6uqUQpGa/cRMnMKubEE0+kdevWdeZrOyZOnMjll19eZxMlgL/97W9stdVWpZbmzNYndWnPkplZndawYUNGjhzJ9OnTazsUAF588cXaDqFcjRs3ZvTo0aU23JutT+rNMpytGy/DrR+8DGdmVjUqswznVN8A2L11c6b6H2IzM7NS6s2eJTMzM7Pq4GTJzMzMLIWTJTMzM7MUTpbMzMzMUniDtwEw58NlFA4cX9th2FrwO+TMzKqXZ5bMzMzMUjhZMjMzM0vhZMnMzMwshZMlMzMzsxROlszMzMxS1LtkSdKPJN0n6R1Jr0l6XFLblPqdJB2S9fwqSRdXU2wDJJ2UNc6HkmYmj0NS2jWUNEPSY5Uc7z5JbdY1bjMzsw1ZvUqWJAl4GJgQETtFRHvgcuCHKc06AWUmKlUYWyPgFOCerOKbI6JT8ng8pfkFwPy1GPYO4NK1aGdmZmaJepUsAQcAqyJiWFFBRMyMiBck3S3p8KJySWMlHQZcAxyTzO4ck5xuL2mCpP9KOj+rzW8kzU0eA5KyQknzJQ2XNE/S05Ka5omtFzA9IlZX5oIkbQP8HBhRxvlGkl6V1DN5foOkwcnpF4DeSaJmZmZma6G+JUsdgGllnBsBnAwgqTmwL/A48Hvg/mR25/6k7i7AwUA3YJCkjSR1SdrvDewDnC5pz6R+G+D2iNgNWAr8Ms/4++WJ7VxJsyWNlNSijLhvITM79F2+k0ny1R+4Q9JPgT7A1cm574C3gY5l9G1mZmblqG/JUpkiYiKws6StgL7AP1NmecZHxMqI+AxYTGYZ7yfAwxHxv4j4GngI6J7UfzciZibH04DCPH0WAJ9mPb8D2InMMuBHwB9zG0g6FFgcEWUlgEXXNg+4G3gUOCUivs06vRjYOl87SWdImipp6prly9KGMDMz22DVt2RpHtAl5fzdwPFkZohGpdRbmXW8hszXwqiS9XOtAJoUPYmITyJiTTL7M5zMLFau/YDDJC0A7gN6SRpTRgy7k5nVyt2f1SQZu5SIsMo/4gAAIABJREFUuDMiukZE14abNC+jWzMzsw1bfUuWngMaSzq9qEDSXpL2T56OBgZA8WwMwFfAphXoexJwhKRNJDUDjiSzJ6ii5gM7Z8VVkHXuSGBuboOIuCwitomIQuBY4LmIOCG3nqRfAFsCPYBbJW2edbotmSTSzMzM1kK9SpYiIsgkHj9NPjpgHnAVsCg5/wmZpCV7Vul5Mhu6szd45+t7Oplk6xXgZWBERMyoRHhPkElmivxB0hxJs8lsTL8QQNLWktLeGVeCpJbAEODUiHgTuA34c3Luh8CKiPioEnGamZlZFmXyiw2DpE2AOUDniKjxTTqSHgYujYi3ami8C4EvI+Lv5dVtXNAmCvrdUgNRWVVbMOTntR2Cmdl6R9K0iOhakbr1amYpjaTewOvAX2ojUUoMJLPRu6YsBe6qwfHMzMzqnQ3m83ci4llgu1qO4Q3gjRocL20Tu5mZmVXABjOzZGZmZrY2nCyZmZmZpdhgluEs3e6tmzPVG4XNzMxK8cySmZmZWQonS2ZmZmYpnCyZmZmZpXCyZGZmZpbCG7wNgDkfLqNw4Pgq68+fKm1mZvWFZ5bMzMzMUjhZMjMzM0vhZMnMzMwshZMlMzMzsxROlszMzMxS1OtkSVJ/SVtXQ78jJLVPjhdIalmJtrdI6pEcnyvpbUmR1oekfpLeSh79Khnrs5JaVKaNmZmZfa9eJ0tAf6DKk6WIOC0iXqtsO0lbAPtExKSk6EWgN/BeOW0GAXsD3YBBlUx+7gbOrmysZmZmllEnkyVJhZJel3SXpNmSHpS0iaQDJT2cVe+nkh6S1FDSaElzJc2RdKGko4CuwFhJMyU1ldRF0kRJ0yQ9Jakg6WeCpKGSXpH0pqTuSXlDSTclfc6WdF5W/a45MTeTNF7SrCSOY/Jc2lHAk0VPImJGRCwo58dxMPBMRCyJiC+AZ4A+OWM3l/SGpHbJ83slnZ6cHgf0LWcMMzMzK0OdTJYS7YA7I2IP4EsysyPPAbtKapXUORkYBXQCWkdEh4jYHRgVEQ8CU4HjI6ITsBr4C3BURHQBRgKDs8ZrFBHdgAFkZnIAzgB2APZM4hibEm8fYFFEdIyIDmQlRVn2A6ZV6qcArYEPsp4vTMqKRcQy4FxgtKRjgRYRMTw59wXQWNKWlRzXzMzMqNvJ0gcR8WJyPAb4SUQEmWWlEyRtDvwYeAL4L7CjpL9I6kMmucrVDugAPCNpJnAlsE3W+YeSP6cBhclxb2BYRKwGiIglKfHOAXonM1TdkwQmVwHwadpF56E8ZVGqIOKZJIbbgdNyTi8mz3KkpDMkTZU0dc3yfOGamZlZXU6WchOCouejgBPILC09EBGrk9mTjsAE4BxgRJ7+BMyLiE7JY/eIOCjr/MrkzzV8/zUwyhNH/mAj3gS6kElYbpD0+zzVVgBNKtJfloXAtlnPtwEW5VaS1ADYNRlji5zTTZLy3JjvjIiuEdG14SbNKxmWmZnZhqEuJ0vbSfpxctwX+A9ARCwikyxcCYwGSN5J1iAi/gn8DuictPsK2DQ5fgNoVdSnpI0k7VZODE8DZ0pqlLTJTUKKJe+6Wx4RY4CbsmLINh/YuZwxcz0FHCSpRbKx+6CkLNeFSf99gZGSNkriEvAjYEElxzUzMzPqdrI0H+gnaTaZmZI7ss6NJbNMV/SOtNbAhGR5bTRwWVI+GhiWlDcks8F6qKRZwExg33JiGAG8D8xO2hyXUnd34JVkrCuA6/LUGQ/0LHoi6XxJC8nMFs2WNCIp71p0nCz9XQu8mjyuyV0OlNSWzNLbRRHxAjCJTDIJmdmuKUVLiWZmZlY5ymwDqlskFQKPJRul852/DZgREX+vybiqgqT/AIdGxNIaGu/PwLiI+HdavcYFbaKg3y1VNu6CIT+vsr7MzMyqmqRpEdG1/Jp1e2YpL0nTgD3IbPpeH10EbFeD480tL1EyMzOzsjUqv0rNSz57KO+sUvK2//VWRLxcw+MNr8nxzMzM6pv1bmbJzMzMrCY5WTIzMzNL4WTJzMzMLEWd3LNkNW/31s2Z6newmZmZleKZJTMzM7MUTpbMzMzMUjhZMjMzM0vhZMnMzMwshZMlMzMzsxR18rvhrOZJ+hRYCiwrp2rzlDplnSurvCXwWUVjrEFp11ib/Va2fUXrl1dvbe55Wed8z6u3fU3c87Tz/l1f937Xpq1/1ysuO9btI6JVhVpFhB9+EBEAd65LnbLOpZRPre1rXtufQ230W9n2Fa1fXr21uedlnfM9X//vedp5/66ve79r09a/69V/b7wMZ9keXcc6ZZ2rSL91SXXFu679VrZ9ReuXV29t7nllxq8LfM+r5vz6dM+hbt73tWnr3/WKW6tYvQxntUbS1IjoWttxWM3xPd8w+b5veOrbPffMktWmO2s7AKtxvucbJt/3DU+9uueeWTIzMzNL4ZklMzMzsxROlszMzMxSOFmyOklST0kvSBomqWdtx2M1Q1IzSdMkHVrbsVjNkLRr8nv+oKSzajseq36SjpA0XNK/JB1U2/FUhJMlq3KSRkpaLGluTnkfSW9IelvSwHK6CeBroAmwsLpitapRRfcc4LfAP6onSqtqVXHfI2J+RJwJHA3Um3dP1VdVdM8fiYjTgf7AMdUYbpXxBm+rcpJ6kEl0/i8iOiRlDYE3gZ+SSX5eBfoCDYEbcro4BfgsIr6T9EPgTxFxfE3Fb5VXRfd8DzKf+tuEzP1/rGait7VVFfc9IhZLOgwYCNwWEffUVPxWeVV1z5N2fwTGRsT0Ggp/rTWq7QCs/omISZIKc4q7AW9HxH8BJN0HHB4RNwBpSy5fAI2rI06rOlVxzyUdADQD2gMrJD0eEd9Va+C2Tqrqdz0ixgHjJI0HnCzVYVX0uy5gCPDE+pAogZMlqzmtgQ+yni8E9i6rsqRfAAcDmwO3VW9oVk0qdc8j4goASf1JZharNTqrLpX9Xe8J/ILMf4oer9bIrLpU6p4D5wG9geaSdo6IYdUZXFVwsmQ1RXnKylwDjoiHgIeqLxyrAZW658UVIkZXfShWgyr7uz4BmFBdwViNqOw9vxW4tfrCqXre4G01ZSGwbdbzbYBFtRSL1Qzf8w2T7/uGp97fcydLVlNeBdpI2kHSxsCxwLhajsmql+/5hsn3fcNT7++5kyWrcpLuBSYD7SQtlHRqRKwGzgWeAuYD/4iIebUZp1Ud3/MNk+/7hmdDvef+6AAzMzOzFJ5ZMjMzM0vhZMnMzMwshZMlMzMzsxROlszMzMxSOFkyMzMzS+FkyczMzCyFkyUzqxWS1kiaKWmupEclbb4OffWUtG/W8zMlnVQ1kVZo/O6S5iXX07QC9UdIar+WYxVKmlvJNi9ltT1ubcZN6fvyfGOZ1Sf+nCUDoGXLllFYWFjbYZiZmdWIadOmfRYRrSpS11+kawAUFhYyderU2g7DzMysRkh6r6J1vQxnZmZmlsLJkpmZmVkKJ0tmZmZmKZwsmZmZmaXwBm8DYM6HyygcOL62wzAzM8trwZCf19rYnlkyMzMzS+FkyczMzCyFkyUzMzOzFE6WzMzMzFI4WTIzMzNLUWeSpZwv1XxA0iYpdUt8GaSk/pJuq6a4jpD0++S4h6TpklZLOiqrTidJk5Mv0pwt6Ziscy8k1zVT0iJJjyTlxyd1Z0t6SVLHMsbfQdLLkt6SdL+kjSsR+7mSTl77qzczM7NqTZYkVeajCVZERKeI6AB8C5yZUrcQqNJvzk5xKfDX5Ph9oD9wT06d5cBJEbEb0Ae4pegb1COie3JdnYDJwENJm3eB/SNiD+Ba4M4yxh8K3BwRbYAvgFMrEftI4PxK1DczM7McqclSMoPzuqS7khmQB4tmfCQtkNQyOe4qaUJyfJWkOyU9DfyfpIaSbpT0atLHrysQ1wvAzpKulXRBVjyDJZ0PDAG6J7M1Fyant5b0ZDID84esNn0lzUlmrIZmlX+d9DdL0hRJP8xz/W2BlRHxGUBELIiI2cB32fUi4s2IeCs5XgQsBlrl9LUp0At4JKn3UkR8kZyeAmyTZ3wlbR5Miu4CjshT79as2a+DJU2S1CAilgMLJHXLbWNmZmYVU5GZpXbAnckMyJfA2RVo0wU4PCKOIzMTsiwi9gL2Ak6XtENZDZPZqJ8Bc4C/A/2S8gbAscBYYCDwQjJjc3PStBNwDLA7cIykbSVtTWZmpldyfi9JRclGM2BKRHQEJgGn5wlnP2B6Ba43O/5uwMbAOzmnjgT+HRFf5ml2KvBEnvItgaURsTp5vhBonafeQDLXfABwK3ByRBQldFOB7pW5BjMzM/teRZKlDyLixeR4DPCTCrQZFxErkuODgJMkzQReJpMAtMnTpmlSZyqZ5a6/R8QC4HNJeyb9zIiIz8sY898RsSwivgFeA7Ynk5xNiIhPk4RjLNAjqf8t8FhyPI3M0l6uAuDTClwvAJIKgLspmawU6Qvcm6fNAWSSpd/m6zJPWZQqyMwgnQ48A9wWEdmJ2mJg6zLiPUPSVElT1yxflq+KmZnZBq8ie4py/3Euer6a75OtJjl1/pd1LOC8iHiqnHFWJPt6co0gs0/oR2T24JRlZdbxGjLXli/ZKLIqIiKnfqmYgOYpfRSTtBkwHrgyIqbknNsS6EZmdim7fA8y1/ezMpLAz4DNJTVKkr1tgEVlhLA78DmlE6MmyXWUEhF3kuyValzQplQSZmZmZhWbWdpO0o+T477Af5LjBWSW2wB+mdL+KeAsSRtBZh+QpGaViPFhMpum90r6AvgK2LQCbV8G9pfUUlLDJP6JlRh7PrBzeZWSd6g9DPxfRDyQp8qvgMeSWa+iNtuR2ex9YkS8ma/fJJl7Hih6510/4F95xt8euAjYE/iZpL2zTrcF5pZ3DWZmZpZfRZKl+UA/SbOBLYA7kvKrgT9LeoHMzExZRpBZFpsuaS7wNyrxBb4R8S2ZhOEfEVE0zmxgdbI5+8KUth8BlyXtZwHTI6JUspFiErBnstEaSXtJWkgm+fmbpHlJvaPJLO/1z/qYgOxZsmMpvQT3ezJLkn9N6k8tOiHp8WS/FWSW534j6e2k/t+zO0li+ztwcbK5/FRghKSi2b79gGcrcc1mZmaWRd+vROU5KRWSmRHpUFMB5YmhAZlN1r8qesdZDY//Z+DRiFjvEo5kr9dvIuLE8uo2LmgTBf1uqYGozMzMKm/BkJ9XaX+SpkVE14rUrTMfSpmPpPbA22Q2b9d4opS4HijzAzLruJbA72o7CDMzs/VZ6nJY8m60WptViojXgB1ra/wkhk+AcbUZw9qKiGdqOwYzM7P1XZ2eWTIzMzOrbU6WzMzMzFJU5rvbrB7bvXVzplbx5jkzM7P6wDNLZmZmZimcLJmZmZmlcLJkZmZmlsLJkpmZmVkKb/A2AOZ8uIzCgeNrOwyr56r6E3jNzGqCZ5bMzMzMUjhZMjMzM0vhZMnMzMwshZMlMzMzsxROlszMzMxSOFmqBEk9JT2WHPeXdFsZ9Y6Q9PvkuIek6ZJWSzoqp14/SW8lj34p454n6Q1J8yT9oRLx7i5pdEXrm5mZWWn+6IAckhpFxOp17OZS4LDk+H2gP3BxzjhbAIOArkAA0ySNi4gvcuodABwO7BERKyVtVdEgImKOpG0kbRcR76/11ZiZmW3AanRmSVKhpPmShiezJE9LapqcmyCpa3LcUtKC5Li/pEckPSrpXUnnSvqNpBmSpiRJR1H7WyS9JGmupG6SGiSzNq2SOg0kvS2pZU5cV0m6U9LTwP9JaiJplKQ5yTgHVOIa2wIrI+IzgIhYEBGzge9yqh4MPBMRS5IE6RmgT54uzwKGRMTKpL/FecY8UtKzyiiQ9KakHyWnHwWOrWj8ZmZmVlJtLMO1AW6PiN2ApcAvK9CmA3Ac0A0YDCyPiD2BycBJWfWaRcS+wNnAyIj4DhgDHJ+c7w3MKkpkcnQBDo+I44BzACJid6AvcJekJhW8vv2A6RWo1xr4IOv5wqQsV1ugu6SXJU2UtFduhYh4GPg4iXs4MCgiPk5OTwW6VzB2MzMzy1EbydK7ETEzOZ4GFFagzfMR8VVEfAosIzNbAjAnp/29ABExCdhM0ubASL5PqE4BRpUxxriIWJEc/wS4O+nrdeA9MklLRRQAn1agnvKURZ6yRkALYB/gEuAfkvK1PQ+4jMys1r1Z5YuBrfMGIJ0haaqkqWuWL6tAyGZmZhue2kiWVmYdr+H7fVOr+T6e3Fmc7DbfZT3/jpL7rnKTjYiID4BPJPUC9gaeKCOu/2Ud50tGKmoFpePPZyGwbdbzbYBFZdR7KDJeIXPNLfPUa52c+6Gk7PvaJImplIi4MyK6RkTXhps0r0DIZmZmG5669G64BWSWwgCOSqmX5hgAST8BlkVE0XTJCDLLcf+IiDUV6GcSydJdsgdpO+CNCsYwH9i5AvWeAg6S1EJSC+CgpCzXI0CvrFg2BkosI0pqRGbG7Lhk/N9knW4LzK1g7GZmZpajLiVLNwFnSXqJ/DMnFfFF0n4YcGpW+TjgB5S9BJfrr0BDSXOA+4H+RRusK2ASsGfRUpmkvSQtBH4F/E3SPICIWAJcC7yaPK5JypA0omizO5llxB0lzQXuA/pFRO4M2uXACxHxAplE6TRJuybnDgD8DblmZmZrSaX/3V0/SZoAXBwRU/Oc6wrcHBE1stFZ0p+BRyPi2ZoYLyWOxsBE4CflfRxC44I2UdDvlpoJzDZYC4b8vLZDMDMDQNK0iOhafs26NbNULSQNBP5JZvNzTbke2KQGxyvLdsDAKvjcKDMzsw1WvflQyojoWUb5EGBIDcfyCZmlv1oVEW8Bb9V2HGZmZuuzej+zZGZmZrYunCyZmZmZpXCyZGZmZpai3uxZsnWze+vmTPU7lczMzErxzJKZmZlZCidLZmZmZimcLJmZmZml8J4lA2DOh8soHOhvRakP/CnZZmZVyzNLZmZmZimcLJmZmZmlcLJkZjXmkUceoUePHmy11VY0bdqU7bffniOOOIInn3yyuM7o0aORxIIFC2o0tvPOO4//9//+X/HzqVOncsYZZ7DLLruwySabsN1223H88cfz7rvvlmr73XffccMNN1BYWEiTJk3o2LEj//znP1PHe+mll2jQoAGSWL265Nc3rlmzhptvvpkOHTrQrFkzCgoKOPLII5k9e/ZaXdu///1vTjjhBHbaaSeaNm3KTjvtxFlnncXixYtL1f3mm2+45JJLKCgooGnTpvz4xz9m0qRJqf3fe++9SGKbbbYpdW758uUMGjSItm3b0rRpU7bddltOOumkUvf38MMP55xzzlmr6zOrboqI2o7B6oDGBW2ioN8ttR2GVYG6umfp1ltv5YILLuCUU07hiCOOoFmzZrzzzjuMHz+etm3b8oc//AGATz/9lHfeeYc999yTxo0b10hs77zzDrvuuisvvfQSXbtmvoT84osvZvLkyRx//PHstttufPjhh1x77bUsXryYmTNnsu222xa3v+KKK7jpppsYPHgwXbp04b777mP48OE89thjHHLIIaXGW7VqFZ07d+azzz7j448/ZtWqVTRq9P0W0ssvv5yhQ4dy2WWX0atXLz777DOuu+46PvzwQ2bNmpU3KUnzq1/9iq+//pqjjz6aHXfckbfeeotBgwbRuHFjZs+ezQ9+8IPiuscffzzjx4/nxhtvZMcdd+T222/niSeeYPLkyXTq1KlU30uXLmWXXXZBEg0bNmThwoUlzh933HE88sgjXH311XTt2pX333+fQYMG0bBhQ2bNmlU89owZM9h7772ZO3cubdu2rdT1ma0NSdMiomuF6jpZMnCyVJ/U1WRpu+22o0uXLjz88MOlzn333Xc0aFB7E93nnXceU6ZM4dVXXy0u+/TTT2nVqlWJeu+99x477LADV155Jddccw0AixcvZtttt2XgwIFcffXVxXUPPPBAPv3007yzQddffz333HMPhx9+ONdff32pZGnrrbemZ8+e3HPPPcVlr7/+OrvuuivDhg3j17/+daWuL9+1TJo0if3335+///3vnHLKKQDMmjWLTp06MXLkSE4++WQAVq9ezW677Ua7du0YN67094OfccYZvPfeexQUFPDss8+WSJZWrFjBpptuyqWXXsr1119fXP7kk0/ys5/9jCeffJKDDz64uLxbt2507dqVv/71r5W6PrO1UZlkyctwZZBUKGluctxV0q1V1O8tknokx6MlvStpZvIo9d82SdtLmpacnyfpzEqO96ykFlURu9m6WLJkCT/60Y/ynstOlHKX4fr374+kvI8JEyYUt5s1axaHHXYYLVq0oGnTpuy333688MIL5ca1cuVKxowZw3HHHVeiPDe5ANh+++1p1aoVH374YXHZU089xbfffssJJ5xQou4JJ5zAnDlzSi3bvfPOOwwePJi//vWvbLTRRnlj+vbbb9lss81KlG2++eZAJrGsrHzXstdeewGUuJZx48ax0UYbccwxxxSXNWrUiGOPPZannnqKlStXlujjxRdfZMyYMdx+++15x129ejVr1qyp8LUce+yxjB07lhUrVlTi6syqn5OlCoiIqRFx/rr2I2kLYJ+IyN4AcElEdEoeM/M0+wjYNyI6AXsDAyVtXYlh7wbOXvuozapGt27duOuuu7jxxht58803K9zud7/7HZMnTy7x2G+//Yr3EQFMnz6dfffdlyVLljB8+HD++c9/suWWW9K7d2+mTZuW2v+UKVNYunQp3bt3LzeW+fPns3jxYnbdddfisnnz5tG4cWN23nnnEnV32203AF577bUS5WeddRZHHXUUPXr0KHOcs88+mzFjxvCvf/2LL7/8kv/+97+cffbZbLPNNiUSmXUxceJEgFLXssMOO7DJJpuUupZvv/2Wt99+u7hs1apVnHHGGVxyySWlrr3Ipptuyoknnsitt97K888/z9dff828efO45JJL6NixIwceeGCJ+j169ODLL79k8uTJVXKNZlWlRpOlZLZmvqThySzJ05KaJucmSOqaHLeUtCA57i/pEUmPJrMw50r6jaQZkqYkCUhR+1skvSRprqRukhpIektSq6ROA0lvS2qZE9f+WbM7MyRtmnO+p6THkuMfSBolaY6k2ZJ+mZQfJGmypOmSHpD0A0o7CngyT3mZIuLbiCj671xj8twzSc0lvSGpXfL8XkmnJ6fHAX0rM6ZZdRg2bBg777wzl156Ke3ataNly5b07duXp59+OrXdTjvtxD777FP8+M9//sPkyZMZO3YsO+64IwCXXHIJ2223Hc899xxHHXUUhxxyCA8//DA77rgj1157bWr/U6ZMQRJ77LFHar3Vq1dz5pln0qpVK0499dTi8iVLlrD55psjqUT9LbbYovh8kTFjxjB16lRuvPHG1LGuueYaLrvsMn7xi1/QvHlzdtppJ+bNm8eECROK+10XX331FQMGDGDXXXfliCOOKHEtLVqUnojOdy1Dhw5l5cqVXHbZZaljjRo1iiOPPJJevXqx6aab0qFDB1atWsUzzzzDxhtvXKJux44dadCgAVOmTFmXyzOrcrUxs9QGuD0idgOWAr+sQJsOwHFAN2AwsDwi9gQmAydl1WsWEfuSmUkZGRHfAWOA45PzvYFZEfFZTv8XA+ckszfdgbQ54N8ByyJi94jYA3guSb6uBHpHRGdgKvCbPG33A3L/mzs4SbpulpR3N6ukbSXNBj4AhkbEouzzEbEMOBcYLelYoEVEDE/OfQE0lrRlyjWZVbu2bdsyY8YMJk6cyBVXXEGnTp14+OGHOfjgg7nuuusq1Mejjz7Kb3/7W4YOHVr8j/yKFSuYOHEiv/rVr2jQoAGrV69m9erVRAS9e/cu951cixYtYrPNNiv1D3euc889l5deeokxY8aUSCgiolSiVFSebcmSJVx00UVcf/31bLXVVqlj3XHHHVx33XVceeWVPP/88zzwwANsuummHHTQQSxatCi1bXlWr15N3759+fDDD7nvvvtK7JWq6LW8/fbbDB48mNtuu40mTZqkjnfllVcyZswYbrrpJiZOnMjdd9/N559/zs9+9jP+97//lai70UYb0bx583W+RrOqVhuf4P1u1nLTNKCwAm2ej4ivgK8kLQMeTcrnANn/HbwXICImSdpM0ubASOBfwC3AKcCoPP2/CPxJ0ljgoYhYmO8vjERv4NiiJxHxhaRDgfbAi0m7jckkcrkKgE+znl8GfJzUvxP4LXBNbqOI+ADYI1l+e0TSgxHxSU6dZyT9Crgd6JjTxWJga+Dz7EJJZwBnADTcrPSeBrOq1rBhQ3r06FG8BLVo0SL69OnD1VdfzTnnnJN3VqPIrFmzOO644zj11FO5+OKLi8uXLFnCmjVruPbaa8ucRUrbQP7NN9+U+667yy67jDvvvJO77rqLgw46qMS5LbbYgi+++KJUovHFF18Un4dM0vDDH/6Qo48+mqVLlxaPDbBs2TKaNGlCs2bNWLJkCRdeeCGXXHJJiQ3jvXr1orCwkBtvvJGbb745Nd6yfPfdd/Tr149nn32W8ePHl5pN22KLLXj//fdLtcu9lvPPP59evXqxzz77FF/Lt99+S0SwdOlSGjduTNOmTZk3bx5DhgxhxIgRJWbj9t57b9q2bcuIESO44IILSozVtGlT71myOqc2kqXsHYJrgKbJ8Wq+n+nK/a9Kdpvvsp5/R8lryH1rX0TEB5I+kdSLzJ6f4yldaYik8cAhwBRJvYFvyohfecYR8ExElLfctYKsa4uIj5LDlZJGkZnhKlNELJI0j8zs14MlApAaALsmY2wBZL9/twl5Zssi4k4ySRqNC9r4bZFW47beemtOO+00LrjgAt566y26deuWt94nn3zCYYcdxj777FPqnVKbb745DRo04JxzzuGkk07K2z7tnXZbbrllcTKQz+DBgxkyZAi33norJ554Yqnzu+22GytXruSdd94psXenaK9S+/bti5/PmTOHLbcsPcnbsmVLDj/8cB555BHefPNNVq5cWbx1nlOoAAAWZUlEQVQBu8gWW2zBTjvtxPz588uMtTxnnnkm999/Pw8++GCp/UJF1/Lwww+zfPnyEvuWXnvtNTbeeOPi63vttdd477338ia3LVq04IILLuCWW25hzpw5AKWupU2bNmy++eZ5r2XJkiW0bNmyVLlZbapLG7wXAF2S46PWso9jACT9hMxS2bKkfASZ5bh/RMSa3EaSdoqIORExlMwS2i4pYzxNZsmrqG0LYAqwn6Sdk7JNJOX7oJD5wM5ZbQuSPwUcAczNE9s2Wfu6WpBZynsjT98XJv33BUZK2iir7x+R+fma1ZoPPvggb/nrr78OUOY75b755hsOP/xwmjVrxgMPPFBi2QigWbNmdO/enVmzZtG5c2e6du1a6pFml112YdWqVaU+Hwgynw115ZVXMnjwYM4777y87fv06cPGG2/M2LFjS5SPGTOGDh06sMMOOwBwyy238Pzzz5d49OvXD4Bnn322eCmy6OfwyiuvlOhvyZIlvP3227Ru3Tr1espy0UUXMWLECEaNGlVin1K2ww47jFWrVvHAAw8Ul61evZr777+fgw46qHgG7r777it1LQcffDAtW7bk+eef59xzz029ljfffJOlS5eWupaPP/6Yb775hnbt2q3VNZpVl7r0Rbo3Af+QdCLw3Fr28YWkl4DNyCy5FRlHZvkt3xIcwABJB5CZ6XoNeILMklk+1wG3Jx8rsAa4OiIektQfuDdr39GVQO5bfsYDvyaTvAGMTTafC5gJnAmZjyoAzoyI08jMFv1RUiT1boqIOdmdJonZaUC3iPhK0qRk/EFkEtApEVHyI4LNaliHDh044IADOPLII9lhhx348ssvefzxxxk2bBhHH3108Tvbcg0YMIDp06czevTo4sSqSPv27dlss83405/+RI8ePTj44IM59dRTKSgo4LPPPmP69OmsWbOGIUOGlBlX0ZLgK6+8UuLDHu+77z4GDBhAnz596NWrV4lNx5tttlnxjNFWW23FhRdeyA033MCmm25K586duf/++3nuuef417/+Vdwm3wc6Fn30wf7771+cBBYWFnLooYdy44030qBBA/bff38+//xz/vCHP7By5UrOOuusEu0POOAARo0aRf/+/cu8xqFDh/KnP/2JU045hTZt2pS4llatWrHTTjsVx3jMMccwYMAAVq1axQ477MAdd9zBu+++WyIZ3GeffUqNMXr0aBo3bkzPnj2Ly7p3707Hjh256KKL+OKLL4o/lPK6666jefPmxclikZdffhkg9Z2CZrWh3nwopaQJwMURMTXPua7AzRFR/nuDq5mk/wCHRsTSGhrvz8C4iPh3Wj1/KGX9UVc/lHLYsGE8/vjjzJo1i08++YSGDRvStm1b+vbty4ABA4o3WI8ePZqTTz6Zd999l8LCQnr27Fn8Nvdczz//fPE/zvPnz+fqq6/mueeeY9myZbRq1YrOnTtz5pln5v0U7Wx777037du3Z9So7/8/1b9/f+6666689ffff/8Sn/G0Zs0abrjhBoYPH87HH39Mu3bt+P3vf89RR6VPkl911VVcffXVpT6Ucvny5fzxj3/k3nvv5b333mOzzTajc+fODBo0qMRS5fjx4zn00EN54okn6NOnT5njpP0M+/Xrx+jRo4ufr1ixgiuuuIJ77rmHpUuX0rFjR4YOHVoiCcqnf//+pT6UEuDzzz/n+uuvZ9y4cSxcuJCWLVuy7777cs0115SaQTr99NOZMWMGU6eW+mvcrMppQ/wE77KSJUkDgbOA4yPiP7URW048ewMrImLtvuSp8uOdXvTOuDROluqPupos1WWjR4/mggsu4KOPPir1GUN12eWXX864ceOYM2dO3nexrU+++eYbCgoKuOmmm0psBjerLpVJlurSnqV1EhE9880qRcSQiNi+LiRKABHxck0lSsl45SZKZhu6E088kdatW693X7MxceJELr/88vU+UQL429/+xlZbbVVqac6sLqhLe5bMzGpFw4YNGTlyJNOnT6/tUCrlxRdfrO0Qqkzjxo0ZPXp0qQ38ZnVBvVmGs3XjZbj6w8twZmbl2yCX4czMzMyqg+c7DYDdWzdnqmckzMzMSvHMkpmZmVkKJ0tmZmZmKZwsmZmZmaVwsmRmZmaWwhu8DYA5Hy6jcOD42g5jg+W3+5uZ1V2eWTIzMzNL4WTJzMzMLIWTJTMzM7MUTpbMzMzMUtTLZEnSjyTdJ+kdSa9JelxS25T6nSQdkvX8KkkXV1NsAySdlBxfK2m2pJmSnpa0dRlt/iBpnqT5km5VJb5iXNKzklpUVfxmZmYbmnqXLCWJxMPAhIjYKSLaA5cDP0xp1gk4JOV8VcXWCDgFuCcpujEi9oiITsBjwO/ztNkX2A/YA+gA7AXsX4lh7wbOXpe4zczMNmT1LlkCDgBWRcSwooKImBkRL0i6W9LhReWSxko6DLgGOCaZ4TkmOd1e0gRJ/5V0flab30iamzwGJGWFyazP8GQG6GlJTfPE1guYHhGrk7i+zDrXDIg8bQJoAmwMNAY2Aj7JriCpuaQ3JLVLnt8r6fTk9Digb3k/NDMzM8uvPiZLHYBpZZwbAZwMmQQD2Bd4nMyMzv0R0Ski7k/q7gIcDHQDBknaSFKXpP3ewD7A6ZL2TOq3AW6PiN2ApcAv84y/X25skgZL+gA4njwzSxExGXge+Ch5PBUR83PqLAPOBUZLOhZoERHDk3NfAI0lbVnGz8TMzMxS1MdkqUwRMRHYWdJWZGZb/lk0y5PH+IhYGRGfAYvJLOP9BHg4Iv4XEV8DDwHdk/rvRsTM5HgaUJinzwLg05yYroiIbYGxZBKeEiTtDOwKbAO0BnpJ6pHn2p4B5gC3A6flnF4MlNoPJekMSVMlTV2zfFmecM3MzKw+JkvzgC4p5+8mM4tzMjAqpd7KrOM1ZD7tPG1jdb76uVaQWVLL5x7yz0YdCUyJiK+TBO0JMrNaJUhqQCapWgFskXO6SVJeQkTcGRFdI6Jrw02alxGWmZnZhq0+JkvPkVl2Ktqzg6S9JBVtih4NDACIiHlJ2VfAphXoexJwhKRNJDUjk8i8UInY5gM7Z8XVJuvcYcDredq8D+wvqZGkjchs7p6fp96FSXlfYGRSt2jD+4+ABZWI08zMzBL1LlmKiCCTxPw0+eiAecBVwKLk/CdkkorsWaXnyWzozt7gna/v6WSSrVeAl4ERETGjEuE9AWQvoQ1JNorPBg4CLgCQ1FXSiKTOg8A7ZJbYZgGzIuLR7E6Tj0U4DbgoIl4gk9RdmZzuQmZmqqzlRjMzM0uhTG6x4ZC0CZnEo3OyMbqmx38YuDQi3qqh8f4MjIuIf6fVa1zQJgr63VITIVke/iJdM7OaJWlaRHStSN16N7OURlJvMktdf6mNRCkxkMxG75oyt7xEyczMzMqWbxNyvRURzwLb1XIMbwBv1OB4w2tqLDMzs/pog5pZMjMzM6ssJ0tmZmZmKZwsmZmZmaXYoPYsWdl2b92cqX5HlpmZWSmeWTIzMzNL4WTJzMzMLIWTJTMzM7MUTpbMzMzMUniDtwEw58NlFA4cX2X9+es7zMysvvDMkpmZmVkKJ0tmZmZmKZwsmZmZmaVwsmRmZmaWwsmSmZmZWYp6nSxJ6i9p62rod4Sk9snxAkktK9H2Fkk9kuNzJb0tKdL6kNRP0lvJo18lY31WUovKtDEzM7Pv1etkCegPVHmyFBGnRcRrlW0naQtgn4iYlBS9CPQG3iunzSBgb6AbMKiSyc/dwNmVjdXMzMwy6mSyJKlQ0uuS7pI0W9KDkjaRdKCkh7Pq/VTSQ5IaShotaa6kOZIulHQU0BUYK2mmpKaSukiaKGmapKckFST9TJA0VNIrkt6U1D0pbyjppqTP2ZLOy6rfNSfmZpLGS5qVxHFMnks7Cniy6ElEzIiIBeX8OA4GnomIJRHxBfAM0Cdn7OaS3pDULnl+r6TTk9PjgL7ljGFmZmZlqJPJUqIdcGdE7AF8SWZ25DlgV0mtkjonA6OATkDriOgQEbsDoyLiQWAqcHxEdAJWA38BjoqILsBIYHDWeI0iohswgMxMDsAZwA7AnkkcY1Pi7QMsioiOEdGBrKQoy37AtEr9FKA18EHW84VJWbGIWAacC4yWdCzQIiKGJ+e+ABpL2jK3Y0lnSJoqaeqa5csqGZaZmdmGoS4nSx9ExIvJ8RjgJxERZJaVTpC0OfBj4Angv8COkv4iqQ+Z5CpXO6AD8IykmcCVwDZZ5x9K/pwGFCbHvYFhEbEaICKWpMQ7B+idzFB1TxKYXAXAp2kXnYfylEWpgohnkhhuB07LOb2YPMuREXFnRHSNiK4NN2leybDMzMw2DHU5WcpNCIqejwJOILO09EBErE5mTzoCE4BzgBF5+hMwLyI6JY/dI+KgrPMrkz/X8P3XwChPHPmDjXgT6EImYblB0u/zVFsBNKlIf1kWAttmPd8GWJRbSVIDYNdkjC1yTjdJys3MzKyS6nKytJ2kHyfHfYH/AETEIjLJwpXAaIDknWQNIuKfwO+Azkm7r4BNk+M3gFZFfUraSNJu5cTwNHCmpEZJm9wkpFjyrrvlETEGuCkrhmzzgZ3LGTPXU8BBklokG7sPSspyXZj03xcYKWmjJC4BPwIWVHJcMzMzo24nS/OBfpJmk5kpuSPr3Fgyy3RF70hrDUxIltdGA5cl5aOBYUl5QzIbrIdKmgXMBPYtJ4YRwPvA7KTNcSl1dwdeSca6ArguT53xQM+iJ5LOl7SQzGzRbEkjkvKuRcfJ0t+1wKvJ45rc5UBJbcksvV0UES8Ak8gkk5CZ7ZpStJRoZmZmlaPMNqC6RVIh8FiyUTrf+duAGRHx95qMqypI+g9waEQsraHx/gyMi4h/p9VrXNAmCvrdUmXjLhjy8yrry8zMrKpJmhYRXcuvWbdnlvKSNA3Yg8ym7/XRRcB2NTje3PISJTMzMytbo/Kr1Lzks4fyziolb/tfb0XEyzU83vCaHM/MzKy+We9mlszMzMxqkpMlMzMzsxR1chnOat7urZsz1ZuyzczMSvHMkpmZmVkKJ0tmZmZmKZwsmZmZmaVwsmRmZmaWwsmSmZmZWYo6+XUnVvMkfQq8V0XdNQeW1XD7irapSL3y6pR1vqzylsBnFYittqzr/arOvv1aqFl+LaxdnbU5V5dfC3X5dbA2fZRVf/uIaFWhHiLCDz+q9AHcWdPtK9qmIvXKq1PW+ZTyqbV9T6rzfvm14NfChv5aWJtzdfm1UJdfB2vTR1WM6WU4qw6P1kL7irapSL3y6pR1fl2vu7ZUZ9x+Laxf/FpYuzpre66uqsuvg7XpY53H9DKcWTWTNDUq+M3WVr/5tWBF/FpYv3hmyaz63VnbAVid4deCFfFrYT3imSUzMzOzFJ5ZMjMzM0vhZMnMzMwshZMlMzMzsxROlsxqkaQjJA2X9C9JB9V2PFZ7JO0o6e+SHqztWKxmSWom6a7k74LjazseK83JktlakjRS0mJJc3PK+0h6Q9Lbkgam9RERj0TE6UB/4JhqDNeqURW9Fv4bEadWb6RWUyr5mvgF8GDyd8FhNR6slcvJktnaGw30yS6Q1BC4HfgZ0B7oK6m9pN0lPZbz2Cqr6ZVJO1s/jabqXgtWP4ymgq8JYBvgg6TamhqM0SqoUW0HYLa+iohJkgpzirsBb0fEfwEk3QccHhE3AIfm9iFJwBDgiYiYXr0RW3WpiteC1S+VeU0AC8kkTDPxJEad5JtiVrVa8/3/ECHzl2DrlPrnAb2BoySdWZ2BWY2r1GtB0paShgF7SrqsuoOzWlHWa+Ih4JeS7mD9/HqUes8zS2ZVS3nKyvzk14i4Fbi1+sKxWlTZ18LngBPm+i3vayIi/gecXNPBWMV5Zsmsai0Ets16vg2wqJZisdrl14Ll8mtiPeVkyaxqvQq0kbSDpI2BY4FxtRyT1Q6/FiyXXxPrKSdLZmtJ0r3AZKCdpIWSTo2I1cC5wFPAfOAfETGvNuO06ufXguXya6J+8RfpmpmZmaXwzJKZmZlZCidLZmZmZimcLJmZmZmlcLJkZmZmlsLJkpmZmVkKJ0tmZmZmKZwsmVmtkLRG0kxJcyU9Kmnzdeirp6R9s56fKemkqom0QuN3lzQvuZ6mFag/Ivm2+bUZq1DS3Eq2eSmr7XFrM25K35fnG8usPvHnLJlZrZD0dUT8IDm+C3gzIgavZV9XAV9HxE1VGGJlxh8GvBwRo2pgrELgsYjosBZtewIXR8ShlWjTMCLWpJwvvo9m9ZVnlsysLphM5tvXi2aJHis6Iek2Sf2T4wWSrpY0XdIcSbskycOZwIXJzE53SVdJujhpM0HSzZImSZovaS9JD0l6S9J1WeOcIOmVpI+/SWqYG6SkAyXNSMYeKamxpNOAo4HfSxqbU7+ZpPGSZiUzaMdkxdQ1Of5a0uCkzhRJP0zKd0qevyrpGklf54mnoaQbkzqzJf063w83q+0QoHtyjReW1T65B89LugeYk5Q9ImlaMoN2RlI2BGia9Dc2eyxl3Jhc95ysa++ZXP+Dkl6XNFZSvi+YNasznCyZWa1KkpIDqfh3ZH0WEZ2BO8jMkiwAhgE3R0SniHghT5tvI6JHUu9fwDlAB6C/pC0l7QocA+wXEZ2ANcDxOXE2AUYDx0TE7kAj4KyIGJHEfklElGgD9AEWRUTHZCboyTyxNQOmRERHYBJwelL+Z+DPEbEXZX/Z6qnAsqTOXsDpknYooy7AQOCF5Od0czntuwFXRETRcuEpEdEF6AqcL2nLiBgIrEj6y732XwCdgI5Ab+BGSQXJuT2BAUB7YEdgv5SYzWqdkyUzqy1NJc0EPge2AJ6pYLuHkj+nAYUVbFOUiM0B5kXERxGxEvgvmW+BPxDoAryaxHQgmX/Es7UD3o2IN5PndwE9yhl3DtBb0lBJ3SNiWZ463wJFM2nZ1/Rj4IHk+J4y+j8IOCmJ+WVgS6BNOTFVtP0rEfFuVt3zJc0CppD5mZU3zk+AeyNiTUR8Akwkk5AV9b0wIr4DZlLx+2hWKxrVdgBmtsFaERGdJDUnkyycA9wKrKbkf+Sa5LRbmfy5hor/HVbU5rus46LnjQABd0XEZSl9VHqpKCLelNQFOAS4QdLTEXFNTrVV8f3m0cpcU1FM50XEU5WNLa19srfpfznPewM/jojlkiZQ+r7k67ss2fegstdsVuM8s2RmtSqZbTkfuFjSRsB7QPtkP1BzMrM85fkK2HQdwvg3cJSkrQAkbSFp+5w6rwOFknZOnp9IZrakTJK2BpZHxBjgJqBzJWKaAvwyOT62jDpPAWclPzcktZXULKXP3J9TRds3B75IEqVdgH2yzq0qap9jEnBMsi+qFZlZuFdSYjOrs5wsmVmti4gZwCzg2Ij4APgHMBsYC8yoQBePAkcWbfBei/FfA64EnpY0m8ySYEFOnW+Ak4EHJM0hMys1rJyudwdeSZa5rgCuK6d+tgHAbyS9ksSSbwlvBPAaMF2ZjxP4G+mzNLOB1clm8gsr0f5JoFHys7mWTCJX5E5gdu7mduDhZLxZwHPApRHxcdoFm9VV/ugAM7M6SNImZJYqQ9KxQN+IOLy24zLbEHmd2MysbuoC3Ja8rX4pcEotx2O2wfLMkpmZmVkK71kyMzMzS+FkyczMzCyFkyUzMzOzFE6WzMzMzFI4WTIzMzNL4WTJzMzMLMX/BwN7hYEmcbWLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_order = ['pystencils', 'Cython', 'numba', 'numpy slice', 'numpy roll', 'pure Python']\n",
    "plot_order = [p for p in plot_order if p in all_implementations]\n",
    "\n",
    "def bar_plot(*shape):\n",
    "    names = plot_order\n",
    "    runtimes = tuple(do_benchmark(all_implementations[name], shape) for name in names)\n",
    "    speedups = tuple(runtime / min(runtimes) for runtime in runtimes)\n",
    "    y_pos = np.arange(len(names))\n",
    "    labels = tuple(f\"{name} ({round(speedup, 1)} x)\" for name, speedup in zip(names, speedups))\n",
    "    \n",
    "    plt.text(0.5, 0.5, f\"Size {shape}\", horizontalalignment='center', fontsize=16,\n",
    "             verticalalignment='center', transform=plt.gca().transAxes)\n",
    "    plt.barh(y_pos, runtimes, log=True)\n",
    "     \n",
    "    plt.yticks(y_pos, labels);\n",
    "    plt.xlabel('Runtime of single iteration')\n",
    "    \n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "plt.subplot(3, 1, 1)\n",
    "bar_plot(32, 32)\n",
    "\n",
    "plt.subplot(3, 1, 2)\n",
    "bar_plot(128, 128)\n",
    "\n",
    "plt.subplot(3, 1, 3)\n",
    "bar_plot(2048, 2048)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All runtimes are plotted logarithmically. Numbers next to the labels show how much slower the version is than the fastest one."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}