demo_benchmark.ipynb

{
 "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 benchmark we compare different ways of implementing a simple stencil kernel in Python.\n",
    "The benchmark kernel computes the average of the four neighbors in 2D and stores 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(input_arr, output_arr):       \n",
    "    for x in range(1, input_arr.shape[0] - 1):\n",
    "        for y in range(1, input_arr.shape[1] - 1):\n",
    "            output_arr[x, y] = (input_arr[x + 1, y] + input_arr[x - 1, y] +\n",
    "                                input_arr[x, y + 1] + input_arr[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(input_arr, output_arr):\n",
    "    neighbors = [np.roll(input_arr, axis=a, shift=s) for a in (0, 1) for s in (-1, 1)]\n",
    "    np.divide(sum(neighbors), 4, out=output_arr)"
   ]
  },
  {
   "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(input_arr, output_arr):\n",
    "    output_arr[1:-1, 1:-1] = input_arr[2:, 1:-1] + input_arr[:-2, 1:-1] + \\\n",
    "                             input_arr[1:-1, 2:] + input_arr[1:-1, :-2]\n",
    "    "
   ]
  },
  {
   "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] input_arr, object[double, ndim=2] output_arr):\n",
    "    cdef int xs, ys, x, y\n",
    "    xs, ys = input_arr.shape\n",
    "    for x in range(1, xs - 1):\n",
    "        for y in range(1, ys - 1):\n",
    "            output_arr[x, y] = (input_arr[x + 1, y] + input_arr[x - 1, y] +\n",
    "                                input_arr[x, y + 1] + input_arr[x, y - 1]) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we also create a *pystencils* version of the same stencil code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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",
    "kernel = ps.create_kernel(update).compile()\n",
    "\n",
    "def avg_pystencils(input_arr, output_arr):\n",
    "    kernel(src=input_arr, dst=output_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_implementations = {\n",
    "    'pure Python': avg_pure_python,\n",
    "    'numpy roll': avg_numpy_roll,\n",
    "    'numpy slice': avg_numpy_slice,\n",
    "    'Cython': None,\n",
    "    'pystencils': avg_pystencils,\n",
    "}\n",
    "if 'avg_cython' in globals():\n",
    "    all_implementations['Cython'] = avg_cython\n",
    "else:\n",
    "    del all_implementations['Cython']"
   ]
  },
  {
   "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": 29,
   "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",
    "    timer = timeit.Timer('f(a, b)', globals={'f': func, 'a': in_arr, 'b': out_arr})\n",
    "    calls, time_taken = timer.autorange()\n",
    "    return time_taken / calls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHnCAYAAABZgEKHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmclWX9//HXG1A2BRRIEcVxARVBEUcwDUJFMTWBoq+BPwUVy1yCEvcVxcK0chcVZUzQyhIDQUNj0wJ12CFcExMxRVlcWBT4/P449+CZmTMbzsJw3s/H4zy4z3Vfy+cej/CZ67ru+ygiMDMzM8tmdWo6ADMzM7Oa5oTIzMzMsp4TIjMzM8t6TojMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrFevpgOw6tOiRYvIycmp6TDMzMyqzZw5cz6OiJZl1XNClEVycnLIz8+v6TDMzMyqjaR3y1PPS2ZmZmaW9ZwQmZmZWdZzQmRmZmZZzwmRmZmZZT1vqs4ii95fS86Vkyq1z2UjT63U/szMzGqCZ4jMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzr1fqESNIgSXtVQb+jJbVPjpdJalGBtndI6p4cXyzpLUlRWh+SBkp6M3kNrGCsL0jarSJtzMzM7Gu1PiECBgGVnhBFxOCI+HdF20naHTg6ImYmRf8EegIlfpdK0uYGoCvQBbihggnOY8CFFY3VzMzMUmosIZKUI+k1SY9KWijpL5IaSTpB0vi0eidKekpSXUl5khZLWiTpF5L6AbnAOEnzJTWUdKSkGZLmSPq7pFZJP9Ml3SrpFUlvSOqWlNeVdHvS50JJl6TVzy0Sc2NJkyQtSOI4I8Ol9QOeK3gTEfMiYlkZP45ewPMRsSoiVgPPAycXGbuppNclHZS8f0LS+cnpCUD/MsYwMzOzEtT0DNFBwIMRcRjwKalZjqnAIZJaJnXOAcYAnYDWEdEhIjoCYyLiL0A+cGZEdAI2AXcD/SLiSOAR4Ja08epFRBdgKKkZGYCfAPsBRyRxjCsl3pOBFRFxeER0IC3xSXMsMKdCPwVoDbyX9n55UrZVRKwFLgbyJP0Y2C0iHkrOrQbqS2pewXHNzMyMmk+I3ouIfybHY4HvRESQWgL6f5KaAd8GngX+A+wv6W5JJ5NKoIo6COgAPC9pPnAtsHfa+aeSP+cAOclxT2BURGwCiIhVpcS7COiZzDR1S5KUoloBK0u76AyUoSyKFUQ8n8RwLzC4yOmPyLB0KOknkvIl5W9elylcMzMzq+mEqOg/+gXvxwD/j9Qy0JMRsSmZBTkcmA5cBIzO0J+AJRHRKXl1jIiT0s5vTP7czNdfW6IMcWQONuIN4EhSScmvJV2fodp6oEF5+kuzHNgn7f3ewIqilSTVAQ5Jxti9yOkGSXnRmB+MiNyIyK3bqGkFwzIzM8sONZ0QtZH07eS4P/ASQESsIJUQXAvkASR3aNWJiL8C1wGdk3afAbsmx68DLQv6lLSTpEPLiGEKcIGkekmboonGVsndbOsiYixwe1oM6ZYCB5YxZlF/B06StFuymfqkpKyoXyT99wcekbRTEpeAPYFlFRzXzMzMqPmEaCkwUNJCUjMe96edG0dqSa3gTq/WwPRkKSwPuCopzwNGJeV1SW1qvlXSAmA+cEwZMYwG/gssTNoMKKVuR+CVZKxrgBEZ6kwCehS8kfRzSctJzfoslDQ6Kc8tOE6W6W4GXk1eNxVdupPUjtQy2aUR8SIwk1TCCKlZq9kFy35mZmZWMUpt2amBgaUc4Jlkc3Km8/cA8yLi4eqMqzJIegk4LSLWVNN4dwITIuIfpdWr36pttBp4R6WO7W+7NzOz7ZmkORGRW1a9mp4hykjSHOAwUhuta6NLgTbVON7ispIhMzMzK1m9sqtUjeTZPBlnh5Jb5mutiHi5msd7qDrHMzMz29FslzNEZmZmZtXJCZGZmZllvRpbMrPq17F1U/K9CdrMzKwYzxCZmZlZ1nNCZGZmZlnPCZGZmZllPSdEZmZmlvW8qTqLLHp/LTlXTqqx8f1UazMz2155hsjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrOeEyMzMzLJerU2IJO0p6Y+S3pb0b0mTJbUrpX4nSaekvb9R0rAqim2opLOT4x9JWiJpi6TcUto8IukjSYu3YbzbJR3/TWI2MzPLZrUyIZIkYDwwPSIOiIj2wNXAHqU06wScUsr5yoqtHnAu8HhStBj4ATCzjKZ5wMnbOOzdwJXb2NbMzCzr1cqECDgO+CoiRhUURMT8iHhR0mOSeheUSxon6XTgJuAMSfMlnZGcbi9puqT/SPp5WptfSlqcvIYmZTmSlkp6KJnxmSKpYYbYjgfmRsSmJK6lEfF6WRcUETOBVaXVkfS3tJmnn0oal7R9F2guac+yxjEzM7PiamtC1AGYU8K50cA5AJKaAscAk4HrgT9FRKeI+FNS92CgF9AFuEHSTpKOTNp3BY4Gzpd0RFK/LXBvRBwKrAF+mGH8Y0uJ7Zv6CXC9pG7ApcAlaefmJmObmZlZBdXWhKhEETEDOFDSt4D+wF8LZmsymBQRGyPiY+AjUktu3wHGR8QXEfE58BTQLan/TkTMT47nADkZ+mwFrKycqyksIj4kldhNAy6NiPQZpY+AvYq2kfQTSfmS8jevW1sVYZmZmdV6tTUhWgIcWcr5x4AzSc30jCml3sa0482kvspEFaxf1HqgQSl9fFMdgU8onvw0SMYuJCIejIjciMit26hpFYZlZmZWe9XWhGgqUF/S+QUFko6S9N3kbR4wFCAiliRlnwG7lqPvmUAfSY0kNQb6Ai9WILalwIEVqF9ukroA3wOOAIZJ2i/tdDtSG7jNzMysgmplQhQRQSpROTG57X4JcCOwIjn/IanEJH12aBqpTdTpm6oz9T2XVEL1CvAyMDoi5lUgvGeB7gVvJPWVtBz4NjBJ0t+T8r0kTU6r9wQwCzhI0nJJ56V3Kqk+8BBwbkSsILWH6BGl7EQqCcuvQJxmZmaWUCq32LFIagQsAjpHRLVvnJE0Hrg8It6spvH6krrW60qrV79V22g18I7qCCkjf9u9mZlVN0lzIqLE5wAWqJUzRKWR1BN4Dbi7JpKhxJWkNldXl3rAb6txPDMzsx1Kpk3BtVpEvAC0qeEYXgfKfPZQJY73ZHWNZWZmtiPa4WaIzMzMzCrKCZGZmZllPSdEZmZmlvV2uD1EVrKOrZuS7zu9zMzMivEMkZmZmWU9J0RmZmaW9ZwQmZmZWdbzHqIssuj9teRcOammwzDzU8vNbLvjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzLZ7Tz/9NN27d+db3/oWDRs2ZN9996VPnz4899xzW+vk5eUhiWXLllVrbJdccgnf//73C5VdffXVnHTSSTRv3hxJ5OXlldh+9erVDB06lDZt2lC/fn323ntvBg0atE2xLF68mJ/+9KcceeSR7Lzzzkgqtf7s2bM5+eSTadasGY0bN6Zjx4788Y9/3Hp+3rx5NGrUiP/+97/bFI9ZbeKEyMy2a3fddRd9+/albdu2PPzww0yaNIlrr70WgKlTp26td+qppzJr1ixatWpVbbG9/fbbPPDAA9xwww2Fyu+++27Wr1/PaaedVmr71atX853vfIcXXniBESNG8Pzzz3P77bez6667blM8c+bMYfLkybRp04bc3NxS606aNInu3buz55578vjjj/O3v/2N888/nw0bNmytc8QRR3DiiSdy3XXXbVM8ZrWJIqKmY7BqUr9V22g18I6aDsOsQneZtWnThiOPPJLx48cXO7dlyxbq1Km53+suueQSZs+ezauvvlqovCCut956i7Zt2zJmzJiMsz4XXHABzz77LIsWLaJJkybfOJ70n8e1117LLbfcQqa/4z/77DMOOOAABgwYwB13lP53wuTJk+nduzfvvvsue+211zeO0ay6SZoTEaX/hkCWzxBJypG0ODnOlXRXJfV7h6TuyXGepHckzU9enUpo00bSFElLJf1bUk4FxvujpLaVEbvZ9mbVqlXsueeeGc+lJ0NFl8wGDRqEpIyv6dOnb223YMECTj/9dHbbbTcaNmzIsccey4svvlhmXBs3bmTs2LEMGDCg1LhK8sUXX/CHP/yBwYMHV0oyVN5xAZ588klWrlzJpZdeWmbdk046iSZNmpS67Ge2I8jqhChdRORHxM+/aT+SdgeOjoiZacWXRUSn5DW/hKZ/AG6LiEOALsBHFRj2fuDybYvYbPvWpUsXHn30UW677TbeeOONcre77rrrmDVrVqHXscceS6NGjWjTpg0Ac+fO5ZhjjmHVqlU89NBD/PWvf6V58+b07NmTOXPmlNr/7NmzWbNmDd26ddum65ozZw7r169njz32oF+/fjRs2JBddtmFPn368M4772xTn+X10ksvsfvuu7No0SI6duxIvXr12GeffRg+fDibN28uVLdevXp8+9vfLrRfy2xHVOkJUTLrslTSQ5KWJLMeDZNz0yXlJsctJC1LjgdJelrSxGQ25WJJv5Q0T9LsJMkoaH+HpH9JWiypi6Q6kt6U1DKpU0fSW5JaFInru2mzNPMk7VrkfA9JzyTHu0gaI2mRpIWSfpiUnyRplqS5kp6UtEuGH0E/oEJ/c0hqD9SLiOcBIuLziFhXpE49Sa9K6pG8/7WkW5LTLwI9JflBm7bDGTVqFAceeCCXX345Bx10EC1atKB///5MmTKl1HYHHHAARx999NbXSy+9xKxZsxg3bhz7778/AJdddhlt2rRh6tSp9OvXj1NOOYXx48ez//77c/PNN5fa/+zZs5HEYYcdtk3XtWLFCgCGDRtG3bp1mTBhAg8++CDz5s2jR48efPbZZ9vUb3nHXrduHQMGDGDQoEG88MILDBw4kJtvvplhw4YVq3/EEUfwyiuvsGXLliqLyaymVdUMUVvg3og4FFgD/LAcbToAA0jNjtwCrIuII4BZwNlp9RpHxDHAhcAjEbEFGAucmZzvCSyIiI+L9D8MuCgiOgHdgPWlxHIdsDYiOkbEYcDUJMG6FugZEZ2BfOCXGdoeCxT91fKWJLH6vaT6Gdq0A9ZIeipJ1m6TVDe9QkRsAgYB90s6ETgZGJ6c2wK8BRxetGNJP5GULyl/87q1pVyy2fapXbt2zJs3jxkzZnDNNdfQqVMnxo8fT69evRgxYkS5+pg4cSJXXHEFt956K3369AFg/fr1zJgxgx/96EfUqVOHTZs2sWnTJiKCnj17MnPmzFL7XLFiBU2aNGHnnXfepusqSC72228//vjHP3LiiScyYMAA/vznP/Pf//6XsWPHblO/5R17w4YNXH/99Vx66aX06NGDESNGcP7553Pvvfeydm3hvytatmzJxo0bWbVqVZXFZFbTqioheidtaWgOkFOONtMi4rOIWAmsBSYm5YuKtH8CIFmSaiKpGfAIXydN5wJjMvT/T+B3kn4ONEsSjJL0BO4teBMRq4GjgfbAPyXNBwYC+2Zo2wpYmfb+KuBg4Chgd+CKDG3qkUrShiX19ieV/BQSEUuAx0j9bM6NiC/TTn8EFNvxGBEPRkRuROTWbdS0hMs1277VrVuX7t27M2LECF544QX+85//0LFjR4YPH87q1atLbbtgwQIGDBjAeeedV2j2Y9WqVWzevJmbb76ZnXbaqdDrnnvuYfXq1aXOiGzYsIH69TP9flM+zZs3B6Bnz56Fbo/v2rUrTZo0Yd68edvcd3nHPvHEEwuVn3TSSXz11VcsWbKkUHnDhg2BVBJptqOqqiWWjWnHm4GGyfEmvk7CGpTSZkva+y0UjrPoLRMREe9J+lDS8UBXvp4tSq80UtIk4BRgtqSewIai9RLKMI6A5yOifwltCqwn7doi4oPkcKOkMaSSnqKWA/Mi4j8Akp4mlYA9nKFuR1KzbnsUKW9A6bNeZjuMvfbai8GDBzNkyBDefPNNunTpkrHehx9+yOmnn87RRx/NfffdV+hcs2bNqFOnDhdddBFnn312xvalbVJu3rx5mclYaQ499FCAEp8VVJV3z5U0dsEdaUXHLpgZatGi0E4Esx1KdW+qXgYcmRz328Y+zgCQ9B1Sy1oFc7ujSS2d/TkiNhdtJOmAiFgUEbeSWu46uJQxpgAXp7XdDZgNHCvpwKSskaR2GdouBQ5Ma9sq+VNAH2BxhjavArsV7IMCjgf+neEafgA0B7oDdyWzYwXaAUuKtjGr7d57772M5a+99hpAiXegbdiwgd69e9O4cWOefPJJ6tUr/Ptf48aN6datGwsWLKBz587k5uYWe5Xm4IMP5quvvmL58uXbcFWw9957k5uby5QpUwrdGj9r1iw+/fRTjjrqqG3qtzwKlg2LbpT++9//ToMGDejQoUOh8nfeeYd99tln60yR2Y6oujfh3g78WdJZwNSyKpdgtaR/AU1ILY8VmEBqqSzTchnAUEnHkZqx+jfwLKnlrUxGAPcmt+RvBoZHxFOSBgFPpO0DuhYoetvLJOCnpBI0gHFJoiNgPnABpG7zBy6IiMERsVnSMOAfSeI0B3govdNkD9NI4IRkRuwe4E5goKQ9gPVps1FmO4wOHTpw3HHH0bdvX/bbbz8+/fRTJk+ezKhRo/i///u/rXeMFTV06FDmzp1LXl7e1uSpQPv27WnSpAm/+93v6N69O7169eK8886jVatWfPzxx8ydO5fNmzczcuTIEuPq3r07AK+88gp77713oXMzZsxg5cqV/O9//wMgPz+fXXZJ3YPRr9/XvwuOHDmSXr160a9fPwYPHszKlSu55pprOPjggwvdzn/jjTcyfPhw3nnnHXJyckqMad26dUyePBn4OmH8y1/+AkBOTs7WJK9Dhw4MGjSI66+/ni1bttC5c2deeOEFRo8ezXXXXbc11gIvv/zy1us121HVqgczSpoODIuI/AzncoHfR8S23QNbiSS9BJwWEWuqabxfAJ9GRKYltq38YEbbXlTkwYyjRo1i8uTJLFiwgA8//JC6devSrl07+vfvz9ChQ7duas7Ly+Occ87ZmjT06NGDGTNmZOxz2rRp9OjRA4ClS5cyfPhwpk6dytq1a2nZsiWdO3fmggsu4JRTTik1tq5du9K+fXvGjCn8e1hpYxf9O/fZZ5/l+uuvZ9GiRTRu3JhTTz2V2267jT32+HpV/LLLLuPuu+/mf//7H82aNSva5VbLli1jv/32y3hu4MCBhZ4l9OWXX3LTTTfx6KOP8uGHH5KTk8NFF13EkCFDCrV777332HfffZkwYUKZT9422x6pnA9m3CESIklXAj8DzoyIl2oitiLxdCU1Y7OwmsY7B3isjI3iTohsu1GRhGh7lpeXx5AhQ/jggw9o1KhRlY1zzDHH0KlTp2L7oKrDrbfeyv3338/bb79N3bp1y25gtp0pb0JUqx7MGBE9Ms0ORcTIiNh3e0iGACLi5epKhpLxxpSVDJlZ5TvrrLNo3bp1lSYq69atY8GCBVxxRaYbVKvWhg0buPPOO7npppucDNkOr1YlRGZm25O6devyyCOPVOnsUKNGjfjiiy/Yd99MT/moWsuWLWPIkCGcddZZ1T62WXXzk43NzL6Bgidh74gOPvhgDj64tBtyzXYcToiySMfWTcnfQfZumJmZVSYvmZmZmVnWc0JkZmZmWc8JkZmZmWU9J0RmZmaW9bypOossen8tOVdOqukwzMxq3I7ycFCrPJ4hMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrOeEKANJPSQ9kxwPknRPCfX6SLq+SFk/SSEpN63sKklvSXpdUq8S+tpP0suS3pT0J0k7VyDeiyWdU976ZmZmVljWJkSSKuORA5cD96X1uSvwc+DltLL2wI+BQ4GTgfsk1c3Q163A7yOiLbAaOK8CcTySjGtmZmbboEoSIkk5kpZKekjSEklTJDVMzk0vmD2R1ELSsuR4kKSnJU2U9E4y6/FLSfMkzZa0e1r7OyT9S9JiSV0k1UlmVlomdeokMzItisR1o6QHJU0B/iCpgaQxkhYl4xxXgWtsB2yMiI/Tim8GfgNsSCvrDfwxIjZGxDvAW0CXIn0JOB74S1L0KNAnw5h3FcxISeolaaakOhGxDlgmqUvRNmZmZla2qpwhagvcGxGHAmuAH5ajTQdgAKmE4RZgXUQcAcwCzk6r1zgijgEuBB6JiC3AWODM5HxPYEGRZKXAkUDviBgAXAQQER2B/sCjkhqU8/qOBeYWvJF0BLBPRDxTpF5r4L2098uTsnTNgTURsamUOgBXAmckidtdwDnJtQPkA92KNpD0E0n5kvI3r1tbviszMzPLMlWZEL0TEfOT4zlATjnaTIuIzyJiJbAWmJiULyrS/gmAiJgJNJHUjNSyUUHSdC4wpoQxJkTE+uT4O8BjSV+vAe8C7coRJ0ArYCWkZqSA3wOXZqinDGWxDXVIZoLOB54H7omIt9NOfwTslaHNgxGRGxG5dRs1zXQdZmZmWa8qE6KNaceb+fprQjaljVt0Nia9zZa091so/DUjRZOFiIj3gA8lHQ90BZ4tIa4v0o4zJSLltZ6v49+V1OzW9GQJ8GhgQrI0uBzYJ63d3sCKIn19DDRL29eUqU6BjsAnFE9+GiQxmZmZWQXVxKbqZaSWrQD6bWMfZwBI+g6wNiIK1oJGk1o6+3NEbC5HPzNJltmSPUFtgNfLGcNS4ECAiFgbES0iIicicoDZwOkRkQ9MAH4sqb6k/UgtJb6S3lFEBDCNr38eA4G/FR1Q0r6kZqGOAL4nqWva6XbA4nLGbmZmZmlqIiG6HfiZpH8BLcqqXILVSftRFL4bawKwCyUvlxV1H1BX0iLgT8CgiNhYRpsCM4Ejkg3RJYqIJcCfgX8DzwEXFSRrkiZLKpjpuQL4paS3SO0peji9n2Sch4FhEbGC1HWPTtvzdCzwQjljNzMzszRKTU7UHpKmk0oK8jOcyyV163qxzcVVFMudwMSIqNFEJNnQ/cuIOKu0evVbtY1WA++opqjMzLZf/rb77CFpTkTkllVvh3kOkaQrgb8CV1XjsL8CGlXjeCVpAVxX00GYmZnVVpXxcMJqFRE9SigfCYys5lg+JLVMV6Mi4vmajsHMzKw222FmiMzMzMy2lRMiMzMzy3q1bsnMtl3H1k3J90ZCMzOzYjxDZGZmZlnPCZGZmZllPSdEZmZmlvWcEJmZmVnW86bqLLLo/bXkXDmppsMwMzPLqCafIO4ZIjMzM8t6TojMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3rbVUIkabOk+ZIWS3pSUqNS6uZIGpD2fpCke6oorj6Srk+Ou0uaK2mTpH5pdfaVNCeJf4mkC9LO9Ze0SNJCSc9JapFhjN7J+fmS8iV9pwLxtZT03De9TjMzs2xV5QmRpIrc2r8+IjpFRAfgS+CCUurmAANKOV+ZLgfuS47/CwwCHi9S5wPgmIjoBHQFrpS0V3L9dwLHRcRhwELg4gxj/AM4PGl/LjC6vMFFxErgA0nHlv+SzMzMrECZCVEyE/OapEeTGYy/FMzcSFpWMNshKVfS9OT4RkkPSpoC/EFSXUm3SXo16eOn5YjtReBASTdLGpIWzy2Sfg6MBLolMyq/SE7vlczAvCnpN2ltCmZoFku6Na3886S/BZJmS9ojw/W3AzZGxMcAEbEsIhYCW9LrRcSXEbExeVufr3+2Sl6NJQloAqwoOk5EfB4RkbxtDETROpKOSn5+DSQ1TmaiOiSnnwbOzPiTNDMzs1KVd4boIODBZIbjU+DCcrQ5EugdEQOA84C1EXEUcBRwvqT9SmqYzKp8D1gEPAwMTMrrAD8GxgFXAi8mM0q/T5p2As4AOgJnSNpH0l7ArcDxyfmjJPVJ6jcGZkfE4cBM4PwM4RwLzC3H9ZKMtxB4D7g1IlZExFfAz5JrWQG0T64pU/u+kl4DJpGaJSokIl4FJgAjgN8AYyNicXI6H+hWnjjNzMyssPImRO9FxD+T47FAefa3TIiI9cnxScDZkuYDLwPNgbYZ2jRM6uSTWpp6OCKWAZ9IOiLpZ15EfFLCmP+IiLURsQH4N7AvqQRsekSsjIhNpJKp7kn9L4FnkuM5pJbhimoFrCzH9RIR7yVJ44HAQEl7SNqJVEJ0BLAXqSWzq0poPz4iDgb6ADeXMMxNwIlALqmkqMBHSf+FSPpJsicpf/O6teW5DDMzs6xT3v09RZdvCt5v4uukqkGROl+kHQu4JCL+XsY465M9NEWNJrVvZ0/gkVLab0w73kzq+lRK/a/SlqkK6heLCWhaSh/FRMQKSUtIzdi8m5S9DSDpz6Rmt0prP1PSAZJaFCzVpdkd2AXYidTPvODn3CCJtWhfDwIPAtRv1bbYMpyZmZmVf4aojaRvJ8f9gZeS42WklsYAflhK+78DP0tmS5DUTlLjCsQ5HjiZ1GxPQVL1GbBrOdq+DHxXUgtJdZP4Z1Rg7KWkZnxKJWlvSQ2T491ILbW9DrwPtJfUMql6YtJn0fYHJnuMkNQZ2BnINBP2IHAdqZmuW9PK2wGLM9Q3MzOzMpR3hmgpqSWgB4A3gfuT8uHAw5KuJpV4lGQ0qeWouck/+itJLQuVS0R8KWkasCYiNifFC4FNkhYAecDqEtp+IOkqYBqp2aLJEfG38o5Nam/RbyUpIkLSUaQStN2A70saHhGHAock9SIZ5/aIWAQgaTgwU9JXpGaMBiXlFyQxjiKVUJ6d1FkPnJE2e0VS/2xgU0Q8niR3/5J0fERMBY4jtffIzMzMKkhF/s0tXkHKAZ5JboWvEclm6rnAjyLizRoY/05gYkS8UN1jl5ekmaQ2sWdMDCG1ZNZq4B3VGJWZmVn5VcW33UuaExG5ZdXbrh7MmImk9sBbpDZMV3sylPgVUOJDImtashz3u9KSITMzMytZmUtmyV1eNTY7FBH/BvavqfGTGD4kdbv7dil5MOPTNR2HmZlZbbXdzxCZmZmZVTUnRGZmZpb1KvI9Y1bLdWzdlPwq2LBmZmZW23mGyMzMzLKeEyIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsl6ZX91hOw5JK0l9l9q2aAqsrcRwtscYKqv/b9LPtrTQxMIBAAAgAElEQVStSJvy1i1PvRbAx+Uct7by5756+tkePvf+zKfsiJ/5fSOiZZm1IsIvv8p8AQ/u6DFUVv/fpJ9taVuRNuWtW556QH5Nfyaq+uXPffX0sz187v2Zr9zPQ22MwUtmVl4TazoAqj6Gyur/m/SzLW0r0qa8dbeH/97bg+3h5+DP/TdvU56628N/6+3B9vBzqJEYvGRmZttEUn6U4xukzXYU/szv2DxDZGbb6sGaDsCsmvkzvwPzDJGZmZllPc8QmZmZWdZzQmRmZmZZzwmRmZmZZT0nRGZWqST1kPSipFGSetR0PGbVRVJjSXMknVbTsVjFOSEys60kPSLpI0mLi5SfLOl1SW9JurKMbgL4HGgALK+qWM0qSyV97gGuAP5cNVFaVfNdZma2laTupJKZP0REh6SsLvAGcCKpBOdVoD9QF/h1kS7OBT6OiC2S9gB+FxFnVlf8Ztuikj73h5H6ao8GpP4feKZ6orfKUq+mAzCz7UdEzJSUU6S4C/BWRPwHQNIfgd4R8WugtKWB1UD9qojTrDJVxude0nFAY6A9sF7S5IjYUqWBW6VyQmRmZWkNvJf2fjnQtaTKkn4A9AKaAfdUbWhmVaZCn/uIuAZA0iCSWdIqjc4qnRMiMyuLMpSVuNYeEU8BT1VdOGbVokKf+60VIvIqPxSrDt5UbWZlWQ7sk/Z+b2BFDcViVl38uc8yTojMrCyvAm0l7SdpZ+DHwIQajsmsqvlzn2WcEJnZVpKeAGYBB0laLum8iNgEXAz8HVgK/DkiltRknGaVyZ97A992b2ZmZuYZIjMzMzMnRGZmZpb1nBCZmZlZ1nNCZGZmZlnPCZGZmZllPSdEZmZmlvWcEJlZlZK0WdJ8SYslTZTU7Bv01UPSMWnvL5B0duVEWq7xu0laklxPw3LUHy2p/TaOlSNpcQXb/Cut7YBtGbeUvq/ONJbZjsLPIcoiLVq0iJycnJoOw8zMrNrMmTPn44hoWVY9f7lrFsnJySE/P7+mwzAzM6s2kt4tTz0vmZmZmVnWc0JkZmZmWc8JkZmZmWU97yHKIoveX0vOlZMqtc9lI0+t1P7MzMxqgmeIzMzMLOs5ITIzM7Os54TIzMzMsp4TIjMzM8t6TojMzMws6zkhMjMzs6xX6xMiSYMk7VUF/W79UkZJyyS1qEDbOyR1T44vlvSWpCitD0kDJb2ZvAZWMNYXJO1WkTZmZmb2tVqfEAGDgEpPiCJicET8u6LtJO0OHB0RM5OifwI9gRK/SyVpcwPQFegC3FDBBOcx4MKKxmpmZmYpNZYQScqR9JqkRyUtlPQXSY0knSBpfFq9EyU9JamupDxJiyUtkvQLSf2AXGCcpPmSGko6UtIMSXMk/V1Sq6Sf6ZJulfSKpDckdUvK60q6PelzoaRL0urnFom5saRJkhYkcZyR4dL6Ac8VvImIeRGxrIwfRy/g+YhYFRGrgeeBk4uM3VTS65IOSt4/Ien85PQEoH8ZY5iZmVkJanqG6CDgwYg4DPiU1CzHVOAQSS2TOucAY4BOQOuI6BARHYExEfEXIB84MyI6AZuAu4F+EXEk8AhwS9p49SKiCzCU1IwMwE+A/YAjkjjGlRLvycCKiDg8IjqQlvikORaYU6GfArQG3kt7vzwp2yoi1gIXA3mSfgzsFhEPJedWA/UlNS/asaSfSMqXlL953doKhmVmZpYdajohei8i/pkcjwW+ExFBagno/0lqBnwbeBb4D7C/pLslnUwqgSrqIKAD8Lyk+cC1wN5p559K/pwD5CTHPYFREbEJICJWlRLvIqBnMtPULUlSimoFrCztojNQhrIoVhDxfBLDvcDgIqc/IsPSYUQ8GBG5EZFbt1HTCoZlZmaWHWo6ISr6j37B+zHA/yO1DPRkRGxKZkEOB6YDFwGjM/QnYElEdEpeHSPipLTzG5M/N/P197gpQxyZg414AziSVFLya0nXZ6i2HmhQnv7SLAf2SXu/N7CiaCVJdYBDkjF2L3K6QVJuZmZmFVTTCVEbSd9OjvsDLwFExApSCcG1QB5AcodWnYj4K3Ad0Dlp9xmwa3L8OtCyoE9JO0k6tIwYpgAXSKqXtCmaaGyV3M22LiLGArenxZBuKXBgGWMW9XfgJEm7JZupT0rKivpF0n9/4BFJOyVxCdgTWFbBcc3MzIyaT4iWAgMlLSQ143F/2rlxpJbUCu70ag1MT5bC8oCrkvI8YFRSXpfUpuZbJS0A5gPHlBHDaOC/wMKkzYBS6nYEXknGugYYkaHOJKBHwRtJP5e0nNSsz0JJo5Py3ILjZJnuZuDV5HVT0aU7Se1ILZNdGhEvAjNJJYyQmrWaXbDsZ2ZmZhWj1JadGhhYygGeSTYnZzp/DzAvIh6uzrgqg6SXgNMiYk01jXcnMCEi/lFavfqt2kargXdU6tjLRp5aqf2ZmZlVJklzIiK3rHo1PUOUkaQ5wGGkNlrXRpcCbapxvMVlJUNmZmZWsnplV6kaybN5Ms4OJbfM11oR8XI1j/dQdY5nZma2o9kuZ4jMzMzMqpMTIjMzM8t6NbZkZtWvY+um5HsTtJmZWTGeITIzM7Os54TIzMzMsp4TIjMzM8t6TojMzMws63lTdRZZ9P5acq6cVNNhbPf89G0zs+zjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsp4TIjMzM8t6tTYhkrSnpD9KelvSvyVNltSulPqdJJ2S9v5GScOqKLahks5Ojm+T9JqkhZLGS2qWoX4DSa9IWiBpiaThFRzvj5LaVlb8ZmZm2aZWJkSSBIwHpkfEARHRHrga2KOUZp2AU0o5X1mx1QPOBR5Pip4HOkTEYcAbwFUZmm0Ejo+Iw5M4T5Z0dAWGvR+4fNujNjMzy261MiECjgO+iohRBQURMT8iXpT0mKTeBeWSxkk6HbgJOEPSfElnJKfbS5ou6T+Sfp7W5peSFievoUlZjqSlkh5KZnGmSGqYIbbjgbkRsSmJa0rBMTAb2Ltog0j5PHm7U/KK9DqS6kl6VVKP5P2vJd2SnH4R6JkkY2ZmZlZBtTUh6gDMKeHcaOAcAElNgWOAycD1wJ8iolNE/CmpezDQC+gC3CBpJ0lHJu27AkcD50s6IqnfFrg3Ig4F1gA/zDD+saXEdi7wbKYTkupKmg98BDwfES+nn0+SqkHA/ZJOBE4GhifntgBvAYdn6PcnkvIl5W9et7aEsMzMzLJbbU2IShQRM4ADJX0L6A/8NW2GpqhJEbExIj4mlYjsAXwHGB8RXySzNk8B3ZL670TE/OR4DpCToc9WwMqihZKuATYB40qIe3NEdCI1g9RFUocMdZYAjwETgXMj4su00x8Be2Vo82BE5EZEbt1GTTMNbWZmlvVqa0K0BDiylPOPAWeSmukZU0q9jWnHm0l9lYkqWL+o9UCD9AJJA4HTgDMjIjK02Soi1gDTSc0AZdKR1OxU0f1SDZKxzczMrIJqa0I0Fagv6fyCAklHSfpu8jYPGApbZ1UAPgN2LUffM4E+khpJagz0JbVHp7yWAgemxXUycAVwekSsy9RAUsuCu8+SfUk9gdcy1PsB0BzoDtxV5I61dqQSRTMzM6ugWpkQJbMsfYETk9vulwA3AiuS8x+SSkzSZ4emkdpEnb6pOlPfc0klVK8ALwOjI2JeBcJ7llTCUuAeUonY88nYowAk7SVpclKnFTBN0kLgVVJ7iJ5J71RSC2AkcF5EvJH0e2dybg9gfUR8UIE4zczMLKEyVnBqJUmNgEVA54io9p3EksYDl0fEm9U03i+ATyPi4dLq1W/VNloNvKM6QqrV/G33ZmY7DklzIiK3rHq1coaoNJIKlpvurolkKHElqVmf6rIGeLQaxzMzM9uh7HDPrYmIF4A2NRzD68Dr1TheaRvHzczMrAw73AyRmZmZWUU5ITIzM7Ost8MtmVnJOrZuSr43DJuZmRXjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsp43VWeRRe+vJefKSTUdRq3hJ1abmWUPzxCZmZlZ1nNCZGZmZlnPCZGZ1Zinn36a7t27861vfYuGDRuy77770qdPH5577rmtdfLy8pDEsmXLqjW2Sy65hO9///uFyq6++mpOOukkmjdvjiTy8vKKtfvggw+46qqryM3NpWnTprRs2ZITTjiBmTNnFqu7efNmfv/739OhQwcaN25Mq1at6Nu3LwsXLtymmD/77DOGDRtGjx49aNKkCZKYPn16sXpvvPEGQ4YM4bDDDmOXXXahVatWnH766SxYsKBY3XXr1nHDDTfQrl07GjZsyD777MPZZ59d7L9H7969ueiii7YpbrPtgRMiM6sRd911F3379qVt27Y8/PDDTJo0iWuvvRaAqVOnbq136qmnMmvWLFq1qr7vS3777bd54IEHuOGGGwqV33333axfv57TTjutxLZz5szhT3/6E7179+Yvf/kLeXl5NGjQgB49evDMM88UqnvdddcxbNgw+vTpw8SJE7nzzjt5++23Oe6441i+fHmF4/7kk0945JFHqFevHieeeGKJ9aZMmcK0adMYOHAgEydO5L777mPlypV07dqVOXPmFKo7ePBgbrvtNs4//3wmT57MiBEjmDlzJieccAKff/751no33ngjDz30EG+88UaF4zbbHigiajoGqyb1W7WNVgPvqOkwag1vqq5abdq04cgjj2T8+PHFzm3ZsoU6dWru97VLLrmE2bNn8+qrrxYqL4jrrbfeom3btowZM4ZBgwYVqrNmzRp22WUX6tX7+p6VTZs2ceihh7LHHnsUminaa6+96NGjB48//vjWstdee41DDjmEUaNG8dOf/rRCcUcEkgB44YUXOPHEE5k2bRo9evQoVO/jjz/eOstVYO3ateTk5PD973+fP/zhDwCsX7+eXXfdlcsvv5xf/epXW+s+99xzfO973+O5556jV69eW8u7dOlCbm4u9913X4XiNqtKkuZERG5Z9bJ6hkhSjqTFyXGupLsqqd87JHVPjsdJel3SYkmPSNqplHZNJL0v6Z4KjveCpN2+adxm1WnVqlXsueeeGc+lJ0NFl8wGDRqEpIyv9OWhBQsWcPrpp7PbbrvRsGFDjj32WF588cUy49q4cSNjx45lwIABpcZVkmbNmhVKhgDq1atHp06deP/99wuVf/nllzRp0qRYe0glXxWVnuCUpkWLFsXqNm3alHbt2hWKcdOmTWzevLncMf74xz9m3LhxrF+/vsKxm9W0rE6I0kVEfkT8/Jv2I2l34OiIKPg1cBxwMNARaAgMLqX5zcCMbRj2MeDCbWhnVmO6dOnCo48+ym233VahZZbrrruOWbNmFXode+yxNGrUiDZt2gAwd+5cjjnmGFatWsVDDz3EX//6V5o3b07Pnj2LLQkVNXv2bNasWUO3bt2+0fWl+/LLL5k1axaHHHJIofILL7yQsWPH8re//Y1PP/2U//znP1x44YXsvffenHHGGZU2fnmsWrWKxYsXF4px11135ayzzuKuu+5i2rRpfP755yxZsoTLLruMww8/nBNOOKFQH927d+fTTz9l1qxZ1Rq7WWWo9IQomXVZKukhSUskTZHUMDk3XVJuctxC0rLkeJCkpyVNlPSOpIsl/VLSPEmzkySjoP0dkv6VzLh0kVRH0puSWiZ16kh6S1KLInF9V9L85DVP0q5FzveQ9ExyvIukMZIWSVoo6YdJ+UmSZkmaK+lJSbtk+BH0A7buCI2IyZEAXgH2LuHndiSwBzClhPNNk5mmg5L3T0g6Pzk9Aehfwn8Ss+3SqFGjOPDAA7n88ss56KCDaNGiBf3792fKlIz/C2x1wAEHcPTRR299vfTSS8yaNYtx48ax//77A3DZZZfRpk0bpk6dSr9+/TjllFMYP348+++/PzfffHOp/c+ePRtJHHbYYZV2rTfeeCPLly/niiuuKFR+0003cdVVV/GDH/yApk2bcsABB7BkyRKmT5/O7rvvXmnjl8cll1xCRDB06NBC5WPGjKFv374cf/zx7LrrrnTo0IGvvvqK559/np133rlQ3cMPP5w6deowe/bs6gzdrFJU1QxRW+DeiDgUWAP8sBxtOgADgC7ALcC6iDgCmAWcnVavcUQcQ2pG5JGI2AKMBc5MzvcEFkTEx0X6HwZcFBGdgG5AaXO61wFrI6JjRBwGTE0SrGuBnhHRGcgHfpmh7bFAsV9Bk6Wys0hLltLO1QF+C1xWUkARsRa4GMiT9GNgt4h4KDm3GqgvqXkp12S2XWnXrh3z5s1jxowZXHPNNXTq1Inx48fTq1cvRowYUa4+Jk6cyBVXXMGtt95Knz59gNS+lxkzZvCjH/2IOnXqsGnTJjZt2kRE0LNnz4x3e6VbsWIFTZo0KfaP/bZ6/PHHGTlyJNddd12xWaf777+fESNGcO211zJt2jSefPJJdt11V0466SRWrFhRKeOXx69//Wsef/xx7rnnHg488MBC56699lrGjh3L7bffzowZM3jsscf45JNP+N73vscXX3xRqO5OO+1E06ZNqzV2s8pSVU+qfici5ifHc4CccrSZFhGfAZ9JWgtMTMoXAem/qj0BEBEzkz03zYBHgL8BdwDnAmMy9P9P4HeSxgFPRcTyUtbbewI/LngTEaslnQa0B/6ZtNuZVLJWVCtgZYby+4CZEZFpE8OFwOSIeK+0PQAR8bykHwH3AocXOf0RsBfwSXqhpJ8APwGo26RliX2b1YS6devSvXt3unfvDqSSkZNPPpnhw4dz0UUXsdtuJW+NW7BgAQMGDOC8885j2LBhW8tXrVrF5s2bufnmm0ucDSpt0/aGDRuoX7/+N7iqr02cOJFBgwZx3nnnMXz48ELnVq1axS9+8Qsuu+yyQueOP/54cnJyuO222/j9739fKXGUZtSoUVx99dWMGDGCc889t9C5JUuWMHLkSEaPHs155523tbxr1660a9eO0aNHM2TIkEJtGjZs6D1EVitVVUK0Me14M6m9MwCb+HpWqkEpbbakvd9C4TiL3hYXSSLxoaTjga58PVuUXmmkpEnAKcBsST2BDSXErwzjCHg+IspamlpPkWuTdAPQEijplpFvA90kXQjsAuws6fOIuLJIP3WAQ5IxdgfS78ttQIZZr4h4EHgQUneZlRG7WY3aa6+9GDx4MEOGDOHNN9+kS5cuGet9+OGHnH766Rx99NHF7mhq1qwZderU4aKLLuLss8/O2L60zdHNmzdn9erV234RiX/84x/86Ec/om/fvjzwwAPFzr/xxhts3LiRo446qlD57rvvzgEHHMDSpUu/cQxleeyxx7jwwgu59NJLueaaa4qdX7RoEUCxGNu2bUuzZs0yxrhq1SpatGhRrNxse1fdm6qXAUcmx/22sY8zACR9h9Sy1tqkfDSppbM/R8Tmoo0kHRARiyLiVlLLXQeXMsYUUstTBW13A2YDx0o6MClrJKldhrZLgQPT2g4GegH9k+W9YiLizIhoExE5pJb2/lA0GUr8Ium/P7D1jjWlppX2JPXzNasV3nvvvYzlr732GkCJd6Bt2LCB3r1707hxY5588slid3Q1btyYbt26sWDBAjp37kxubm6xV2kOPvhgvvrqq216DlCBWbNm0bt3b0444QTGjh2bMQEruL5XXnmlUPmqVat46623aN269TaPXx7jx4/nnHPOYfDgwdx+++0Z65QU4xtvvMGaNWuKxfi///2PDRs2cNBBB1VN0GZVqLq/3PV24M+SzgKmllW5BKsl/QtoQmp5rMAEUktlmZbLAIZKOo7UjNW/gWdJLW9lMgK4N7klfzMwPCKekjQIeEJSwXz6tUDR22MmkZoJGp28HwW8C8xKlsOeioibks3lF0REaXedbZUkX4OBLhHxmaSZyfg3kEoyZ0fEpvL0ZbY96NChA8cddxx9+/Zlv/3249NPP2Xy5MmMGjWK//u//9t6x1hRQ4cOZe7cueTl5W1Nngq0b9+eJk2a8Lvf/Y7u3bvTq1cvzjvvPFq1asXHH3/M3Llz2bx5MyNHjiwxroLlu1deeYW99y58D8SMGTNYuXIl//vf/wDIz89nl11S91b065f6He+1117j1FNPpUWLFlx22WXF7mo7+uijAcjJyeG0007jtttuo06dOnz3u9/lk08+4Te/+Q0bN27kZz/72dY206dP57jjjsv43KOinn32Wb744outszszZszg448/pnHjxnzve98DYObMmfTv35/DDjuMQYMGFdoEXb9+fY444ggAunXrxuGHH86ll17K6tWryc3N5b///S8jRoygadOmDBw4sNDYL7/8cqGfoVltUqsezChpOjAsIvIznMsFfh8RlXev7DaS9BJwWkSsqabx7gQmRMQ/SqvnBzNWjB/MWLVGjRrF5MmTWbBgAR9++CF169alXbt29O/fn6FDh27d1JyXl8c555zDO++8Q05ODj169GDGjMxPp0h/COHSpUsZPnw4U6dOZe3atbRs2ZLOnTtzwQUXcMopp5QaW9euXWnfvj1jxhT+/aq0sQv+Li2ItyTpf+euW7eO3/72tzzxxBO8++67NGnShM6dO3PDDTcUWi6cNGkSp512Gs8++ywnn3xyqbHn5OTw7rvvFivfd999tz7L6cYbbyy2pylTPUg9/fpXv/oVEyZMYPny5bRo0YJjjjmGm266qdhM0Pnnn8+8efPIzy/2V7RZjSnvgxl3iIRI0pXAz4AzI+KlmoitSDxdgfURsW1fSFTx8c4vuOOsNE6IKsYJUfbKy8tjyJAhfPDBBzRq1Kimw+Hqq69mwoQJLFq0qNwPX6xuGzZsoFWrVtx+++2FNmCb1bTyJkS16sGMEdEj0+xQRIyMiH23h2QIICJerq5kKBmvzGTIzMrvrLPOonXr1tvNV1DMmDGDq6++ertNhgAeeOABvvWtbxVbRjOrLap7D5GZ2Xavbt26PPLII8ydO7emQwHgn//8Z02HUKb69euTl5dXbJO7WW1Rq5bM7JvxklnFeMnMzKz2K++SmVP5LNKxdVPy/Y+8mZlZMbVqD5GZmZlZVXBCZGZmZlnPCZGZmZllPSdEZmZmlvW8qTqLLHp/LTlXTqrpMGw74bvozMy+5hkiMzMzy3pOiMzMzCzrOSEyMzOzrOeEyMzMzLKeEyIzMzPLek6IKkBSD0nPJMeDJN1TQr0+kq5Pq7dS0vzkNbiENv0lLZK0UNJzklpUIK7bJR2/LddkZmZmToiKkVQZjyK4HLgv7f2fIqJT8hpdwph3AsdFxGHAQuDiCox3N3DlNwnYzMwsm1VrQiQpR9JSSQ9JWiJpiqSGybnpknKT4xaSliXHgyQ9LWmipHckXSzpl5LmSZotafe09ndI+pekxZK6SKoj6U1JLZM6dSS9VXT2RdKNkh6UNAX4g6QGksYkMzbzJB1XgWtsB2yMiI8r8qNJXo0lCWgCrMjQ998knZ0c/1TSOICIeBdoLmnPCoxpZmZmiZqYIWoL3BsRhwJrgB+Wo00HYADQBbgFWBcRR/z/9u49Tus5///446nohKJYyWFYhZRCck5WTrt+YvElFq1YZ6tdS+yyrPWVZb/rsEiR7BbLsmw5HxO2MKWD5BCycozIKVJevz+u99Q111xzzUzNzNV0Pe+329z6XJ/36fWZ+dxmXr3f7891AROBY7PqtYmIXYFTgZER8T0wGjg6lfcDplWTrOwA9I+Io4DTACKiOzAAuFVSy1pe327AlJxzh6alsLskbZzbICK+A04BZpBJhLoCN+fp+xfAhZL2AH4NnJFVNiWNbWZmZnVUjITorYiYmo4nA2W1aPNkRHwREfOABcC4dH5GTvvbASJiArC2pHbASJYlTccDt1QzxtiIWJiOdwf+nvp6BXgb6FKLOAE6AvOyXo8DytJS2GPArbkNJK1OJiHaDtiQzJLZebn1IuJD4ELgSeDXETE/q/ij1Da3719IKpdUvuTrBbW8BDMzs9JSjITo26zjJSz7+JDFLIsndzYmu833Wa+/p/LHj0ROu4iId4AP06bjnYAHq4nrq6xjVRt9zRaSFX9EfBIRFfGOIDMTlatnqvtGRARwJ7BrNf13Bz6havLTMo1dSUQMj4heEdGrWeu2dboQMzOzUrEybaqew7Jk4bDl7OMIAEm7AwsiomJK5CYyS2d3RsSSWvQzgbTMlvYEbQK8WssYZgFbVLyQ1DGr7KBUnutdoGvFXidgn3z1JPUGDiAzk3S2pM2yirsAL9UyRjMzM8uyMiVEVwKnSPoPUOtHznN8mtoPAwZlnR8LrEn1y2W5rgeaSZoB3AEMzJrlqckEYLu0ORrgzLSBfBpwJjCwoqKkqQAR8R5wMTBB0nQyM0b/m92ppBZkZpiOT/V/DYxUxupkkrDyWsZoZmZmWZRZoWn6JI0Hzo6IKklBenrtLxGxRyPFcjUwLiIea6TxDgG2j4gLCtVr0bFzdDzuqsYIyZoAf9q9mZUCSZMjoldN9VamGaIGIWkIcDd5Nik3oP8FWjfieM2BPzfieGZmZquU+ngTwpVCRPSt5vxQYGgjx/IhmWW6xhrvn401lpmZ2apolZ8hMjMzM6uJEyIzMzMreavMkpnVrHuntpR7I62ZmVkVniEyMzOzkueEyMzMzEqeEyIzMzMreU6IzMzMrOR5U3UJmfHuAsqG3F/sMMzMzPIq5jvoe4bIzMzMSp4TIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3krTUIkaYmkqZJekvRPSa0L1C2TdFTW64GS/tpAcR0s6cKscealOKdKOiGr3uUp9pckHZGnn2slfVngehZm9TusjjE+Jmmdul6bmZmZZTRoQiSpLu9ztDAiekZEN2ARcHKBumXAUQXK69M5wPVZr+9IcfaMiJsAJP0E2B7oCewE/EbS2hUNJPUC2tUwzhtZ/Ra69nz+DpxaxzZmZmaWFEyI0szFK5JulTRd0l0VMzeS5kjqkI57SRqfji+SNFzSI8DfJDWTdOd0l3UAACAASURBVIWkF1IfJ9UirqeBLSRdIumXWfFcKulMYCiwR5pNGZyKN5T0kKTXJf0pq80ASTPSzM3lWee/TP1NkzRJ0g/yXH8X4NuI+LiGeLsCT0XE4oj4CpgG7J/6aAZcQSaxWm6S2kp6VdKW6fXtkk5MxWOBASvSv5mZWSmrzQzRlsDwiNgW+JzazUTsAPSPiKOAQcCCiNgR2BE4UdJm1TVMs0oHADOAm4Hj0vnVgCOBMcAQ4Ok0m/KX1LQncATQHThC0saSNgQuB36UyneUdHCq3waYFBE9gAlARXKRbTdgSs65Q7OSw43TuWnAAZJapyRxL6Ci7HRgbES8X+gbBmwm6UVJT0naI7cwIhakvkZJOhJYJyJGpLJPgRaS2ue2k/QLSeWSypd8vaCGEMzMzEpTbRKidyLi2XQ8Gti9Fm3GRsTCdLwvcKykqcBzQHugc542rVKdcuC/wM0RMQf4RNJ2qZ8XI+KTasZ8PCIWRMQ3wMvApmQSsPERMS8iFpNJpvqk+ouA+9LxZDLLcLk6AvOyXo8DylJy+BhwK0BEPAI8APwHuB2YCCxOCdnhwLXVxFzhfWCTiNgO+BVwW/aSW4WIeJRMongdcEJO8UfAhnnaDI+IXhHRq1nrtjWEYWZmVppqs8cnqnm9mGUJVcucOl9lHQs4IyIermGchRHRM8/5m4CBwAbAyALtv806XkLm2lSg/ncRETn1q8QELM0icpKxEWRmnyrKLgUuBZB0G/A6sB2wBTBbEkBrSbMjYovsQSLi24r4I2KypDeALmSSw6XSLNnWKa51gblZxS3TeTMzM6uj2swQbSJpl3Q8AHgmHc8hszQGcGiB9g8Dp0haHTL7ciS1qUOM95DZj7Nj6gvgC2CtWrR9DthTUoe0l2cA8FQdxp5FJqEBQFLHrLKDUjlpn1T7dLwtsC3wSETcHxEbRERZRJQBX+cmQ6nNeik+JG1OZgbtzTzxDE5jDgBGZn1PRSZhnFOHazMzM7OkNjNEs4DjJN1IZtbjhnT+YuBmSeeTSTyqcxOZ5agp6Q/3PODgAvUriYhFkp4EPouIJen0dDJLUtOAUcCn1bR9X9J5wJNkZoseiIh/13ZsMnuL/ixJaTbpTEkHkZkdm09m5gpgdeDpNAv0OfCztERXrdRPr4i4kMwy3h8kLSYzW3VyRMzPqd+FzDJZ74j4QtIE4HfA78kkppNqGtPMzMzy07JVozyFUhlwX3oUvijSMtEU4PCIeL0I418NjIuIxxp77NpKMY6NiMcL1WvRsXN0PO6qRorKzMysbuYM/Um99ylpckT0qqneSvPGjPlI6grMJrNhutGToeR/gWrfJHIl8VJNyZCZmZlVr+CSWXrKq2izQxHxMrB5scZPMXxI5n1+VloVj9+bmZnZ8lmpZ4jMzMzMGoMTIjMzMyt5dfmsMWviundqS3kDbFgzMzNr6jxDZGZmZiXPCZGZmZmVPCdEZmZmVvKcEJmZmVnJc0JkZmZmJa/gR3fYqkXSPODtFeymLbCgEdrVtn5N9QqVFyrrAHxci/FXFsv7cynWOKvSfVSovKndR9A491J9jtEY91J93Uc11fHvpIYZZ9OIWK/GniLCX/6q9RcwvDHa1bZ+TfUKlddQVl7s73Vj/FyKNc6qdB8VKm9q91F9/owba4zGuJfq6z6qqY5/JxV3HC+ZWV2Na6R2ta1fU71C5ct7LSujxrqW+hpnVbqP6jJOU9AY11KfYzTGvVRf91FNdXwfFXEcL5mZ5SGpPGrx6chmhfg+svrie6nheYbILL/hxQ7AVgm+j6y++F5qYJ4hMjMzs5LnGSIzMzMreU6IzMzMrOQ5ITIzM7OS54TIbDlIaiNpsqQDix2LNU2StpY0TNJdkk4pdjzWdEk6WNIISf+WtG+x42mqnBBZSZE0UtJHkl7KOb+/pFclzZY0pBZdnQvc2TBR2squPu6jiJgVEScD/wP4ceoSVU/30r0RcSIwEDiiAcNdpfkpMyspkvoAXwJ/i4hu6Vwz4DVgH2Au8AIwAGgGXJbTxfHAtmTeRr8l8HFE3Nc40dvKoj7uo4j4SNJBwBDgrxFxW2PFbyuP+rqXUrs/A2MiYkojhb9KaV7sAMwaU0RMkFSWc7o3MDsi3gSQ9A+gf0RcBlRZEpO0F9AG6AoslPRARHzfoIHbSqU+7qPUz1hgrKT7ASdEJaieficJGAo86GRo+TkhMoNOwDtZr+cCO1VXOSJ+CyBpIJkZIidDBnW8jyT1BX4KtAAeaNDIrKmp070EnAH0A9pK2iIihjVkcKsqJ0RmoDznalxLjohR9R+KNWF1uo8iYjwwvqGCsSatrvfSNcA1DRdOafCmarPM/742znq9EfBekWKxpsv3kdUX30tF4ITILLNhsbOkzSStARwJjC1yTNb0+D6y+uJ7qQicEFlJkXQ7MBHYUtJcSYMiYjFwOvAwMAu4MyJmFjNOW7n5PrL64ntp5eHH7s3MzKzkeYbIzMzMSp4TIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3n+6I4S0qFDhygrKyt2GGZmZo1m8uTJH0fEejXVc0JUQsrKyigvLy92GGZmZo1G0tu1qeclMzMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnjdVl5AZ7y6gbMj99drnnKE/qdf+zMzMisEzRGZmZlbynBCZmZlZyXNCZGZmZiXPCZGZmZmVPCdEZmZmVvKafEIkaaCkDRug35skdU3HcyR1qEPbqyT1ScenS5otKQr1Iek4Sa+nr+PqGOtjktapSxszMzNbpsknRMBAoN4Toog4ISJerms7SesCO0fEhHTqWaAfUO1nqaQ2vwd2AnoDv69jgvN34NS6xmpmZmYZRUuIJJVJekXSrZKmS7pLUmtJe0u6J6vePpL+JamZpFGSXpI0Q9JgSYcBvYAxkqZKaiVpB0lPSZos6WFJHVM/4yVdLul5Sa9J2iOdbybpytTndElnZNXvlRNzG0n3S5qW4jgiz6UdBjxU8SIiXoyIOTV8O/YDHo2I+RHxKfAosH/O2G0lvSppy/T6dkknpuKxwIAaxjAzM7NqFHuGaEtgeERsC3xOZpbjCWBrSeulOj8HbgF6Ap0ioltEdAduiYi7gHLg6IjoCSwGrgUOi4gdgJHApVnjNY+I3sBZZGZkAH4BbAZsl+IYUyDe/YH3IqJHRHQjK/HJshswuU7fBegEvJP1em46t1RELABOB0ZJOhJYJyJGpLJPgRaS2tdxXDMzM6P4CdE7EfFsOh4N7B4RQWYJ6GeS2gG7AA8CbwKbS7pW0v5kEqhcWwLdgEclTQV+B2yUVf6v9O9koCwd9wOGRcRigIiYXyDeGUC/NNO0R0pScnUE5hW66DyU51xUORHxaIrhOuCEnOKPyLN0KOkXksollS/5Ol+4ZmZmVuyEKPePfsXrW4CfkVkG+mdELE6zID2A8cBpwE15+hMwMyJ6pq/uEbFvVvm36d8lLPvYEuWJI3+wEa8BO5BJSi6TdGGeaguBlrXpL8tcYOOs1xsB7+VWkrQasHUaY92c4pbpfG7MwyOiV0T0ata6bR3DMjMzKw3FTog2kbRLOh4APAMQEe+RSQh+B4wCSE9orRYRdwMXANundl8Aa6XjV4H1KvqUtLqkbWqI4RHgZEnNU5vcRGOp9DTb1xExGrgyK4Zss4Atahgz18PAvpLWSZup903ncg1O/Q8ARkpaPcUlYANgTh3HNTMzM4qfEM0CjpM0ncyMxw1ZZWPILKlVPOnVCRiflsJGAeel86OAYel8MzKbmi+XNA2YCuxaQww3Af8Fpqc2RxWo2x14Po31W+CPeercD/SteCHpTElzycz6TJd0Uzrfq+I4LdNdAryQvv6Qu3QnqQuZZbJfR8TTwAQyCSNkZq0mVSz7mZmZWd0os2WnCANLZcB9aXNyvvK/Ai9GxM2NGVd9kPQMcGBEfNZI410NjI2IxwvVa9Gxc3Q87qp6Hdufdm9mZiszSZMjoldN9Yo9Q5SXpMnAtmQ2WjdFvwY2acTxXqopGTIzM7PqNa+5SsNI782Td3YoPTLfZEXEc4083ojGHM/MzGxVs1LOEJmZmZk1JidEZmZmVvKKtmRmja97p7aUexO0mZlZFZ4hMjMzs5LnhMjMzMxKnhMiMzMzK3lOiMzMzKzkeVN1CZnx7gLKhtxf7DCKyu+sbWZm+XiGyMzMzEqeEyIzMzMreU6IzMzMrOQ5ITIzM7OS54TIzMzMSl6TTYgkbSDpH5LekPSypAckdSlQv6ekH2e9vkjS2Q0U21mSjk3HV0h6RdJ0SfdIaldNm/0lvSpptqQhdRzvH5I610fsZmZmpahJJkSSBNwDjI+IH0ZEV+B84AcFmvUEflygvL5iaw4cD9yWTj0KdIuIbYHXgPPytGkGXAccAHQFBkjqWodhbwDOWZG4zczMSlmTTIiAvYDvImJYxYmImBoRT0v6u6T+FecljZF0EPAH4AhJUyUdkYq7Shov6U1JZ2a1+ZWkl9LXWelcmaRZkkZIminpEUmt8sT2I2BKRCxOcT1ScQxMAjbK06Y3MDsi3oyIRcA/gP7ZFSQ1l/SCpL7p9WWSLk3FTwP9UjJmZmZmddRUE6JuwORqym4Cfg4gqS2wK/AAcCFwR0T0jIg7Ut2tgP3IJCS/l7S6pB1S+52AnYETJW2X6ncGrouIbYDPgEPzjL9bgdiOBx7Mc74T8E7W67np3FIpqRoI3CBpH2B/4OJU9j0wG+hRzbhmZmZWQFNNiKoVEU8BW0haHxgA3J01Q5Pr/oj4NiI+Bj4is+S2O3BPRHwVEV8C/wL2SPXfioip6XgyUJanz47AvNyTkn4LLAbG5GmjfJeS59pmAn8HxgHHp9mkCh8BG+YZ9xeSyiWVL/l6QZ5hzMzMrKkmRDOBHQqU/x04msxMzy0F6n2bdbyEzEeZ5EtOCtXPtRBomX1C0nHAgcDREVEl0SEzI7Rx1uuNgPeqiaE7mdmp3P1SLdPYlUTE8IjoFRG9mrVuW02XZmZmpa2pJkRPAC0knVhxQtKOkvZML0cBZ8HSWRWAL4C1atH3BOBgSa0ltQEOIbNHp7ZmAVtkxbU/cC5wUER8XU2bF4DOkjaTtAZwJDA2t5KknwLtgT7ANTlPrHUhkyiamZlZHTXJhCjNshwC7JMeu58JXESaVYmID8kkJtmzQ0+S2USdvak6X99TyCRUzwPPATdFxIt1CO9BMglLhb+SScQeTWMPA5C0oaQH0piLgdOBh1Pcd2YlcqT6HYChwKCIeC31e3Uq+wGwMCLer0OcZmZmlij/Ck7TJqk1MAPYPiIafeOMpHuAcyLi9UYabzDweUTcXKhei46do+NxVzVGSCstf9q9mVlpkTQ5InrVVK9JzhAVIqkf8ApwbTGSoWQImc3VjeUz4NZGHM/MzGyVssq9b01EPAZsUuQYXgVebcTxCm0cNzMzsxqscjNEZmZmZnXlhMjMzMxKnhMiMzMzK3mr3B4iq173Tm0p91NWZmZmVXiGyMzMzEqeEyIzMzMreU6IzMzMrOR5D1EJmfHuAsqG3F/sMFZJfgdsM7OmzTNEZmZmVvKcEJmZmVnJc0JkZmZmJc8JkZnVu3vvvZc+ffqw/vrr06pVKzbddFMOPvhgHnrooaV1Ro0ahSTmzJnTqLGdccYZ/L//9/8qnTv//PPZd999ad++PZIYNWpUte1HjBjBVlttRYsWLdhyyy0ZNmxYpfLPP/+cP/zhD+y66660b9+edu3aseuuu3LvvfcWjOuzzz5jgw02QBKPPfbYcl3bF198wdlnn03fvn1Ze+21kcT48ePz1v3++++57LLLKCsro2XLlvTo0YO77767Up3333+f8847j169etG2bVvWW2899t57byZMmFAwjjfffJPWrVsjidmzZ1cq69+/P6eddtpyXZ9ZQ3JCZGb16pprruGQQw6hc+fO3Hzzzdx///387ne/A+CJJ55YWu8nP/kJEydOpGPHjo0W2xtvvMGNN97I73//+0rnr732WhYuXMiBBx5YsP2IESM46aSTOPTQQ3nooYc4/PDDOfXUU7nhhhuW1vnvf//L9ddfz5577sno0aO544476NKlC4cccgjXXXddtX2fe+65SFqh6/vkk08YOXIkzZs3Z5999ilY94ILLuCiiy7i9NNP58EHH2TnnXfm8MMP54EHHlhaZ/Lkydxxxx3079+fu+66i1GjRtGyZUv69u3LfffdV23fp556Km3bts1bdtFFFzFixAhee+215btIswaiiCh2DNZIWnTsHB2Pu6rYYayS/JTZMptssgk77LAD99xzT5Wy77//ntVWK97/w8444wwmTZrECy+8UOl8RVyzZ8+mc+fO3HLLLQwcOLBSncWLF7PhhhtywAEHcOutty49f/zxxzN27Fjef/99Vl99db766isk0bp160rt9957b15//XX++9//Vonr2WefZd999+Xaa69l0KBBPProo/Tr16/O1xcRS5Oqxx57jH322Ycnn3ySvn37Vqr30UcfsfHGGzNkyBAuvvjiSjHOmzeP6dOnA5lZqzXXXJPmzZc9kLx48WK22WYbfvCDH+SdKbrtttsYPHgw5513HoMHD+b1119niy22qFSnd+/e9OrVi+uvv77O12hWV5ImR0SvmuqV9AyRpDJJL6XjXpKuqad+r5LUJx2PkfSqpJckjZS0ejVt/iRppqRZkq5RHf6rKOkxSevUR+xmK2r+/PlssMEGecuyk6HcJbOBAwciKe9X9rLPtGnTOOigg1hnnXVo1aoVu+22G08//XSNcX377beMHj2ao446qmBc1Zk4cSLz5s3jZz/7WaXzxxxzDJ988gnPPPMMAG3atKmSDAH06tWL9957r8r57777jpNOOokhQ4aw+eab1xhHIbX9tfHwww+zaNGiKtfys5/9jBkzZvDWW28B0K5du0rJEEDz5s3p2bMn7777bpV+P/30U371q19x5ZVX0q5du2rHP/LIIxkzZgwLFy6sVbxmjaGkE6JsEVEeEWeuaD+S1gV2joiK/zqNAbYCugOtgBPytNkV2A3YFugG7AjsWYdh/w6cugJhm9Wb3r17c+utt3LFFVfUaVnkggsuYOLEiZW+dtttN1q3bs0mm2wCwJQpU9h1112ZP38+I0aM4O6776Z9+/b069ePyZMnF+x/0qRJfPbZZ+yxxx7LdV0zZ84EoFu3bpXOb7PNNgC8/PLLBdtPmDCBrbbaqsr5P/3pTyxatIhzzjlnueJaHjNnzqRFixZVZm5qcy2LFi1i4sSJbL311lXKzjnnHLbaaiuOOeaYguP36dOHzz//nIkTJy5H9GYNo94TojTrMkvSiDTj8YikVqlsvKRe6biDpDnpeKCkeyWNk/SWpNMl/UrSi5ImpSSjov1Vkv6TZlx6S1pN0uuS1kt1VpM0W1KHnLj2lDQ1fb0oaa2c8r6S7kvHa0q6RdIMSdMlHZrO7ytpoqQpkv4pac0834LDgKU7RyPigUiA54GN8rQJoCWwBtACWB34MCe+tmmmacv0+nZJJ6biscCAgj8Ys0YybNgwtthiC8455xy23HJLOnTowIABA3jkkUcKtvvhD3/IzjvvvPTrmWeeYeLEiYwZM2bpzMlvfvMbNtlkE5544gkOO+wwfvzjH3PPPfew+eabc8kllxTsf9KkSUhi2223Xa7rmj9/PgDrrFN5MnbdddetVJ7P8OHDmTRpEuedd16l87Nnz+aPf/wj1113HS1atFiuuJbH/PnzadeuXZUZpdpcy0UXXcTcuXM599xzK51/5pln+Nvf/larZbAePXqw2mqrMWnSpOWI3qxhNNQMUWfguojYBvgMOLQWbboBRwG9gUuBryNiO2AicGxWvTYRsSuZGZGREfE9MBo4OpX3A6ZFxMc5/Z8NnBYRPYE9gEJztRcACyKie0RsCzyREqzfAf0iYnugHPhVnra7AVX+q5qWyo4hK1mqEBETgSeB99PXwxExK6fOAuB0YJSkI4F1ImJEKvsUaCGpfZ5xfyGpXFL5kq8XFLhks/rRpUsXXnzxRZ566il++9vf0rNnT+655x72228//vjHP9aqj3HjxnHuuedy+eWXc/DBBwOwcOFCnnrqKQ4//HBWW201Fi9ezOLFi4kI+vXrV+OTT++99x5rr702a6yxxnJdV8V+y7pufB4/fjxnnnkmxxxzDEcffXSlslNOOYX+/fvXuAG6vmXvNco9X8htt93G0KFDueCCCyrNtC1atIiTTjqJwYMH07Vr1xrHX3311Wnbtm3eJUSzYmmoj+54KyKmpuPJQFkt2jwZEV8AX0haAIxL52eQWUqqcDtAREyQtLakdsBI4N/AVcDxwC15+n8W+D9JY4B/RcTcAr/Y+gFHVryIiE8lHQh0BZ5N7dYgk6zl6gjMy3P+emBCRFTZ7CBpC2Brls0ePSqpT9ayW0Ucj0o6HLgO6JHTzUfAhsAnOW2GA8Mhs6k679Wa1bNmzZrRp08f+vTpA2SSkf3335+LL76Y0047rcosS7Zp06Zx1FFHMWjQIM4+++yl5+fPn8+SJUu45JJLqp0NKrRp+5tvvlmhWZjs2ZPsJ+MqZlMqyrO98MILHHTQQfzoRz/i5ptvrlR255138uyzz1JeXs5nn30GwJdffgnAV199xYIFC6p9UmtFrbvuunz66adVEqNPP/202msZN24cAwcOZNCgQZU2YgNcddVVzJ8/nzPPPHPptXz99ddA5q0AvvjiC9Zaq9KkPK1atfIeIlupNFRC9G3W8RIye2cAFrNsVqplgTbfZ73+nspx5v5Rj4h4R9KHkn4E7MSy2aLsSkMl3Q/8GJgkqR/wTTXxK884Ah6NiJqWphaSc22Sfg+sB5xUTZtDgEkR8WWq/yCwM1ApIZK0GpnEaSGwLjA3q7glhWe9zIpmww035IQTTuCXv/wlr7/+Or17985b78MPP+Sggw5i5513rrL00q5dO1ZbbTVOO+00jj322LztC22Obt++/dI/+MujYn/NzJkzKyVEFfttcmdGZsyYwX777UfPnj25++67WX31ys9TvPzyyyxcuHBpv9kOPvhg2rZtuzS5qG/bbLMN3377LW+88UalfUTVXcvjjz/O4YcfziGHHMKNN95Ypb+XX36ZDz74gE6dOlUp23777enRowdTp06tdH7+/Pl06NChSn2zYmnsD3edA+xAZi/NYcvZxxHAk5J2J7OsVbEOdBOZpbO/R8SS3EaSfhgRM4AZknYhs9F5am695BEyy1NnpbbrAJOA6yRtERGzJbUGNoqI3F2js4AtgPGp7QnAfsDeaXkvn/8CJ0q6jEzitSeZ2a5cg1P/5wMjJe0SEd+lJ9I2IPP9NSuqd955h4033rjK+VdeeQWg2ifQvvnmG/r370+bNm345z//WeXppjZt2rDHHnswbdo0tt9++zo/vr/VVlvx3XffMXfuXDbaKN9WvsJ22WUXOnTowJgxYyo9Ej969GjWXXdddtttt6XnXn/9dfbZZx8233xz7rvvPlq1alWlv4EDB1Z5HH7q1KkMHjyYK6+8kp122qnOMdbW/vvvzxprrMGYMWMqvSfT6NGj6datG5ttttnScxMnTqR///7svffejB49Ou/3fciQIVXepuChhx7i8ssvZ/To0Wy55ZaVyj744AO++eabKufNiqmxE6IrgTslHQM8UVPlanwq6T/A2mSWxyqMJbNUlm+5DOAsSXuRmbF6GXiQzPJWPn8kk/y8lOpfHBH/kjQQuF1Sxbz774DchOh+MjNBN6XXw4C3gYlpavpfEfGHtLn85Ig4AbgL+BGZ5cEAHoqIcdmdSupC5gm13hHxhaQJafzfk0kyJ0XE4mqux6zRdOvWjb322otDDjmEzTbbjM8//5wHHniAYcOG8T//8z9LnxjLddZZZzFlyhRGjRq1NHmq0LVrV9Zee23+7//+jz59+rDffvsxaNAgOnbsyMcff8yUKVNYsmQJQ4cOrTauiuW7559/vkpC9NRTTzFv3jw++OADAMrLy1lzzcwzE4cdlvm/2+qrr84ll1zCqaeeSqdOnejXrx9PPPEEI0eO5Nprr126N+mjjz5in332YdGiRVx88cVVntjabrvtaNGiBWVlZZSVleWNtUePHuy+++5LX48fP5699tor7/sj5XrwwQf56quvmDFjxtJr+/jjj2nTpg0HHHAAAOuvvz6DBw/msssuY6211mL77bfnjjvu4IknnuDf//730r5eeeUVfvKTn9ChQwd+85vfVHmSb+eddwYyyWbuE3QVb6ew0047VXma7bnnngOW/UzMVgZN6o0ZJY0Hzo6I8jxlvYC/RMTyPVNbjyQ9AxwYEQ0z3111vKuBsRHxeKF6fmPGhuM3Zlxm2LBhPPDAA0ybNo0PP/yQZs2a0aVLFwYMGMBZZ521NHEYNWoUP//5z3nrrbcoKyujb9++PPXUU3n7zH5zwVmzZnHxxRfzxBNPsGDBAtZbbz223357Tj75ZH784x8XjG2nnXaia9eu3HJL5f83FRo793fkjTfeyJ///GfefvttNtlkEwYPHsyppy5714uK5KU6FdebT0Xb3DdmvP/++znwwAN58MEH2X///QteY1lZGW+//XaV85tuummlj0lZsmQJl112GSNGjOCDDz5gyy235MILL1yaAMKyn1F1Cv39qGib740ZTzzxRF588UXKy6v8Kjerd6rlGzOuEgmRpCHAKcDREfFMMWLLiWcnYGFETG+k8U6seOKsECdEDccJUdMwatQofvnLX/L+++/nffPEldX555/P2LFjmTFjxgp/vEexffPNN3Ts2JErr7ySQYMGFTscKwG1TYia1BszRkTffLNDETE0IjZdGZIhgIh4rrGSoTRejcmQmWXeVbpTp05N7iMjnnrqKc4///wmnwxBZoZt/fXX57jjjit2KGaVNPYeIjOzomnWrBkjR45kypQpV6y2UgAAEQ9JREFUxQ6lTp599tlih1BvWrRowahRo6psmjcrtia1ZGYrxktmDcdLZmZmK6faLpk5RS8h3Tu1pdx/uM3MzKpoUnuIzMzMzBqCEyIzMzMreU6IzMzMrOQ5ITIzM7OS503VJWTGuwsoG3J/scOwEucn8sxsZeQZIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3lOiOpAUl9J96XjgZL+Wk29gyVdmI43lfS4pOmSxkvaqJo2l0p6R9KXyxHX6ZJ+Xtd2ZmZmluGEKIek+ngrgnOA69PxlcDfImJb4A/AZdW0GQf0Xs7xRgJnLmdbMzOzkteoCZGkMkmzJI2QNFPSI5JapbLxknql4w6S5qTjgZLulTRO0ltpNuRXkl6UNEnSulntr5L0H0kvSeotaTVJr0taL9VZTdJsSR1y4rpI0nBJjwB/k9RS0i2SZqRx9qrDNXYBvo2Ij9OprsDj6fhJoH++dhExKSLer6Hva7JmnvaTNEHSahHxNTBH0vImVGZmZiWtGDNEnYHrImIb4DPg0Fq06QYcRWYG5VLg64jYDpgIHJtVr01E7AqcCoyMiO+B0cDRqbwfMC0rWcm2A9A/Io4CTgOIiO7AAOBWSS1reX27AVOyXk9j2TUeAqwlqX0t+8o1BDgiJWjXAD9P1whQDuyR20DSLySVSypf8vWC5RzWzMxs1VaMhOitiJiajicDZbVo82REfBER84AFZJaXAGbktL8dICImAGtLakdmOakiaToeuKWaMcZGxMJ0vDvw99TXK8DbQJdaxAnQEZiX9fpsYE9JLwJ7Au8Ci2vZVyVpJuhE4FHgrxHxRlbxR8CGedoMj4heEdGrWeu2yzOsmZnZKq8YH93xbdbxEqBVOl7MsgQtdzYmu833Wa+/p/I1RE67iIh3JH0o6UfATiybLcr1Vdaxqg+/RguBpZlHRLwH/BRA0prAoRGxIlM13YFPqJr8tExjm5mZWR2tTJuq55BZtgI4bDn7OAJA0u7AgqzE4yYyS2d3RsSSWvQzgZQ4pT1BmwCv1jKGWcAWFS/SfqiK7/N5ZGasloukTYFfA9sBB0jaKau4C/DS8vZtZmZWylamhOhK4BRJ/wE61FS5Gp+m9sOAQVnnxwJrUv1yWa7rgWaSZgB3AAMj4tsa2lSYAGwnqWKWqS/wqqTXgB+Q2QMFgKSpWcd/kjQXaC1prqSLsjtN/d0MnJ1mnQYBN2XtbdoNeKyWMZqZmVkWReSuMjVNksaTSRbK85T1Av4SEVU2HTdQLFcD4yKiURIUSdsBv4qIYwrVa9Gxc3Q87qrGCMmsWv60ezNrTJImR0SvmuqtTDNEDULSEOBuMstVjeV/gdaNOF4H4IJGHM/MzGyVUoxN1Q0iIvpWc34oMLSRY/mQzDJdY433aGONZWZmtipa5WeIzMzMzGrihMjMzMxK3iqzZGY1696pLeXe0GpmZlaFZ4jMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3neVF1CZry7gLIh9xc7DDMzs7yK+U72niEyMzOzkueEyMzMzEqeEyIzMzMreU6IzMzMrOQ5ITIzM7OSt9IkRJKWSJoq6SVJ/5TUukDdMklHZb0eKOmvDRTXwZIuzBpnXopzqqQT8sQ/VdLYrPOnS5otKSR1KDDOJpIekTRL0suSyuoQ4z8kdV6+KzQzM7MGTYgk1eWx/oUR0TMiugGLgJML1C0DjipQXp/OAa7Pen1HirNnRNyUdX5h1vmDss4/C/QD3q5hnL8BV0TE1kBv4KM6xHhDitPMzMyWQ8GEKM3EvCLpVknTJd1VMXMjaU7FjIekXpLGp+OLJA2X9AjwN0nNJF0h6YXUx0m1iOtpYAtJl0j6ZVY8l0o6ExgK7JFmYwan4g0lPSTpdUl/ymozQNKMNPN0edb5L1N/0yRNkvSDPNffBfg2Ij6uRcx5RcSLETGnUB1JXYHmEfFoavNlRHydU6d5+h72Ta8vk3RpKn4a6FfHBNTMzMyS2swQbQkMj4htgc+BU2vRZgegf0QcBQwCFkTEjsCOwImSNquuYfqjfgAwA7gZOC6dXw04EhgDDAGeTrMxf0lNewJHAN2BIyRtLGlD4HLgR6l8R0kHp/ptgEkR0QOYAJyYJ5zdgCk55w7NSg43zjrfUlJ5Sq4Opm66AJ9J+pekF1MC2Sy7QkQsBgYCN0jaB9gfuDiVfQ/MBnrkdizpFymu8iVfL6hjWGZmZqWhNgnROxHxbDoeDexeizZjI2JhOt4XOFbSVOA5oD2Qb79Lq1SnHPgvcHOaWflE0napnxcj4pNqxnw8IhZExDfAy8CmZBKw8RExLyUUY4A+qf4i4L50PJnMMlyujsC8rNfjgLKUHD4G3JpVtklE9CKzlHeVpB9WE2c+zYE9gLNTzJuTSX4qiYiZwN9THMdHxKKs4o+ADfO0GR4RvSKiV7PWbesQkpmZWemozRJLVPN6McsSqpY5db7KOhZwRkQ8XMM4CyOiZ57zN5FJDjYARhZo/23W8RIy16YC9b+LiMipXyUmYGkWkZOMjSAz+1RR9l769820fLgd8EaB8bPNJZPsvQkg6V5gZzIzZLm6A58BuUt8LVO8ZmZmVke1mSHaRNIu6XgA8Ew6nkNmaQzg0ALtHwZOkbQ6ZPblSGpThxjvIbM8tGPqC+ALYK1atH0O2FNSh7QENQB4qg5jzwK2qHghqWNW2UGpHEnrSGqRjjuQWWp7uQ7jvACsI2m99PpH+dpL+imZGbY+wDWS2mUVdwFm1mFMMzMzS2qTEM0CjpM0HViXzBNNkNm/crWkp8nMsFTnJjJ/3KdIegm4kTp8qGxaFnoSuDMiKsaZDixOG6IHF2j7PnBeaj8NmBIR/67t2GT2Fm0nqWKm6UxJMyVNA85k2bLW1kB5Ov8kMDQiXgaQdKakucBGwHRJN6XzvSqO03WdDTwuaQaZma0R2YGkRGsoMCgiXgP+Clydyn5AZobt/Tpcm5mZmSVatmqUpzDzXjj3pUfhiyJtpp4CHB4Rrxdh/KuBcRHxWGOPXVspKfw8IvItsS3VomPn6HjcVY0UlZmZWd00xKfdS5qc9vgWtNK8MWM+6XH02WQ2TDd6MpT8L1Dtm0SuJD6j8gZvMzMzq4OCS1fpKa+izQ6lZafNizV+iuFDYGyNFYsoIm4pdgxmZmZN2Uo9Q2RmZmbWGJwQmZmZWcnzRz2UkO6d2lLeABvWzMzMmjrPEJmZmVnJc0JkZmZmJc8JkZmZmZU8J0RmZmZW8pwQmZmZWckr+NEdtmqRNA94ux66agssqId+VqSvurarbf2a6i1veQfg41qMvzKpz59zY41TKvcTNL17yvfTitf3/bRMXb7Hm0bEejXWigh/+atOX8DwYvdV13a1rV9TveUtB8qL/XMr5s+5scYplfsplTWpe8r304rX9/3UMD/nii8vmdnyGLcS9FXXdrWtX1O9FS1vShrrWnw/rfg4TYHvpxWv7/tpmXq/Fi+ZmTUCSeVRi09bNqst31NWn3w/eVO1WWMZXuwAbJXje8rqU8nfT54hMjMzs5LnGSIzMzMreU6IzMzMrOQ5ITIzM7OS54TIrMgkHSxphKR/S9q32PFY0yZpc0k3S7qr2LFY0ySpjaRb0++lo4sdT2NxQmS2AiSNlPSRpJdyzu8v6VVJsyUNKdRHRNwbEScCA4EjGjBcW8nV0/30ZkQMathIramp4731U+Cu9HvpoEYPtkicEJmtmFHA/tknJDUDrgMOALoCAyR1ldRd0n05X+tnNf1damelaxT1dz+ZZRtFLe8tYCPgnVRtSSPGWFTNix2AWVMWERMkleWc7g3Mjog3AST9A+gfEZcBB+b2IUnAUODBiJjSsBHbyqw+7iezfOpybwFzySRFUymhiZOSuVCzRtSJZf+7gswvl04F6p8B9AMOk3RyQwZmTVKd7idJ7SUNA7aTdF5DB2dNWnX31r+AQyXdwKr1cR8FeYbIrP4pz7lq3wE1Iq4Brmm4cKyJq+v99AngxNpqI++9FRFfAT9v7GCKzTNEZvVvLrBx1uuNgPeKFIs1fb6frKH43srihMis/r0AdJa0maQ1gCOBsUWOyZou30/WUHxvZXFCZLYCJN0OTAS2lDRX0qCIWAycDjwMzALujIiZxYzTmgbfT9ZQfG/VzB/uamZmZiXPM0RmZmZW8pwQmZmZWclzQmRmZmYlzwmRmZmZlTwnRGZmZlbynBCZmZlZyXNCZGYNStISSVMlvSRpnKR2K9BXX0m7Zr0+WdKx9RNprcbfQ9LMdD2talH/pvTp4cszVpmkl+rY5j9ZbY9annEL9H1+vrHMVhV+HyIza1CSvoyINdPxrcBrEXHpcvZ1EfBlRFxZjyHWZfxhwHMRcUsjjFUG3BcR3ZajbV/g7Ig4sA5tmkXEkgLlS3+OZqsizxCZWWOaSPqk9jTbc19FgaS/ShqYjudIuljSFEkzJG2VEoSTgcFphmYPSRdJOju1GS/pL5ImSJolaUdJ/5L0uqQ/Zo3zM0nPpz5ulNQsN0hJe0t6MY09UlILSScA/wNcKGlMTv02ku6XNC3NhB2RFVOvdPylpEtTnUmSfpDO/zC9fkHSHyR9mSeeZpKuSHWmSzop3zc3q+1QYI90jYOra59+Bk9Kug2Ykc7dK2lymgn7RTo3FGiV+huTPZYyrkjXPSPr2vum679L0iuSxkjK92GiZisFJ0Rm1ihS4rE3tf+spI8jYnvgBjKzHXOAYcBfIqJnRDydp82iiOiT6v0bOA3oBgyU1F7S1sARwG4R0RNYAhydE2dLYBRwRER0B5oDp0TETSn230REpTbA/sB7EdEjzeg8lCe2NsCkiOgBTABOTOevBq6OiB2p/oM1BwELUp0dgRMlbVZNXYAhwNPp+/SXGtr3Bn4bERVLe8dHxA5AL+BMSe0jYgiwMPWXe+0/BXoCPYB+wBWSOqay7YCzgK7A5sBuBWI2KyonRGbW0FpJmgp8AqwLPFrLdv9K/04GymrZpiLZmgHMjIj3I+Jb4E0yn+q9N7AD8EKKaW8yf6izbQm8FRGvpde3An1qGHcG0E/S5ZL2iIgFeeosAipmxLKvaRfgn+n4tmr63xc4NsX8HNAe6FxDTLVt/3xEvJVV90xJ04BJZL5nNY2zO3B7RCyJiA+Bp8gkXRV9z42I74Gp1P7naNbomhc7ADNb5S2MiJ6S2pJJCE4DrgEWU/k/ZS1z2n2b/l1C7X9XVbT5Puu44nVzQMCtEXFegT7qvKwTEa9J2gH4MXCZpEci4g851b6LZZs263JNFTGdEREP1zW2Qu3TXqOvcl73A3aJiK8ljafqzyVf39XJ/hnU9ZrNGpVniMysUaRZkzOBsyWtDrwNdE37c9qSma2pyRfAWisQxuPAYZLWB5C0rqRNc+q8ApRJ2iK9PobMrEe1JG0IfB0Ro4Erge3rENMk4NB0fGQ1dR4GTknfNyR1kdSmQJ+536fatm8LfJqSoa2AnbPKvqton2MCcETap7Qemdm05wvEZrZSckJkZo0mIl4EpgFHRsQ7wJ3AdGAM8GItuhgHHFKxqXo5xn8Z+B3wiKTpZJbvOubU+Qb4OfBPSTPIzC4Nq6Hr7sDzaUnqt8Afa6if7SzgV5KeT7HkW267CXgZmKLMo/g3Uni2ZTqwOG3gHlyH9g8BzdP35hIyyVqF4cD03A3lwD1pvGnAE8A5EfFBoQs2Wxn5sXszsyKS1JrMsmJIOhIYEBH9ix2XWanxeq6ZWXHtAPw1PZL+GXB8keMxK0meITIzM7OS5z1EZmZmVvKcEJmZmVnJc0JkZmZmJc8JkZmZmZU8J0RmZmZW8pwQmZmZWcn7/4EI+iOdz5+1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def bar_plot(*shape):\n",
    "    names = tuple(all_implementations.keys())\n",
    "    runtimes = tuple(do_benchmark(all_implementations[name], shape) for name in names)\n",
    "    for runtime, name in zip(runtimes, names):\n",
    "        assert runtime >= runtimes[names.index('pystencils')], runtimes\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(16, 16)\n",
    "\n",
    "plt.subplot(3, 1, 2)\n",
    "bar_plot(128, 128)\n",
    "\n",
    "plt.subplot(3, 1, 3)\n",
    "bar_plot(1024, 1024)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All runtimes are plotted logarithmically. Next number next to the labels shows how much slower the version is than the fastest one. For small arrays Cython produces faster code than *pystencils*. The larger the arrays, the better pystencils gets."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}