{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ba44776-ef27-4684-9ce4-1a34ce2e3d57",
   "metadata": {
    "tags": []
   },
   "source": [
    "## **Notebook 3: Introduction to Numba**\n",
    "\n",
    "This tutorial introduces Numba basics and common usages before diving into our main nearest neighbors example. Numba is a just-in-time compiler optimized for numerically-focused Python functions that can work on the CPU or the GPU. Numba can compile a large number of NumPy functions and can create ufuncs (universal functions), parallelize loops, and generate GPU-accelerated code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e24a3ca0-70a1-49f5-9880-836e0216b733",
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "libjpeg.so.62: cannot open shared object file: No such file or directory",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_19282/29068482.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcupy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mcp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/conda/envs/rapids/lib/python3.9/site-packages/matplotlib/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    111\u001b[0m \u001b[0;31m# cbook must import matplotlib only within function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    112\u001b[0m \u001b[0;31m# definitions, so it is safe to import from it here.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 113\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0m_api\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_version\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_docstring\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrcsetup\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    114\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msanitize_sequence\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    115\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_api\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mMatplotlibDeprecationWarning\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/envs/rapids/lib/python3.9/site-packages/matplotlib/rcsetup.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0m_api\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcbook\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mls_mapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolors\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mColormap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_color_like\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     28\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fontconfig_pattern\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mparse_fontconfig_pattern\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_enums\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mJoinStyle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mCapStyle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/envs/rapids/lib/python3.9/site-packages/matplotlib/colors.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumbers\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mNumber\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 51\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mPIL\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mImage\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     52\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mPIL\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPngImagePlugin\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mPngInfo\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.9/site-packages/PIL/Image.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     98\u001b[0m     \u001b[0;31m# Also note that Image.core is not a publicly documented interface,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     99\u001b[0m     \u001b[0;31m# and should be considered private and subject to change.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 100\u001b[0;31m     \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0m_imaging\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    102\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0m__version__\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcore\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"PILLOW_VERSION\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: libjpeg.so.62: cannot open shared object file: No such file or directory"
     ]
    }
   ],
   "source": [
    "from numba import (cuda,  \n",
    "                   float32,\n",
    "                   jit, \n",
    "                   vectorize,\n",
    "                   guvectorize)\n",
    "\n",
    "import numpy as np\n",
    "import cupy as cp\n",
    "import time\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecb4665a-ad8c-4033-9169-d1912cb1f2e1",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Numba Vectorize/Guvectorize**\n",
    "\n",
    "The Numba `vectorize` decorator compiles a Python function that takes a scalar input into a [NumPy ufuncs](https://numpy.org/doc/stable/reference/ufuncs.html). These ufuncs or universal functions work elementwise over NumPy ndarrays. This is an example of a dynamic unfunc (DUFunc), as we are not passing in any signatures to our vectorize decorator. When a different input type is passed, we will dynamically compile a new kernel.\n",
    "\n",
    "We are running the relative difference ufunc below on the CPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0b1682a-8177-4736-8497-df96faabd588",
   "metadata": {},
   "outputs": [],
   "source": [
    "@vectorize()\n",
    "def rel_diff_cpu(x, y):\n",
    "    return 2 * (x - y) / (x + y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b963224d-3bf2-453a-a041-17d5fa1d11f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.56.2\n",
      "1.21.6\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numba\n",
    "import numpy\n",
    "\n",
    "@numba.vectorize(nopython=True, cache=True)\n",
    "def v_floor(a):\n",
    "    return numpy.floor(a)\n",
    "\n",
    "print(numba.__version__)\n",
    "print(numpy.__version__)\n",
    "\n",
    "v_floor(1.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4adb740-eb53-4144-b0bf-6feaaa0359fb",
   "metadata": {
    "tags": []
   },
   "source": [
    "In the next cell we run our relative difference function on increasing NumPy array sizes. We are timing a NumPy ufunc implementation and our Numba ufunc implementation. Note we are ignoring divide by zero errors we might get with this method of random number generation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af2eaf45-0004-4f3d-8aee-3361181d3201",
   "metadata": {},
   "outputs": [],
   "source": [
    "size_list = [10, 10000, 100000, 1000000, 10000000, 100000000]\n",
    "numpy_cpu_times = []\n",
    "numba_cpu_times = []\n",
    "\n",
    "np.seterr(divide='ignore')\n",
    "\n",
    "for size in size_list:\n",
    "    x=np.random.randn(size)#.astype(np.float32)\n",
    "    y=np.random.randn(size)#.astype(np.float32)\n",
    "\n",
    "    start_time_numpy = time.monotonic()\n",
    "    2 * (x - y) / (x + y)\n",
    "    numpy_cpu_times += [(time.monotonic() - start_time_numpy)]\n",
    "\n",
    "    start_time_numba_cpu = time.monotonic()\n",
    "    rel_diff_cpu(x, y)\n",
    "    numba_cpu_times += [(time.monotonic() - start_time_numba_cpu)]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5457736e-289b-4198-8e18-9d7590ac4f7e",
   "metadata": {},
   "source": [
    "Next, lets plot these timings to show how we see a bigger performance benefit as we move to larger input arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1210e20b-f3c4-498b-8a55-81fc4ee6b8aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkSklEQVR4nO3dd3xUVfrH8c+kF1JIIfSe0DuKIARUiqBIEcXVVWkqq64iq66oPxF1QWRBFCkWirroojTLIoKN0JWmIiihSCgBkgBJSCBl5v7+uDA4JkACk9xk8n2/XvNaz713Zp65m2QeznnOOTbDMAxEREREPISX1QGIiIiIuJOSGxEREfEoSm5ERETEoyi5EREREY+i5EZEREQ8ipIbERER8ShKbkRERMSj+FgdQGlzOBwcPnyYkJAQbDab1eGIiIhIERiGQWZmJtWrV8fL6+J9MxUuuTl8+DC1atWyOgwRERG5DAcOHKBmzZoXvabCJTchISGAeXNCQ0MtjkZERESKIiMjg1q1ajm/xy+mwiU354aiQkNDldyIiIiUM0UpKVFBsYiIiHgUJTciIiLiUZTciIiIiEepcDU3RWW328nLy7M6DLGYn5/fJacciohI2aLk5k8Mw+DIkSOcPHnS6lCkDPDy8qJevXr4+flZHYqIiBSRkps/OZfYVKlShaCgIC30V4GdW/AxOTmZ2rVr62dBRKScUHLzB3a73ZnYREZGWh2OlAHR0dEcPnyY/Px8fH19rQ5HRESKQMUEf3CuxiYoKMjiSKSsODccZbfbLY5ERESKSslNITT8IOfoZ0FEpPxRciMiIiIeRcmNiIiIeBQlN+IWv//+OzabjW3btlkdioiIVHBKbjzEkCFDsNlsvPzyyy7Hly5d6hF1I4sWLaJbt26EhYVRqVIlWrZsyQsvvMDx48cBmDdvHjabzfmoVq0at99+O/v27XO+hs1mY+nSpQVee9SoUXTr1q2UPomIiIdL+Q3S9lgagpIbDxIQEMDEiRM5ceKE1aG41TPPPMPgwYO56qqr+OKLL9i+fTuTJ0/mxx9/5P3333deFxoaSnJyMocPH+aDDz5g27Zt3HLLLZrpJCJSWrZ9AG91g4/uhbwzloWh5OYSDMMgOzffkodhGMWKtXv37lStWpUJEyZc8Jrnn3+e1q1buxybOnUqdevWdbaHDBlC//79GT9+PDExMYSHhzNu3Djy8/N54okniIiIoGbNmsyZM6fA6//666906tSJgIAAmjVrxnfffec8Z7fbGT58OPXq1SMwMJBGjRrx2muvXfQzff/994wfP57JkyczadIkOnXqRN26denRoweLFi3i3nvvdV5rs9moWrUq1apV47rrrmPs2LFs376d3bt3X/zGiYjIlcnNgiV/g6V/g7xsCKps/q9FtIjfJZzOs9P0uS8tee8dL/QiyK/o/xd5e3szfvx47rzzTh555BFq1qx52e/9zTffULNmTRISEli7di3Dhw9n/fr1xMfHs3HjRhYsWMDIkSPp0aMHtWrVcj7viSeeYOrUqTRt2pQpU6Zwyy23sG/fPiIjI3E4HNSsWZOPPvqIqKgo1q1bx/333+8cQirM/PnzqVSpEg8++GCh58PDwy/4GQIDAwG0R5iISEk6+gt8PARSd4HNC7qNgS7/AC9vy0JSz42HGTBgAK1bt2bs2LFX9DoRERG8/vrrNGrUiGHDhtGoUSOys7N5+umniY2NZcyYMfj5+bF27VqX5z388MPceuutNGnShJkzZxIWFsbs2bMB8PX1Zdy4cVx11VXUq1ePu+66iyFDhvDRRx9dMI7ExETq169f7NWBDx48yKRJk6hZsyZxcXHFvwEiInJxhgGb34W3rzcTm5BqcO9n0PVJSxMbUM/NJQX6erPjhV6WvfflmDhxItdffz3/+Mc/Lvu9mzVr5rIbdkxMDM2bN3e2vb29iYyM5NixYy7P69ixo/O/fXx8aN++PTt37nQemzVrFu+88w779+/n9OnT5ObmFhgm+yPDMIpcEJ2enk6lSpXMocTsbNq2bcvixYu16aWIiLvlZMJno2D7QrPdsDsMeBOCoywN6xwlN5dgs9mKNTRUFsTHx9OrVy+efvpphgwZ4nLOy8urQC1PYcM2f+4psdlshR5zOByXjOdccvLRRx/x2GOPMXnyZDp27EhISAiTJk1i48aNF3xuXFwca9asIS8v75K9NyEhIWzZsgUvLy9iYmIIDg4ucD49Pb3A806ePElYWNglP4eIiADJP5rDUMf3gs0bbvg/6PQoeJWdwaCyE4m41csvv8xnn33GunXrXI5HR0dz5MgRlwTHnWvTbNiwwfnf+fn5bN68mcaNGwOwevVqOnXqxIMPPkibNm1o2LAhe/ZcfLrgnXfeyalTp5gxY0ah50+ePOn8by8vLxo2bEj9+vULJDYAjRs35ocffnA5ZhgGmzdvplGjRkX9iCIiFZNhwPdvwzs9zMQmtAYMXQadHytTiQ2o58ZjtWjRgrvuuotp06a5HO/WrRspKSm88sorDBo0iOXLl/PFF18QGhrqlvedPn06sbGxNGnShFdffZUTJ04wbNgwABo2bMh7773Hl19+Sb169Xj//ff54YcfqFev3gVfr0OHDjz55JP84x//4NChQwwYMIDq1auze/duZs2aRefOnXn00UeLFNvjjz/OvffeS+PGjenZsyenT5/mrbfeYs+ePTz00ENu+fwiIh7pTDp8+nfY8YnZjrsR+s+EoAhr47qAspVqiVu9+OKLBYagmjRpwowZM5g+fTqtWrXi+++/5/HHH3fbe7788stMnDiRVq1asXr1aj755BOioswx2JEjRzJw4EAGDx5Mhw4dSEtLu+AsqD+aOHEiH3zwARs3bqRXr140a9aM0aNH07JlS5ep4Jdy++23M2/ePN59912uuuoqevbsyZ49e1i9ejV16tS57M8sIuLRDm2GWV3MxMbLB3r+C/7y3zKb2ADYjOIuplLOZWRkEBYWRnp6eoHeijNnzrBv3z7q1atHQECARRFKWaKfCRGpsAwDNsyElc+BIw/Ca8OgeVCznSXhXOz7+88s7blJSEigb9++VK9e/YJL4//ZqlWraNeuHQEBAdSvX59Zs2aVfKAiIiIVSfZx+O9d8OUYM7Fp0hceWG1ZYlNcliY3WVlZtGrVijfeeKNI1+/bt48+ffrQpUsXtm7dytNPP80jjzzCokWLSjhSERGRCuLAD/BmPPz2P/D2g96T4Pb3ITDc6siKzNKC4t69e9O7d+8iXz9r1ixq167N1KlTAbN+ZNOmTfz73//m1ltvLfQ5OTk55OTkONsZGRlXFLOIiIhHcjhg/TT4+gVw5EPlenDbXKjexurIiq1cFRSvX7+enj17uhzr1asXmzZtuuAS+xMmTCAsLMz5+ONWASIiIgJkpcGHg8/W1+RDs4HwQEK5TGygnCU3R44cISYmxuVYTEwM+fn5pKamFvqcMWPGkJ6e7nwcOHCgNEIVEREpH/avg1mdIXEFePvDzVNh0BwIcM8SIVYod+vc/Hkp/nOTvS60RL+/vz/+/v4lHpeIiEi54nDAmsnw7XgwHBAZC7fNg6rNL/nUsq5cJTdVq1blyJEjLseOHTuGj48PkZGRFkUlIiJSzpw6Bovvh73fmu2Wd8BNk8G/krVxuUm5Sm46duzIZ5995nJsxYoVtG/fvti7RouIiFRIe1fB4vvg1FHwCYSb/g2t74IiblJcHlhac3Pq1Cm2bdvm3Nto3759bNu2jaSkJMCsl7nnnnuc148cOZL9+/czevRodu7cyZw5c5g9e7ZbV9gVERHxSA47fDsB3utnJjbRjeH+b6HNXz0qsQGLk5tNmzbRpk0b2rQxq7FHjx5NmzZteO655wBITk52JjoA9erVY9myZXz33Xe0bt2aF198kddff/2C08ArkiFDhmCz2Xj55Zddji9duvSC9UjuZLPZnI+QkBDat2/P4sWLS/x9RUSkCDKSzaRm1cuAYSY0930LVZpYHVmJsHRYqlu3bgX2PvqjefPmFTjWtWtXtmzZUoJRlV8BAQFMnDiRBx54gMqVK5f6+8+dO5cbb7yRkydPMmnSJG677TbWrFlDx44dSz0WERE5a/fXZn1Ndir4BsPNr0KrwVZHVaLK1VRwubju3btTtWpVJkyYUOj5559/ntatW7scmzp1KnXr1nW2hwwZQv/+/Rk/fjwxMTGEh4czbtw48vPzeeKJJ4iIiKBmzZrMmTOnwOuHh4dTtWpVGjduzKxZswgICODTTz8lISEBX1/fAsXg//jHP4iPj7/izy0iIoWw58NX4+A/A83EJqY5PLDK4xMbUHJzaYYBuVnWPIq5p6m3tzfjx49n2rRpHDx48LI/8jfffMPhw4dJSEhgypQpPP/889x8881UrlyZjRs3MnLkSEaOHHnRNYN8fX3x8fEhLy+P+Ph46tevz/vvv+88n5+fz3/+8x+GDh162XGKiMgFpB+Cd2+GNVPMdvthMOIriIq1Nq5SUq5mS1kiLxvGV7fmvZ8+DH7BxXrKgAEDaN26NWPHjmX27NmX9bYRERG8/vrreHl50ahRI1555RWys7N5+umnAbPQ++WXX2bt2rXccccdBZ6fk5PDpEmTyMjI4IYbbgBg+PDhzJ07lyeeeAKA//3vf2RnZ3P77bdfVowiInIBu76EJSPh9HHwC4FbXoPmFas2VT03HmjixIm8++677Nix47Ke36xZM7y8zv9oxMTE0KJFC2fb29ubyMhIjh075vK8v/zlL1SqVImgoCCmTJnCv//9b+feYUOGDGH37t1s2LABgDlz5nD77bcTHFy85E1ERC7AngcrnoUPbjcTm2qtzGGoCpbYgHpuLs03yOxBseq9L0N8fDy9evXi6aefZsiQIc7jXl5eBQq4C9uT689rBtlstkKPORwOl2Ovvvoq3bt3JzQ0lCpVqricq1KlCn379mXu3LnUr1/fOetNRETc4MR+WDgMDm0y21c/AD1fBJ+KuUK/kptLsdmKPTRUFrz88su0bt2auLg457Ho6GiOHDmCYRjO6eHn1hhyh6pVq9KwYcMLnh8xYgR33HEHNWvWpEGDBlx77bVue28RkQpr5+fwyYNwJh0CwqDfdGjS1+qoLKVhKQ/VokUL7rrrLqZNm+Y81q1bN1JSUnjllVfYs2cP06dP54svvii1mHr16kVYWBgvvfSSColFRK5Ufg588U9YcJeZ2NRoBw+srvCJDSi58WgvvviiyzBUkyZNmDFjBtOnT6dVq1Z8//33pbq6s5eXF0OGDMFut7usPC0iIsV0fC/M7gkbZ5ntjg/D0OVQuY61cZURNuNiq+h5oIyMDMLCwkhPTyc01HU79zNnzrBv3z7q1atHQECARRF6tvvuu4+jR4/y6aefWh1KkehnQkTKnF+WwKePQE4GBFaG/jOhUW+roypxF/v+/jPV3EipSE9P54cffmD+/Pl88sknVocjIlL+5J2BL5+GTWeX+ah1DQyaDWE1rY2rDFJyI6WiX79+fP/99zzwwAP06NHD6nBERMqX1N3w8RA4+rPZ7jwarnsavH0v+rSKSsmNlApN+xYRuUw/fQSfjYK8LAiKgoFvQsPuVkdVpim5ERERKYtys+GLJ2Hr2a1r6nSGW9+B0GrWxlUOKLkpRAWrsZaL0M+CiFji2K/mMFTKTsAGXZ+E+CfBW1/bRaG79AfnVuHNzs4mMDDQ4mikLMjNzQXMLSdERErF1vmw7HFzb8PgKnDr21C/m9VRlStKbv7A29ub8PBw555JQUFBzpV8peJxOBykpKQQFBSEj49+VUSkhOWcMpOaHz802/W7wcC3oVKViz5NCtJf7D+pWrUqQIFNIaVi8vLyonbt2kpyRaRkHdkOC4dC6i6weZkzoTqPBi/1Gl8OJTd/YrPZqFatGlWqVCl0U0mpWPz8/Fx2SBcRcSvDgM3zYPlTkH8GQqrBrbOhrvbeuxJKbi7A29tbdRYiIlJyzmTA56Ng+yKz3bA7DHgTgqMsDcsTKLkREREpbck/mrOhju8Fmzfc8Bx0egTUU+wWSm5ERERKi2HAD++Y2yjYcyG0JgyaA7U7WB2ZR1FyIyIiUhpOn4RP/w47z24c3KgP9JsOQRGWhuWJlNyIiIiUtEOb4eOhcHI/ePlCjxfgmr+BZmKWCCU3IiIiJcUwYMNMWPkcOPIgvDbcNg9qtLM6Mo+m5EZERKQkZB+HTx6C35aZ7SZ94ZY3IDDc0rAqAiU3IiIi7pa0ERYOg4yD4O0HvcbDVSM0DFVKlNyIiIi4i8MB616Hr18Aww4R9c1hqGqtrI6sQlFyIyIi4g5ZqbBkJOxeabab3wo3T4WAUEvDqoiU3IiIiFyp39fCouGQmQw+AdB7IrS9V8NQFlFyIyIicrkcdlg9Bb4bD4YDImPNYaiqza2OrEJTciMiInI5Th2DxffB3u/Mdss74KbJ4F/J0rBEyY2IiEjx7f0OFt0HWcfANwj6/Bva3GV1VHKWkhsREZGicthh1URY9QpgQHQTcxiqSmOrI5M/UHIjIiJSFBnJsGgE7F9jttveAzdOBL8ga+OSApTciIiIXMrur2DxA5CdCn6VzCneLW+zOiq5ACU3IiIiF2LPh29fgjWvmu2YFuYwVFRDS8OSi1NyIyIiUpj0g7BwOBzYYLbbDze3UfANsDYuuSQlNyIiIn/223JYOhJOnwD/UOj7GjQfaHVUUkRKbkRERM7Jz4Wvx8H6N8x2tdZw21xzjygpN5TciIiIAJzYDwuHwqHNZrvD36DHOPDxtzYuKTYlNyIiIjs/g08egjPpEBAG/WZAk5utjkouk5IbERGpuPJzYMX/wfdvmu0a7WHQHKhcx9q45IoouRERkYopbY85DJX8o9nu9He4YSx4+1obl1wxJTciIlLxbF8Mnz4CuZkQGAEDZkFcL6ujEjdRciMiIhVH3mlYPgY2zzXbtTvCrbMhrIa1cYlbKbkREZGKITURPh4CR7eb7c6j4bpnwFtfhZ5G/4+KiIjn+3EBfP4Y5GVBUBQMfBMadrc6KikhSm5ERMRz5WbDF0/A1v+Y7bpdYODbEFrN2rikRCm5ERERz3RspzkMlfIrYIOu/4SuT4KXt9WRSQlTciMiIp7FMGDbfPjf45B/GirFmL019btaHZmUEiU3IiLiOXJOwf/+AT/912zXvw4GvgWVqlgbl5QqJTciIuIZjmw3h6HSEsHmZc6E6jwavLysjkxKmZIbEREp3wzDXLfmi6fAngMh1WHQbKjTyerIxCKWp7MzZsygXr16BAQE0K5dO1avXn3R6+fPn0+rVq0ICgqiWrVqDB06lLS0tFKKVkREypQzGbBwmDnN254DDXvAyDVKbCo4S5ObBQsWMGrUKJ555hm2bt1Kly5d6N27N0lJSYVev2bNGu655x6GDx/OL7/8wscff8wPP/zAiBEjSjlyERGx3OFt8GY8/LIYvHygxwtw50cQHGl1ZGIxS5ObKVOmMHz4cEaMGEGTJk2YOnUqtWrVYubMmYVev2HDBurWrcsjjzxCvXr16Ny5Mw888ACbNm0q5chFRMQyhgEb34LZPeDEPgirBUO/gGsfVX2NABYmN7m5uWzevJmePXu6HO/Zsyfr1q0r9DmdOnXi4MGDLFu2DMMwOHr0KAsXLuSmm2664Pvk5OSQkZHh8hARkXLq9En46G5zYT57LjTqAw8kQK2rrY5MyhDLkpvU1FTsdjsxMTEux2NiYjhy5Eihz+nUqRPz589n8ODB+Pn5UbVqVcLDw5k2bdoF32fChAmEhYU5H7Vq1XLr5xARkVJycDO82QV2fgZevnDjy3DHBxAUYXVkUsZY3n9ns9lc2oZhFDh2zo4dO3jkkUd47rnn2Lx5M8uXL2ffvn2MHDnygq8/ZswY0tPTnY8DBw64NX4RESlhhgHr3oA5PeFkEoTXgeFfwjV/gwt8X0jFZtlU8KioKLy9vQv00hw7dqxAb845EyZM4Nprr+WJJ54AoGXLlgQHB9OlSxdeeuklqlUruFeIv78//v7+7v8AIiJS8rKPw9K/wa7lZrtpP7hlGgSEWRuXlGmW9dz4+fnRrl07Vq5c6XJ85cqVdOpU+BS+7OxsvP5ULObtbe4RYhhGyQQqIiLWSNoIs7qYiY23P9w0GW57V4mNXJKli/iNHj2au+++m/bt29OxY0feeustkpKSnMNMY8aM4dChQ7z33nsA9O3bl/vuu4+ZM2fSq1cvkpOTGTVqFFdffTXVq1e38qOIiIi7OByw7jX4+kUw7BDRAG6bB9VaWh2ZlBOWJjeDBw8mLS2NF154geTkZJo3b86yZcuoU6cOAMnJyS5r3gwZMoTMzEzeeOMN/vGPfxAeHs7111/PxIkTrfoIIiLiTlmpsOQB2P2V2W4+CPpOBf8QS8OS8sVmVLDxnIyMDMLCwkhPTyc0NNTqcERE5Jzf18CiEZCZDD4B0PsVaHuPioYFKN73t/aWEhERaznssHoyfDcBDAdExZnDUDHNrI5MyiklNyIiYp3Mo7D4Pti3ymy3uhNu+jf4BVsbl5RrSm5ERMQae7+DRfdB1jHwDTJnQ7W+0+qoxAMouRERkdJlz4dVEyFhEmBAlabmMFR0I6sjEw+h5EZEREpPxmGzaHj/WrPd9h64cSL4BVkbl3gUJTciIlI6Er+CJfdDdhr4VYKbp0LL26yOSjyQkhsRESlZ9jz45iVYO9VsV20Bg+ZBVEMroxIPpuRGRERKzskDsGg4HNhotq8aAT3/Bb4B1sYlHk3JjYiIlIzfvjA3vTx9AvxDzQ0vm/W3OiqpAJTciIiIe+XnwtfjYP0bZrt6Gxg0FyLqWRuXVBhKbkRExH1O/A4Lh8GhzWb7mgeh+/Pg429lVFLBKLkRERH32PEpfPIw5KRDQBj0nwmNb7I6KqmAlNyIiMiVyTsDK/8Pvn/LbNe8CgbNgfDa1sYlFZaSGxERuXxpe+DjIXDkJ7N97aNw/f+Bt6+lYUnFpuRGREQuz/ZF8OmjkJsJgREw4E2I62l1VCJKbkREpJjyTsPyp2DzPLNduyPcOhvCalgalsg5Sm5ERKToUnaZw1DHfgFs0OUf0G0MeOvrRMoO/TSKiEjR/Phf+Hw05GVBcDQMfAsaXG91VCIFKLkREZGLy82CZU/Ctv+Y7bpd4NZ3IKSqtXGJXICSGxERubBjO81hqJRfweYFXZ+C+MfBy9vqyEQuSMmNiIgUZBiw9T+w7AnIPw2VYsyi4XpdrI5M5JKU3IiIiKucTLO25uePzHaD62HAW1Ap2tq4RIpIyY2IiJx35GdzGCptN9i84fpn4NrHwMvL6shEikzJjYiImMNQm+bA8jFgz4GQ6uYWCnU6Wh2ZSLEpuRERqejOpMNnj8IvS8x2bC9z08vgSGvjErlMSm5ERCqyw1vh46FwYh94+cANY6HjwxqGknJNyY2ISEVkGOYu3iueBXsuhNU2h6FqXWV1ZCJXTMmNiEhFc/oEfPIw/Pq52W58M/R7AwIrWxuXiJsouRERqUgObjKHodKTwMsXer4EHR4Am83qyETcRsmNiEhFYBiw/g346nlw5EPlujBoLtRoa3VkIm6n5EZExNNlH4elf4Ndy8120/5wy+sQEGZpWCIlRcmNiIgnS9oAC4dBxiHw9ocbJ0D7YRqGEo+m5EZExBM5HLD2VfjmX2DYIaIB3DYPqrW0OjKREqfkRkTE05xKgSX3w55vzHaL2+DmV8E/xNq4REqJkhsREU+ybzUsGgGnjoBPIPR5BdrcrWEoqVCU3IiIeAKHHRL+DateBsMBUY3MYaiYplZHJlLqlNyIiJR3mUdh8QjYl2C2W98FfSaBX7C1cYlYRMmNiEh5tudbWHwfZKWAbxDcNAVa/8XqqEQspeRGRKQ8sufDdxNg9WTAgCrNzGGo6DirIxOxnJIbEZHyJv2QWTSctM5stxsCN74MvoGWhiVSVhQpuRk9enSRX3DKlCmXHYyIiFxC4kpYfD+cPg5+laDva9BikNVRiZQpRUputm7d6tLevHkzdrudRo0aAbBr1y68vb1p166d+yMUERGw58E3L8La18x21ZbmMFRkA0vDEimLipTcfPvtt87/njJlCiEhIbz77rtUrlwZgBMnTjB06FC6dOlSMlGKiFRkJw+YWygc/N5sX3WfuZu3b4C1cYmUUTbDMIziPKFGjRqsWLGCZs2auRzfvn07PXv25PDhw24N0N0yMjIICwsjPT2d0NBQq8MREbm4X5eZm16eOQn+YdBvGjTtZ3VUIqWuON/fXpfz4kePHi1w/NixY2RmZhb35UREpDD5ubB8DPz3L2ZiU70tPLBKiY1IERR7ttSAAQMYOnQokydP5pprrgFgw4YNPPHEEwwcONDtAYqIVDjH98HCoXD4bL3jNQ9B9+fBx8/SsETKi2InN7NmzeLxxx/nr3/9K3l5eeaL+PgwfPhwJk2a5PYARUQqlB2fwCcPQ04GBIRD/5nQuI/VUYmUK8WuuTknKyuLPXv2YBgGDRs2JDi4fCzzrZobESmT8s7Aimfhh7fNds2rYdAcCK9lbVwiZURxvr8vexG/5ORkkpOTiY+PJzAwEMMwsGnXWRGR4kvbAx8PgSM/me1rH4Xr/w+8fS0NS6S8KnZyk5aWxu233863336LzWYjMTGR+vXrM2LECMLDw5k8eXJJxCki4pl+XgifPQq5pyAoEga8CbE9rI5KpFwr9mypxx57DF9fX5KSkggKCnIeHzx4MMuXL3drcCIiHivvNHz6CCwabiY2tTvByDVKbETcoNg9NytWrODLL7+kZs2aLsdjY2PZv3+/2wITEfFYKbvg43vh2A7ABvGPQ9enwFvb/Ym4Q7F/k7Kyslx6bM5JTU3F39/fLUGJiHisbR/C/0ZDXjYEV4GBb0GD66yOSsSjFHtYKj4+nvfee8/ZttlsOBwOJk2axHXX6RdURKRQuVmw9EFYOtJMbOrFm8NQSmxE3K7Yyc2kSZN488036d27N7m5uTz55JM0b96chIQEJk6cWOwAZsyYQb169QgICKBdu3asXr36otfn5OTwzDPPUKdOHfz9/WnQoAFz5swp9vuKiJSaozvgretg23yweUG3p+HupRASY3VkIh6p2MNSTZs25aeffmLmzJl4e3uTlZXFwIEDeeihh6hWrVqxXmvBggWMGjWKGTNmcO211zqTph07dlC7du1Cn3P77bdz9OhRZs+eTcOGDTl27Bj5+fnF/RgiIiXPMGDLe/DFk5B/BipVhVvfgXraZFikJF32In7u0KFDB9q2bcvMmTOdx5o0aUL//v2ZMGFCgeuXL1/OHXfcwd69e4mIiLis99QifiJSKnIy4fPH4OePzXaDG8xp3pWirY1LpJwq0Y0zly9fzpo1a5zt6dOn07p1a+68805OnDhR5NfJzc1l8+bN9OzZ0+V4z549WbduXaHP+fTTT2nfvj2vvPIKNWrUIC4ujscff5zTp09f8H1ycnLIyMhweYiIlKjkn+CtbmZiY/OGG8bCXQuV2IiUkmInN0888YQzQfj5558ZPXo0ffr0Ye/evYwePbrIr5OamordbicmxnXMOSYmhiNHjhT6nL1797JmzRq2b9/OkiVLmDp1KgsXLuShhx664PtMmDCBsLAw56NWLS1lLiIlxDDgh3fgne6QthtCa8DQZdBlNHgV+8+tiFymYtfc7Nu3j6ZNmwKwaNEi+vbty/jx49myZQt9+hR/c7c/b9lwsW0cHA4HNpuN+fPnExYWBsCUKVMYNGgQ06dPJzAwsMBzxowZ45J0ZWRkKMEREfc7k24uyrdjqdmOu9Hc9DLo8obQReTyFTu58fPzIzs7G4CvvvqKe+65B4CIiIhiDflERUXh7e1doJfm2LFjBXpzzqlWrRo1atRwJjZg1ugYhsHBgweJjY0t8Bx/f3+tvyMiJevQFlg4FE78Dl4+0H0cdHwItN+eiCWK3U/auXNnRo8ezYsvvsj333/PTTfdBMCuXbsKrFp8MX5+frRr146VK1e6HF+5ciWdOnUq9DnXXnsthw8f5tSpU85ju3btwsvLq1jvLSLiFoYBG2bC7J5mYhNWG4Z9CZ0eVmIjYqFiJzdvvPEGPj4+LFy4kJkzZ1KjRg0AvvjiC2688cZivdbo0aN55513mDNnDjt37uSxxx4jKSmJkSNHAuaQ0rmeIYA777yTyMhIhg4dyo4dO0hISOCJJ55g2LBhhQ5JiYiUmNMnYMFfYflT4MiDxjfDyASo2d7qyEQqvGIPS9WuXZvPP/+8wPFXX3212G8+ePBg0tLSeOGFF0hOTqZ58+YsW7aMOnXqAJCcnExSUpLz+kqVKrFy5Ur+/ve/0759eyIjI7n99tt56aWXiv3eIiKX7cAPsHAYpCeBtx/0/BdcfZ96a0TKiMta58Zut7NkyRJ27tyJzWajcePG9O/fHx+fsr/pm9a5EZHL5nDA+jfg63HgyIfK9eC2uVC9jdWRiXi84nx/Fzsb2b59O7fccgtHjx6lUaNGgFn3Eh0dzaeffkqLFi0uL2oRkbIsKw2W/g0SvzTbzQZA39cgIOzizxORUlfsmpsRI0bQvHlzDh48yJYtW9iyZQsHDhygZcuW3H///SURo4iItfavg1mdzcTG2x9ufhUGzVViI1JGFbvn5scff2TTpk1UrlzZeaxy5cr861//4qqrrnJrcCIilnI4YM0U+HY8GHaIbAi3zYOq6qEWKcuK3XPTqFEjjh49WuD4sWPHaNiwoVuCEhGx3KkUmH8rfPOimdi0HAz3r1JiI1IOFLvnZvz48TzyyCM8//zzXHPNNQBs2LCBF154gYkTJ7os5KeCXREpl/YlwKIRcOoo+ARCn0nQ5q+aDSVSThR7tpTXH/ZHObdNwrmX+GPbZrNht9vdFafbaLaUiFyQww4Jk2DVRDAcEN3YHIaq0sTqyEQqvBKdLfXtt99edmAiImVW5hGzt+b31Wa79V+hzyvgF2xtXCJSbMVObrp27VoScYiIWGf317D4fshOBd9guHkKtLrD6qhE5DJd9qp72dnZJCUlkZub63K8ZcuWVxyUiEipsOfDd+Nh9RTAgJjm5hTv6DirIxORK1Ds5CYlJYWhQ4fyxRdfFHq+LNbZiIgUkH7IHIZKWme22w2FGyeAr/apEynvij0VfNSoUZw4cYINGzYQGBjI8uXLeffdd4mNjeXTTz8tiRhFRNxr1wpzUb6kdeAXAoPmQN+pSmxEPESxe26++eYbPvnkE6666iq8vLyoU6cOPXr0IDQ0lAkTJnDTTTeVRJwiIlfOngdfvwDrXjfb1VqZw1CRDayNS0Tcqtg9N1lZWVSpUgWAiIgIUlJSAGjRogVbtmxxb3QiIu5yMgnm9j6f2Fz9AAxfqcRGxAMVu+emUaNG/Pbbb9StW5fWrVvz5ptvUrduXWbNmkW1atVKIkYRkSvz6//MTS/PpIN/GPR7A5reYnVUIlJCip3cjBo1iuTkZADGjh1Lr169mD9/Pn5+fsybN8/d8YmIXL78XFj5HGycabZrtDPrayrXtTQsESlZxV6h+M+ys7P59ddfqV27NlFRUe6Kq8RohWKRCuL4Plg4FA5vNdsdH4YbxoKPn7VxichlKdEViv8sKCiItm3bXunLiIi4zy9L4dO/Q04GBITDgFnQqLfVUYlIKSlScjN69Ogiv+CUKVMuOxgRkSuSdwa+fBo2zTbbtTrArbMhvJa1cYlIqSpScrN169YivZhNO+aKiFVSd8PHQ+Doz2a782Nw3TPg7WtpWCJS+oqU3GizTBEp0376GD4fBbmnICgSBrwFsd2tjkpELHLZNTe7d+9mz549xMfHExgYiGEY6rkRkdKVmw3L/wlb3jPbdTrDre9AqJalEKnIir2IX1paGjfccANxcXH06dPHOS18xIgR/OMf/3B7gCIihUr5Dd654WxiY4P4J+GeT5TYiEjxk5vHHnsMX19fkpKSCAoKch4fPHgwy5cvd2twIiKF2vYBvNUNju2A4Cpw9xK4/hnwvuIJoCLiAYr9l2DFihV8+eWX1KxZ0+V4bGws+/fvd1tgIiIF5JyCZY/Djx+a7XpdYeDbEBJjbVwiUqYUO7nJyspy6bE5JzU1FX9/f7cEJSJSwNFfzNlQqbvA5gXdnoYuo8HL2+rIRKSMKfawVHx8PO+9956zbbPZcDgcTJo0ieuuu86twYmIYBiw+V14+3ozsQmpBvd+Bl2fUGIjIoUqds/NpEmT6NatG5s2bSI3N5cnn3ySX375hePHj7N27dqSiFFEKqqcTPhsFGxfaLYbdocBb0Jw2d/qRUSsU+yem6ZNm/LTTz9x9dVX06NHD7Kyshg4cCBbt26lQYMGJRGjiFREyT/Cm/FmYmPzhu7Pw50fK7ERkUu64o0zyxttnClSxhkG/PCOuY2CPRdCa5o7edfuYHVkImKhUt04U0TEbU6fhM8egR2fmO243tB/BgRFWBqWiJQvSm5EpGw4tBk+Hgon94OXL/QYB9c8CFr5XESKScmNiFjLMGDDTFj5HDjyILw2DJoHNdtZHZmIlFNKbkTEOtnH4ZOH4LdlZrtJX7jlDQgMtzQsESnfij1bCiA/P5+vvvqKN998k8zMTAAOHz7MqVOn3BqciHiwA9+bs6F+WwbeftB7Etz+vhIbEblixe652b9/PzfeeCNJSUnk5OTQo0cPQkJCeOWVVzhz5gyzZs0qiThFxFM4HLDudfj6BTDsULke3DYPqre2OjIR8RDF7rl59NFHad++PSdOnCAwMNB5fMCAAXz99dduDU5EPExWGnxwO3w11kxsmg2EBxKU2IiIWxW752bNmjWsXbsWPz8/l+N16tTh0KFDbgtMRDzM/nWwcDhkHgafALjxZWg3RLOhRMTtip3cOBwO7HZ7geMHDx4kJCTELUGJiAdxOGDNZPh2PBgOiIw1h6GqNrc6MhHxUMUelurRowdTp051tm02G6dOnWLs2LH06dPHnbGJSHl36hj8ZyB885KZ2LS8A+7/TomNiJSoYm+/cPjwYa677jq8vb1JTEykffv2JCYmEhUVRUJCAlWqVCmpWN1C2y+IlJK9q2DxfXDqKPgEwk3/htZ3aRhKRC5LiW6/UL16dbZt28aHH37Ili1bcDgcDB8+nLvuusulwFhEKiiHHVZNhFWvAAZENzGHoao0tjoyEakgtHGmiLhPRrLZW/P7arPd5m7o/Qr4BVkbl4iUeyW+ceahQ4dYu3Ytx44dw+FwuJx75JFHLuclRaS82/01LL4fslPBNxj6ToWWt1sdlYhUQMVObubOncvIkSPx8/MjMjIS2x/Gz202m5IbkYrGng/f/gvWTDHbMS3MYaiohpaGJSIVV7GHpWrVqsXIkSMZM2YMXl6XtXuDpTQsJeJG6QfNtWsObDDb7YdBr/Hgq/o7EXGvEh2Wys7O5o477iiXiY2IuNGuL2HJA3D6BPiFwC2vQ/OBVkclIlL8dW6GDx/Oxx9/XBKxiEh5YM+DL58xt1E4fQKqtYaRCUpsRKTMKPawlN1u5+abb+b06dO0aNECX19fl/NTpkxxa4DupmEpkStwYj8sHAaHNpntDiOhxwvg429tXCLi8Up0WGr8+PF8+eWXNGrUCKBAQbGIeKidn8MnD8KZdAgIg37ToUlfq6MSESmg2MnNlClTmDNnDkOGDCmBcESkzMnPgZXPwcZZZrtGOxg0FyrXsTYuEZELKHZy4+/vz7XXXlsSsYhIWXN8L3w8FJK3me2OD8MNY8HHz9KwREQuptgFxY8++ijTpk0riVhEpCzZvhhmxZuJTWBl+MsC6PUvJTYiUuYVu+fm+++/55tvvuHzzz+nWbNmBQqKFy9e7LbgRMQCeWfgyzGwaY7ZrnUNDJoNYTWtjUtEpIiKndyEh4czcKCmfIp4pNTd8PEQOPqz2e48Gq57Brwva6cWERFLXNb2C+40Y8YMJk2aRHJyMs2aNWPq1Kl06dLlks9bu3YtXbt2pXnz5mzbts2tMYlUSD99BJ+NgrwsCIqCgW9Cw+5WRyUiUmyWLjO8YMECRo0axTPPPMPWrVvp0qULvXv3Jikp6aLPS09P55577uGGG24opUhFPFhuNnzysLmbd14W1O0CI9cosRGRcqtIi/i1bduWr7/+msqVK9OmTZuLrmezZcuWIr95hw4daNu2LTNnznQea9KkCf3792fChAkXfN4dd9xBbGws3t7eLF269KI9Nzk5OeTk5DjbGRkZ1KpVS4v4iQAc+9UchkrZCdig65PQ9Z/g5W11ZCIiLty+iF+/fv3w9zdXIO3fv/8VBwiQm5vL5s2beeqpp1yO9+zZk3Xr1l3weXPnzmXPnj385z//4aWXXrrk+0yYMIFx48ZdcbwiHmfrfPjfPyD/NFSKgYFvQ/2uVkclInLFipTcjB07lmHDhvHaa68xduxYt7xxamoqdrudmJgYl+MxMTEcOXKk0OckJiby1FNPsXr1anx8ilYuNGbMGEaPHu1sn+u5Eamwck7Bssfhxw/Ndv1uZmJTqYqlYYmIuEuRa27effddTp8+7fYA/jzEZRhGocNedrudO++8k3HjxhEXF1fk1/f39yc0NNTlIVJhHdkOb19nJjY2L7j+WfjrEiU2IuJRijxbqpj7a15SVFQU3t7eBXppjh07VqA3ByAzM5NNmzaxdetWHn74YQAcDgeGYeDj48OKFSu4/vrr3RqjiMcwDNg8D5Y/BflnIKQa3Dob6mq1cRFxj8wzeazbk8aqXSl422y82L+5ZbEUayq4OzfG9PPzo127dqxcuZIBAwY4j69cuZJ+/foVuD40NJSff/7Z5diMGTP45ptvWLhwIfXq1XNbbCIe5UwGfD4Kti8y2w17wIBZEBxlaVgiUr45HAa/HM5g1a5jJOxKZUvSCfIdZkdIsJ83/3dzU/x8rJmUXazkJi4u7pIJzvHjx4v8eqNHj+buu++mffv2dOzYkbfeeoukpCRGjhwJmPUyhw4d4r333sPLy4vmzV2zwCpVqhAQEFDguIiclfyjORvq+F6weUP3sdDx7+Bl6SoQIlJOpWTmsDoxhYRdKaxOTCUtK9flfL2oYOJjo4iPi8aN/SHFVqzkZty4cYSFhbntzQcPHkxaWhovvPACycnJNG/enGXLllGnjrnbcHJy8iXXvBGRQhgG/PAOfPk02HMhrBYMmgO1rrY6MhEpR3LzHWzef4KEswnNL4czXM4H+3nTqaGZzHSNjaZ2ZJBFkboq0jo3AF5eXhw5coQqVcp34WFx5smLlEunT8Knf4edn5rtRn2g33QIirA0LBEpH/anZZGwK4VVu1JZvyeVrFy7y/nmNUKJj40mPi6atrUrl9rQk9vXuQH31tuISAk5uBkWDoGTSeDlCz1egGv+hqX9wyJSpmXl5LN+T5qzd+b3tGyX85HBfsTHRRMfF0XnhtFEh/hbFGnRWTZbSkTcyDBgwwxYORYceRBeB26bCzXaWR2ZiJQxhmGwMzmThMQUVv2Wwqb9x8mzn/+O9/Gy0bZOZbrGRdM1Lpqm1ULx8ipf/0AqcnLjcDhKMg4RuVzZx2Hpg7DrC7Pd5Ba4ZRoEhlsaloiUHcezcs8WAqeSkJhCSmaOy/laEYHEx5rJTMcGkYQE+FoUqXsUe1dwESlDkjbCwmGQcRC8/aDXeLhqhIahRCq4fLuDrQdOnq2dSeHnQ+n8cQAm0Nebjg0iiY+NomujKtSNDPKo8hMlNyLlkcMB616Hr18Aww4R9eG2eVCtldWRiYhFDp7INntmdqWwdncqmTn5LucbVw2ha5xZCNy+bmX8fTx3g1wlNyLlTVYqLBkJu1ea7eaDoO9U8A+xNCwRKV2nc+1s3GeuCJywK4U9KVku58ODfOkSG+1cdyYmNMCiSEufkhuR8uT3tbBoOGQmg08A9J4Ibe/VMJRIBWAYBonHTjmHmjbuO05u/vl6WC8btKld2dk706JGGN7lrBDYXZTciJQHDjusngLfjQfDAVFx5jBUTDOrIxOREpSencea3ams2nWM1YmpJKefcTlfPSzAXEAvLppODaMICyzfhcDuouRGpKw7dQwWjYB9q8x2q79An3+DfyVr4xIRt7M7DH48eL4Q+McDJ3H8oRDY38eLDvXNQuBujaJpEF3JowqB3UXJjUhZtvc7WHQfZB0D3yC4aTK0vtPqqETEjY6knzGTmcQU1iSmkn46z+V8bJVKZxfRi6ZDvQgCfD23ENhdlNyIlEUOO6yaCKteAQyo0hQGzYUqja2OTESu0Jk8Oz/8fpyEXea6M78dzXQ5HxrgQ+fYKOcWB9XDAy2KtPxSciNS1mQkm8NQ+9eY7bb3wI0Twa9sbEgnIsVjGAZ7U7NY9VsKCYkpbNibxpm884XANhu0rBl+dkXgKFrVDMfHu3T2a/JUSm5EypLEr2DJ/ZCdBn6V4Oap0PI2q6MSkWLKOJPHut3np2kfOnna5XyVEH9nIXDnhlFUDvazKFLPpORGpCyw58G3/4I1r5rtmBbmbKiohpaGJSJF43AYbD+c7iwE3pJ0EvsfKoH9vL24ql5lc4uDRtE0iglRIXAJUnIjYrX0g+YWCgc2mu2rRkDPf4FvxVlwS6Q8OpZ5htVn92panZjK8axcl/P1o4KdvTMd6kcQ5Kev3NKiOy1ipd+Ww9KRcPoE+IfCLa9DswFWRyUihcjNd7Bp/3HnFgc7kjNczlfy96FTg0hnQlMrQnVyVlFyI2KF/Fz4ehysf8NsV2sNt80194gSkTJjf1qWs25m/Z40snLtLueb1wg1VwSOjaZtncr4qhC4TFByI1LaTvxuDkMd2my2O/wNeowDH39LwxIRyMrJZ/2es4XAiSnsT8t2OR9VyY8usWcLgWOjiKqk39uySMmNSGna8Sl88jDkpENAGPSbAU1utjoqkQrLMAx2JGeQsMvc4mDz/hPk2c8XAvt42WhXpzJdG5m9M02rheJVQfdrKk+U3IiUhvwcWPEsfP+W2a55FQyaA+G1rY1LpAJKO5Vzdr8mcxG91FM5LudrRQSeXXOmCh0bRFLJX1+V5Y3+HxMpaWl7YOFQSP7RbHd6BG54Dry1wZ1IacizO9iaZO7XlJCYws+H0jH+sF9ToK+3sxA4Pi6aupFBmqZdzim5ESlJ2xfDp49AbiYERsCAWRDXy+qoRDzegePZJCSahcDrdqeRmZPvcr5JtVDi46LoGhtNu7qV8ffRfk2eRMmNSEnIOw3Lx8DmuWa7dke4dTaE1bA2LhEPdTrXzoZ9ac4tDvamZLmcrxzkS5ezezXFx0ZRJVTrSHkyJTci7paaCB8PgaPbARt0GQ3dngZv/bqJuIthGOw6esq5IvD3vx8nN//8fk3eXjba1DL3a4qPi6Z5jTC8VQhcYeivrYg7/bgAPn8M8rIgKAoGvgUNb7A6KhGPcDI7lzW7U527aR/JOONyvkZ4IPFx5m7anRpGERaouraKSsmNiDvkZsOyJ2Dbf8x23S5w6zsQUtXauETKMbvDYNuB84XAPx44yR+2a8Lfx4tr6p9bETiKBtGVVAgsgJIbkSt3bKc5DJXyK2CDbk9B/BPgpQJFkeJKTj/t7JlZszuV9NN5LufjYioRf7Z25up6EQT46vdMClJyI3K5DAO2zYf/PQ75p6FSjNlbUy/e6shEyo0zeXZ++P24sxB419FTLudDA3zOFgJH0SU2murhgRZFKuWJkhuRy5FzCv43Gn5aYLbrX2fW11SqYm1cImWcYRjsSclyFgJv3JfGmbzzhcA2G7Sqeb4QuFXNMHy0X5MUk5IbkeI68rM5DJW2G2xecN0z0Hk0eOkPsEhhMs7ksW53KqvO7qZ96ORpl/Mxof7Ex0bTtVE0nRtGER7kZ1Gk4imU3IgUlWGY69Z88RTYcyCkOgyaDXU6WR2ZSJnicBj8fCjdWQi8Jekk9j9UAvt5e3F1vQhzEb24KsTFqBBY3EvJjUhRnMmAzx6FXxab7die0H8WBEdaG5dIGXEs4wwJiWbPzJrdqRzPynU5Xz862OydiYumQ/0Igvz09SMlRz9dIpdyeJs5DHViH3j5wA1joePDGoaSCi0338Gm/cedm0/uTM5wOV/J34drG57dryk2mloRQRZFKhWRkhuRCzEMcxfvFc+CPRfCasGguVDrKqsjE7HE76lZJCSmsOq3FNbvTSM71+5yvkWNMGchcJva4fiqEFgsouRGpDCnT8AnD8Ovn5vtRjdBvzcgKMLauERK0amcfNbvSXPObEo6nu1yPqqSP/GxUc5C4MhK/hZFKuJKyY3Inx3cBAuHwskk8PKFni9Ch5HmHFURD+ZwGOxIznDupr15/wny7OcLgX29bbSrU/nsisDRNKkaipf2a5IySMmNyDmGAeunw1djwZEPleuaw1A12lodmUiJSTuVw+qzhcAJiamknspxOV87Isg51NSxQSSV/PW1IWWffkpFALKPw9K/wa7lZrtpP7hlGgSEWRuXiJvl2R1s2X/ibO9MKj8fSnc5H+TnTacG5wuB60YFWxSpyOVTciOStAEWDoeMg+DtDzeOh/bDNQwlHuPA8WxnIfC6PWmcysl3Od+kWujZ3pko2tWpjL+P9muS8k3JjVRcDgesnQrfvASGHSIawG3zoFpLqyMTuSLZufls3HtumnYKe1OzXM5HBPvRuWEUXeOi6RIXRZWQAIsiFSkZSm6kYjqVAksegD1fm+0Wt8HNr4J/iLVxiVwGwzD47Wimczft7/cdJ9d+fr8mby8bbWuHO7c4aF49TIXA4tGU3EjF8/sacxjq1BHwCYA+k6DN3RqGknLlRFYua3anOrc4OJrhWghcIzzw7KymKDo1jCI0wNeiSEVKn5IbqTgcdlg9Gb6bAIYDohqZw1AxTa2OTOSS8u0OfjyY7hxq+vHgSYzzs7QJ8PXimvqRxMeaM5saRAdrvyapsJTcSMWQeRQW3wf7Vpnt1neZPTZ+mgkiZVdy+mnnAnprElPJOONaCBwXU8k5TfuquhEE+KoQWASU3EhFsOdbWHw/ZB0D3yC4aQq0/ovVUYkUcCbPzvf7zhcCJx475XI+LNCXzrFRdI01C4GrhQVaFKlI2abkRjyXPR9WvQwJ/wYMqNLUHIaKbmR1ZCKAWQi8J+UUq3alsmpXChv3ppGTf74Q2MsGrWqdLwRuVTMcbxUCi1ySkhvxTBmHYdEI2L/WbLe9F3pPBF/9S1eslX46j3W7U53rzhxOP+NyvmpoAPFxUcTHmfs1hQf5WRSpSPml5EY8T+JXsOR+yE4Dv0rQ9zVoMcjqqKSCsjsMth86Xwi89cBJ7I7zlcB+Pl50qBfhLASOi6mkQmCRK6TkRjyHPc9ckG/tVLNdtQXc9i5ENrA0LKl4jmWcMZOZxFTWJKZwIjvP5XyD6GBze4O4aK6pF0mgnwqBRdxJyY14hpMHYNFwOLDRbF91H/R8CXy18qqUvJx8O5t/P8GqszObfj2S6XI+xN+HTg0j6RpXhfi4KGpWDrIoUpGKQcmNlH+/LjM3vTxzEvxDzQ0vm/W3OirxYIZh8HtatnOa9vo9aZzOszvP22zQokaYsxC4da1wfL29LIxYpGJRciPlV34ufPU8bJhutqu3gUFzIaKepWGJZzqVk+8sBE7YlUrS8WyX81GV/ImPM/dr6twwishK/hZFKiJKbqR8OvE7fDwUDm8x29c8CN3HgY9mloh7OBwGO5IznIXAm/efIP8PhcC+3jba14k4u8VBNE2qhagQWKSMUHIj5c+OT+GThyEnHQLCof9MaNzH6qjEA6SeymH12Z6Z1YkppJ7KdTlfJzLIXBE4NpqODSIJ9tefUJGyyPLfzBkzZjBp0iSSk5Np1qwZU6dOpUuXLoVeu3jxYmbOnMm2bdvIycmhWbNmPP/88/Tq1auUoxZL5J2Blf8H379ltmteDYNmQ3hta+OScivP7mDL/hNnZzalsP1Qhsv5ID9vOjWIdG5xUCdS23WIlAeWJjcLFixg1KhRzJgxg2uvvZY333yT3r17s2PHDmrXLviFlZCQQI8ePRg/fjzh4eHMnTuXvn37snHjRtq0aWPBJ5BSk7YHPh4CR34y29c+Ctf/H3hrp2MpngPHs52zmtbvSeNUjut+TU2rhTqHmtrVqYyfjwqBRcobm2H8cV/Z0tWhQwfatm3LzJkznceaNGlC//79mTBhQpFeo1mzZgwePJjnnnuu0PM5OTnk5OQ42xkZGdSqVYv09HRCQ0Ov7ANI6fh5IXw2CnIzITACBrwJcT2tjkrKiezcfDbsTSPh7BYH+1KzXM5HBPvRJfZsIXBsFFVCtHyASFmUkZFBWFhYkb6/Leu5yc3NZfPmzTz11FMux3v27Mm6deuK9BoOh4PMzEwiIiIueM2ECRMYN27cFcUqFsk7Dcufgs3zzHbtTnDrOxBWw9KwpGwzDINfj2SScHao6Yd9J8i1n9+vydvLRrvalc/ObKpCs+qheGm/JhGPYllyk5qait1uJyYmxuV4TEwMR44cKdJrTJ48maysLG6//fYLXjNmzBhGjx7tbJ/ruZEyLmWXOQx17BfABvGPQ9enwNvyMjEpg05k5bJ6dyoJu1JYnZjC0Ywcl/M1wgPp2sgsBO7UMJLQAA1ningyy78p/jx10jCMIk2n/PDDD3n++ef55JNPqFKlygWv8/f3x99f602UKz/+Fz4fDXlZEBwNA9+CBtdbHZWUIfl2Bz8ePMmq31JYlZjKTwdP8scB9gBfLzrWj3RucVA/KljTtEUqEMuSm6ioKLy9vQv00hw7dqxAb86fLViwgOHDh/Pxxx/TvXv3kgxTSlNuFix7ArbNN9t1u5jDUCFVrY1LyoTDJ087VwReuzuVjDOuhcCNYkKcvTPt61YmwFf7NYlUVJYlN35+frRr146VK1cyYMAA5/GVK1fSr1+/Cz7vww8/ZNiwYXz44YfcdNNNpRGqlIZjO81hqJRfweZlDkHFPw5e+oKqqM7k2dm477gzodl97JTL+bBAXzqfLQSOj42mapgKgUXEZOmw1OjRo7n77rtp3749HTt25K233iIpKYmRI0cCZr3MoUOHeO+99wAzsbnnnnt47bXXuOaaa5y9PoGBgYSFhVn2OeQKGAZsfR+WPQn5p6FSVbO3pl7hax2J5zIMg93HTjl30964N42c/POFwF42aF0r3DlNu2XNcLxVCCwihbA0uRk8eDBpaWm88MILJCcn07x5c5YtW0adOnUASE5OJikpyXn9m2++SX5+Pg899BAPPfSQ8/i9997LvHnzSjt8uVI5mWZtzc8fme0G18OAt6BStLVxSalJP53H2rOFwAm7UjicfsblfLWwAOJjzbqZzg2jCAtSIbCIXJql69xYoTjz5KUEHfnZHIZK2w02b7j+Wbh2FHhpwTRPZncY/HwonVW/mdO0tx04if0P+zX5+XjRoV6Ec0Xg2CqVVAgsIkA5WedGKijDgE1zYPkYsOdAaA24dTbU6Wh1ZFJCjmaccdbNrNmdysnsPJfzDaKD6RpXhfi4KDrUiyTQT3VWInJllNxI6TmTDp89Cr8sMduxvWDALAi68CKMUv7k5NvZ9PsJZ0Lz65FMl/Mh/j5c2zDKnNkUF02N8ECLIhURT6XkRkrHoS2wcCic+B28fKD783DNQxqG8gCGYbAvNevsisCprN+Txuk8u/O8zQYta4Q5C4Fb1wrHx1v/v4tIyVFyIyXLMGDjm7DiWXDkQVhtGDQHal1ldWRyBTLP5LFuT5pzi4MDx0+7nI8O8T9bCBxFl9hoIoL9LIpURCoiJTdSck6fgE8ehl8/N9uNb4Z+b0BgZWvjkmJzOAx2JGc4d9Pesv8E+X8oBPb1tnFV3QhzReDYaJpUC1EhsIhYRsmNlIyDm+DjoZCeBN5+0PMluPp+c4xCyoXUUzmsTkxh1W8prE5MJS0r1+V83cgg56yma+pHEuyvPyciUjbor5G4l8MBG6bDV8+DIx8q14Xb5kH1NhYHJpeSm+9gS9L5QuBfDme4nA/286ZjA7MQuGtsNLUjgyyKVETk4pTciPtkH4clIyHxS7PdtD/c8joEaPXosiopLZtVieYCeut2p5KVa3c536x6qLMQuG3tyvj5qBBYRMo+JTfiHvvXw6LhkHEIvP3hxgnQfpiGocqYrJx8NuxNc85s2pea5XI+MtiPLrFRxMdF0yU2mugQf4siFRG5fEpu5Mo4HLD2VfjmX2DYIbKhOQxVtYXVkQnmNO1fj2Sa+zXtSmHT7yfItZ/fr8nHy0bbOpWdm082qx6Kl/ZrEpFyTsmNXL5TKbDkftjzjdlucTvcPAX8Q6yNq4I7kZXL6t2pZwuBUziWmeNyvmblQGchcKcGkYQEaL8mEfEsSm7k8uxbDYtGwKkj4BMIfSZBm79qGMoC+XYH2w6cdBYC/3QonT/uGBfo68019c/v11QvKljTtEXEoym5keJx2CFhEqyaCIYDohrB7e9ClSZWR1ahHDp52rmT9prdqWSeyXc537hqiLMQuH3dyvj7aL8mEak4lNxI0WUegcX3wb4Es936r9DnFfALtjauCuBMnv1sIXAqq3YdY0+KayFweJAvnRtGORfRqxoWYFGkIiLWU3IjRbPnG1h8P2SlgG+wWVvT6g6ro/JYhmGw+9gp54rAG/cdJzf/fCGwlw3a1K5MfGw0XRtF06JGGN4qBBYRAZTcyKXY8+G7CbB6MmBAlWbmbKjoOKsj8zjp2Xms3WMWAickppCcfsblfLWwAGfdzLUNoggLUiGwiEhhlNzIhaUfMouGk9aZ7XZDzfVrfAOtjctD2B0GPx086Rxq2nbgJH/Yrgl/Hy+urmcWAneNi6ZhlUoqBBYRKQIlN1K4XStgyQNw+jj4hUDfqdBikNVRlXtHM844h5rW7k7lZHaey/mGVSo5e2c61IsgwFeFwCIixaXkRlzZ8+DrF2Dd62a7aktzGCqygaVhlVc5+XZ+2HeChLMbUP52NNPlfEiAD50bRtE1LpoucdHUCFevmIjIlVJyI+edTIKFw+DgD2b76vuhx4vgq5k3RWUYBvtSs5wrAm/Ye5zTeef3a7LZoGXNcLqe3eKgda1wfLy1X5OIiDspuRHTr/+DpQ/CmZPgHwb9pkHTflZHVS5knslj7e40Es5uQHnwxGmX81VC/M0p2nHRdG4YRUSwn0WRiohUDEpuKrr8XPhqLGyYYbart4Xb5kLlupaGVZY5HAa/HM5g1a5jJOxKZUvSCfL/UAns5+3FVfXMadrxcdE0rhqiQmARkVKk5KYiO/E7fDwUDm8x2x0fhhvGgo96Fv4sJTOH1Wd7ZlYnppKWletyvl5U8NlC4CiuqR9JkJ9+tURErKK/wBXVzs9g6UOQkw4B4TBgFjTqbXVUZUZuvoPN+084h5p+OZzhcj7Yz5tOZ1cE7hobTe3IIIsiFRGRP1NyU9H8eRiq5lUwaC6E17I2rjJgf1rW2c0nU1m/J5WsXLvL+eY1Qp1DTW1rV8bPR4XAIiJlkZKbiuTEflg4FA5tNtsdH4buz4N3xVzpNisnn/V7zhcC/56W7XI+qpIfXWLNoabODaOJDvG3KFIRESkOJTcVxa//g6V/gzNnh6H6z4TGfayOqlQZhsHO5EznmjOb9h8nz36+ENjHy0bbOpWdKwI3rRaKl/ZrEhEpd5TceLr8XPjqedgw3WzXaG/OhgqvbWlYpeV4Vu7ZQuBUEhJTSMnMcTlfKyLQ3HwyLpqODSIJCaiYvVgiIp5EyY0nO5lkzoY6tMlsV4DZUPl2B1sPnDxbO5PCz4fSMf6wX1OgrzcdG0Q6tzioGxmkadoiIh5GyY2n+u0LWDLSXJQvIOzsMNRNVkdVIg6eyDZ7Zs7u15SZk+9yvnHVEOdQU7u6lfH30X5NIiKeTMmNm5zJs/PZj4ex2WwMalfTukDseeYw1Po3zHaNduZsqMp1rIvJzU7n2tm4L825xcGelCyX85WDfOkcG0382S0OYkK1fYSISEWi5MZNMs7k8cTCn/CyYV1yc/KAORvq3N5Q1zwI3ceV+2EowzBIPHbKOdS0cd9xcvMdzvNeNmhbu7Jzi4MWNcLwViGwiEiFpeTGzYxLX1IyflsOSx44vzdU/xnQ5Garorli6dl5rNmdyqpdx1idmEpy+hmX89XDAujaKJr42Gg6NYwiLFCFwCIiYlJy4yY2LOopsOfB1+Ng3TSzXb0N3Dav3O0NZXcY/HjwfCHwjwdO8oftmvD38aJD/ciztTNRNIiupEJgEREplJKb8iz9oDkb6uD3ZrvD36DHOPApH4vNHUk/YyYziSmsSUwl/XSey/nYKpXM7Q3iorm6XgQBvioEFhGRS1Ny42ZGaY1L7VoBS+6H0yfMYah+b0DTW0rpzS/PmTw7P/x+nIRd5rozvx3NdDkfGuBD59go5xYH1cMDLYpURETKMyU3blJqIyT2PPjmJVg71WxXa20OQ0XUK6UAis4wDPamZrHqtxQSElPYsDeNM3nnC4FtNmhVM/xs70wUrWqG4+Ot/ZpEROTKKLkpT9IPwsLhcGCD2b76Aej5Ypkahso4k8e63eenaR86edrlfJUQf+cCep0bRlE5uHzP5BIRkbJHyY2blHjHTeJKWHw/nD4O/qFwyzRo1r+k3/WSHA6D7YfTnYXAW5JOYv9DJbCftxdX1avsTGgaxYSoEFhEREqUkpsSYBiG+77A7fnw7Uuw5lWzXa3V2WGo+u55/ctwLPMMq8/u1bQ6MZXjWbku5+tHBTsLgTvUjyDITz9mIiJSevSt4yYl0huRcRgWDoOk9Wb76vuh50ulPgyVm+9g0/7jzi0OdiRnuJyv5O9DpwaRzoSmVkRQqcYnIiLyR0puSoBhuKHAOPErczZUdhr4hUC/adBsgFviK4r9aVnOupn1e9LIyrW7nG9RI4z4OHNmU9s6lfFVIbCIiJQRSm7cxG39NvZ8+G48rJ5stqu2NIehIhu46x0KlZWTz/o9ZwuBE1PYn5btcj6qkp9zinbn2CiiKpWdImYREZE/UnJTAi57qZuMw+ZsqKR1ZvuqEdDzX+Dr/o0fHQ6DHckZrE40h5o27T9Onv185D5eNtrVqezc4qBptVC8tF+TiIiUA0pu3OSKh6F2f23OhspONYehbnkNmt/qltjOST2Vw+pEcwG91YmppJ7KcTlfOyKI+LgousZVoWODSCr568dDRETKH317lQDDMCjyQJU9H76bcHYYyoCqLeC2d90yDJWb72Dz/hMkJJq1M78cdi0EDvLzpmP984XAdaOCr/g9RURErKbkxk0ua+PMjGRYNAL2rzHb7YdBrwlXNAz1e2qWM5kprBC4abVQ4uOiiY+Lon2dCPx8VAgsIiKeRcmNVfZ8A4vuOzsMVQn6vgYtBhX7ZU7l5LNud+rZhCaVpOMFC4G7xEbTJTaKLrHRRIeoEFhERDybkhs3utq2kxE+yzCO1IQarQq/yGGH716GhEmAATEtzNlQUQ2L9B4Oh8EvhzNISDy7IvD+E+T/YUVgX2+zEDg+ToXAIiJSMSm5cRcb/NXnK3p6b8aYPxCGLYeoWNdrMo+Yw1C/rzbb7YbAjS+D78V3v/7jisBrElNJ+9OKwHUjg5zJzDUqBBYRkQpO34JuFIo5JGTLToX3+psJTngt8+Te78zEJivFHIa6eSq0vK3Q18nJt7Pp9xMk7EohITGVnYWsCNzx3IrAsdHUjtSKwCIiIucouXETmw2CbeYO2Ia3P7aMg/B+fxjyP9g0F1ZNxByGan52GOp8r45hGOxNzTKTmV0pbNh7nNN5dpfXbl5dKwKLiIgUhZIbN6qEmdzk956E7+pJkLYbXm8DeWeLfNveC70ngm8gGWfyWLc7lVVn92s6dPK0y2tFh/jTJTaKrnHRdG4YRaRWBBYRESkSy5ObGTNmMGnSJJKTk2nWrBlTp06lS5cuF7x+1apVjB49ml9++YXq1avz5JNPMnLkyFKMuHA2oBJnADCiG8PdS2HujeYwlG8wjptf5aeIXiQkHCRhVwpbD5zE/odCYD9vL66qV5n42Gi6xEbTpFpIyWzGKSIi4uEsTW4WLFjAqFGjmDFjBtdeey1vvvkmvXv3ZseOHdSuXbvA9fv27aNPnz7cd999/Oc//2Ht2rU8+OCDREdHc+ut7l3N93JUOjcs5VsJohqSOmgxqavn8LGjG4s+CeJk9lqX6+tHBxMfay6g16F+BEF+lueaIiIi5Z7NMJfTtUSHDh1o27YtM2fOdB5r0qQJ/fv3Z8KECQWu/+c//8mnn37Kzp07ncdGjhzJjz/+yPr164v0nhkZGYSFhZGenk5oaOiVf4izTp3Jw39CDL42O1OaL2XFQW9+PZLpck1IgA/XNohyLqJXs7IKgUVERIqiON/flnUV5ObmsnnzZp566imX4z179mTdunWFPmf9+vX07NnT5VivXr2YPXs2eXl5+Pr6FnhOTk4OOTnn91DKyMgocI072NJ24Wszi4DnbkolkyBsNmhZM5yusWZC07pWOD4qBBYRESlRliU3qamp2O12YmJiXI7HxMRw5MiRQp9z5MiRQq/Pz88nNTWVatWqFXjOhAkTGDdunPsCv4CggECO+tYkOT+Unm1i6dq4Cp0bRhER7Ffi7y0iIiLnWV7k8eeiWcMwLlpIW9j1hR0/Z8yYMYwePdrZzsjIoFatWpcb7oXjiqxPzGOriXHk07pSFbe/voiIiBSNZclNVFQU3t7eBXppjh07VqB35pyqVasWer2Pjw+RkZGFPsff3x9//1KaRh0UUTrvIyIiIhdkWQGIn58f7dq1Y+XKlS7HV65cSadOnQp9TseOHQtcv2LFCtq3b19ovY2IiIhUPJZWt44ePZp33nmHOXPmsHPnTh577DGSkpKc69aMGTOGe+65x3n9yJEj2b9/P6NHj2bnzp3MmTOH2bNn8/jjj1v1EURERKSMsbTmZvDgwaSlpfHCCy+QnJxM8+bNWbZsGXXq1AEgOTmZpKQk5/X16tVj2bJlPPbYY0yfPp3q1avz+uuvl4k1bkRERKRssHSdGyuU1Do3IiIiUnKK8/2tRVdERETEoyi5EREREY+i5EZEREQ8ipIbERER8ShKbkRERMSjKLkRERERj6LkRkRERDyKkhsRERHxKEpuRERExKNYuv2CFc4tyJyRkWFxJCIiIlJU5763i7KxQoVLbjIzMwGoVauWxZGIiIhIcWVmZhIWFnbRayrc3lIOh4PDhw8TEhKCzWZz62tnZGRQq1YtDhw4oH2rSpDuc+nQfS4dus+lR/e6dJTUfTYMg8zMTKpXr46X18Wraipcz42Xlxc1a9Ys0fcIDQ3VL04p0H0uHbrPpUP3ufToXpeOkrjPl+qxOUcFxSIiIuJRlNyIiIiIR1Fy40b+/v6MHTsWf39/q0PxaLrPpUP3uXToPpce3evSURbuc4UrKBYRERHPpp4bERER8ShKbkRERMSjKLkRERERj6LkRkRERDyKkptimjFjBvXq1SMgIIB27dqxevXqi16/atUq2rVrR0BAAPXr12fWrFmlFGn5Vpz7vHjxYnr06EF0dDShoaF07NiRL7/8shSjLb+K+/N8ztq1a/Hx8aF169YlG6CHKO59zsnJ4ZlnnqFOnTr4+/vToEED5syZU0rRll/Fvc/z58+nVatWBAUFUa1aNYYOHUpaWlopRVs+JSQk0LdvX6pXr47NZmPp0qWXfI4l34OGFNl///tfw9fX13j77beNHTt2GI8++qgRHBxs7N+/v9Dr9+7dawQFBRmPPvqosWPHDuPtt982fH19jYULF5Zy5OVLce/zo48+akycONH4/vvvjV27dhljxowxfH19jS1btpRy5OVLce/zOSdPnjTq169v9OzZ02jVqlXpBFuOXc59vuWWW4wOHToYK1euNPbt22ds3LjRWLt2bSlGXf4U9z6vXr3a8PLyMl577TVj7969xurVq41mzZoZ/fv3L+XIy5dly5YZzzzzjLFo0SIDMJYsWXLR6636HlRyUwxXX321MXLkSJdjjRs3Np566qlCr3/yySeNxo0buxx74IEHjGuuuabEYvQExb3PhWnatKkxbtw4d4fmUS73Pg8ePNh49tlnjbFjxyq5KYLi3ucvvvjCCAsLM9LS0kojPI9R3Ps8adIko379+i7HXn/9daNmzZolFqOnKUpyY9X3oIaliig3N5fNmzfTs2dPl+M9e/Zk3bp1hT5n/fr1Ba7v1asXmzZtIi8vr8RiLc8u5z7/mcPhIDMzk4iIiJII0SNc7n2eO3cue/bsYezYsSUdoke4nPv86aef0r59e1555RVq1KhBXFwcjz/+OKdPny6NkMuly7nPnTp14uDBgyxbtgzDMDh69CgLFy7kpptuKo2QKwyrvgcr3MaZlys1NRW73U5MTIzL8ZiYGI4cOVLoc44cOVLo9fn5+aSmplKtWrUSi7e8upz7/GeTJ08mKyuL22+/vSRC9AiXc58TExN56qmnWL16NT4++tNRFJdzn/fu3cuaNWsICAhgyZIlpKam8uCDD3L8+HHV3VzA5dznTp06MX/+fAYPHsyZM2fIz8/nlltuYdq0aaURcoVh1fegem6KyWazubQNwyhw7FLXF3ZcXBX3Pp/z4Ycf8vzzz7NgwQKqVKlSUuF5jKLeZ7vdzp133sm4ceOIi4srrfA8RnF+nh0OBzabjfnz53P11VfTp08fpkyZwrx589R7cwnFuc87duzgkUce4bnnnmPz5s0sX76cffv2MXLkyNIItUKx4ntQ//wqoqioKLy9vQv8K+DYsWMFstJzqlatWuj1Pj4+REZGllis5dnl3OdzFixYwPDhw/n444/p3r17SYZZ7hX3PmdmZrJp0ya2bt3Kww8/DJhfwoZh4OPjw4oVK7j++utLJfby5HJ+nqtVq0aNGjUICwtzHmvSpAmGYXDw4EFiY2NLNOby6HLu84QJE7j22mt54oknAGjZsiXBwcF06dKFl156ST3rbmLV96B6borIz8+Pdu3asXLlSpfjK1eupFOnToU+p2PHjgWuX7FiBe3bt8fX17fEYi3PLuc+g9ljM2TIED744AONmRdBce9zaGgoP//8M9u2bXM+Ro4cSaNGjdi2bRsdOnQordDLlcv5eb722ms5fPgwp06dch7btWsXXl5e1KxZs0TjLa8u5z5nZ2fj5eX6Fejt7Q2c71mQK2fZ92CJlit7mHNTDWfPnm3s2LHDGDVqlBEcHGz8/vvvhmEYxlNPPWXcfffdzuvPTYF77LHHjB07dhizZ8/WVPAiKO59/uCDDwwfHx9j+vTpRnJysvNx8uRJqz5CuVDc+/xnmi1VNMW9z5mZmUbNmjWNQYMGGb/88ouxatUqIzY21hgxYoRVH6FcKO59njt3ruHj42PMmDHD2LNnj7FmzRqjffv2xtVXX23VRygXMjMzja1btxpbt241AGPKlCnG1q1bnVPuy8r3oJKbYpo+fbpRp04dw8/Pz2jbtq2xatUq57l7773X6Nq1q8v13333ndGmTRvDz8/PqFu3rjFz5sxSjrh8Ks597tq1qwEUeNx7772lH3g5U9yf5z9SclN0xb3PO3fuNLp3724EBgYaNWvWNEaPHm1kZ2eXctTlT3Hv8+uvv240bdrUCAwMNKpVq2bcddddxsGDB0s56vLl22+/vejf27LyPWgzDPW/iYiIiOdQzY2IiIh4FCU3IiIi4lGU3IiIiIhHUXIjIiIiHkXJjYiIiHgUJTciIiLiUZTciIiIiEdRciMiIiJukZCQQN++falevTo2m42lS5cW+zW+/PJLrrnmGkJCQoiOjubWW29l3759xXoNJTciIn8wZMgQ+vfvb3UYIuVSVlYWrVq14o033ris5+/du5d+/fpx/fXXs23bNr788ktSU1MZOHBgsV5HKxSLiNutW7eOLl260KNHD5YvX251OMWSnp6OYRiEh4dbHYpIuWaz2ViyZInLPxZyc3N59tlnmT9/PidPnqR58+ZMnDiRbt26AbBw4UL+8pe/kJOT49zY9LPPPqNfv37k5OQUebNN9dyIiNvNmTOHv//976xZs4akpKSLXmsYBvn5+QWO5+bmllR4FxUWFqbERqSEDB06lLVr1/Lf//6Xn376idtuu40bb7yRxMREANq3b4+3tzdz587FbreTnp7O+++/T8+ePYu1i7iSGxFxq6ysLD766CP+9re/cfPNNzNv3jyX89999x02m40vv/yS9u3b4+/vz+rVq+nWrRsPP/wwo0ePJioqih49egAwZcoUWrRoQXBwMLVq1eLBBx/k1KlTzvcKDQ1l4cKFLu/x2WefERwcTGZmZqExLly4kBYtWhAYGEhkZCTdu3cnKysLcB2W+v3337HZbAUe5/6VCWYvVXx8PIGBgdSqVYtHHnnE+Voict6ePXv48MMP+fjjj+nSpQsNGjTg8ccfp3PnzsydOxeAunXrsmLFCp5++mn8/f0JDw/n4MGD/Pe//y3Weym5ERG3WrBgAY0aNaJRo0b89a9/Ze7cuRQ2+v3kk08yYcIEdu7cScuWLQF499138fHxYe3atbz55psAeHl58frrr7N9+3beffddvvnmG5588kkAgoODueOOO5x/GM+ZO3cugwYNIiQkpMD7Jicn85e//IVhw4axc+dOvvvuOwYOHFhojLVq1SI5Odn52Lp1K5GRkcTHxwPw888/06tXLwYOHMhPP/3EggULWLNmDQ8//PCV3UQRD7RlyxYMwyAuLo5KlSo5H6tWrWLPnj0AHDlyhBEjRnDvvffyww8/sGrVKvz8/Bg0aFChv6MXVOL7jotIhdKpUydj6tSphmEYRl5enhEVFWWsXLnSef7bb781AGPp0qUuz+vatavRunXrS77+Rx99ZERGRjrbGzduNLy9vY1Dhw4ZhmEYKSkphq+vr/Hdd98V+vzNmzcbgPH7778Xev7ee+81+vXrV+D46dOnjQ4dOhg333yzYbfbDcMwjLvvvtu4//77Xa5bvXq14eXlZZw+ffqSn0XEkwHGkiVLnO3//ve/hre3t/Hrr78aiYmJLo/k5GTDMAzj2WefNdq1a+fyOgcOHDAAY/369UV+b/XciIjb/Pbbb3z//ffccccdAPj4+DB48GDmzJlT4Nr27dsX6di3335Ljx49qFGjBiEhIdxzzz2kpaU5h36uvvpqmjVrxnvvvQfA+++/T+3atZ29K3/WqlUrbrjhBlq0aMFtt93G22+/zYkTJy752YYPH05mZiYffPCBs9Bx8+bNzJs3z+Vfob169cLhcBR76qqIp2vTpg12u51jx47RsGFDl0fVqlUByM7Oxtvb2+V559oOh6PI76XkRkTcZvbs2eTn51OjRg18fHzw8fFh5syZLF68uEACERwcXOD5fz62f/9++vTpQ/PmzVm0aBGbN29m+vTpAOTl5TmvGzFihHNoau7cuQwdOhSbzVZojN7e3qxcuZIvvviCpk2bMm3aNBo1anTRZOSll15i+fLlfPrppy5DXQ6HgwceeIBt27Y5Hz/++COJiYk0aNDgEndLxPOcOnXK+bsAsG/fPrZt20ZSUhJxcXHcdddd3HPPPSxevJh9+/bxww8/MHHiRJYtWwbATTfdxA8//MALL7xAYmIiW7ZsYejQodSpU4c2bdoUPRC39D2JSIWXl5dnxMTEGJMnTzZ+/vlnl0dcXJwxbdo0wzDOD0udOHHC5fldu3Y1Hn30UZdjCxcuNHx8fJzDQIZhGC+++GKB5x8/ftwICAgwXnvtNcPLy8s4cOBAkePOz883atSoYUyePNkwjILDUgsXLjR8fX2Nr776qsBz77zzTuP6668v8nuJeLpzv99/ftx7772GYRhGbm6u8dxzzxl169Y1fH19japVqxoDBgwwfvrpJ+drfPjhh0abNm2M4OBgIzo62rjllluMnTt3FisOH3dmbCJScX3++eecOHGC4cOHExYW5nJu0KBBzJ49u9iFtg0aNCA/P59p06bRt29f1q5dy6xZswpcV7lyZQYOHMgTTzxBz549qVmz5gVfc+PGjXz99df07NmTKlWqsHHjRlJSUmjSpEmBa7dv384999zDP//5T5o1a8aRI0cA8PPzIyIign/+859cc801PPTQQ9x3330EBwezc+dOVq5cybRp04r1WUU8Qbdu3S5a+Ovr68u4ceMYN27cBa+54447nEPbl0vDUiLiFrNnz6Z79+4FEhuAW2+9lW3btrFly5ZivWbr1q2ZMmUKEydOpHnz5syfP58JEyYUeu3w4cPJzc1l2LBhF33N0NBQEhIS6NOnD3FxcTz77LNMnjyZ3r17F7h206ZNZGdn89JLL1GtWjXn49xqqS1btmTVqlUkJibSpUsX2rRpw//93/9RrVq1Yn1OEXEvrVAsIh5h/vz5PProoxw+fBg/Pz+rwxERC2lYSkTKtezsbPbt28eECRN44IEHlNiIiIalRKR8e+WVV2jdujUxMTGMGTPG6nBEpAzQsJSIiIh4FPXciIiIiEdRciMiIiIeRcmNiIiIeBQlNyIiIuJRlNyIiIiIR1FyIyIiIh5FyY2IiIh4FCU3IiIi4lH+H+O8Z0gQVM0iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(size_list, numba_cpu_times, label=\"Numba CPU\")\n",
    "plt.plot(size_list, numpy_cpu_times, label=\"NumPy\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Array size\")\n",
    "plt.ylabel(\"Time elapsed\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de3c56d9-a450-4070-91cc-213f6baa10b1",
   "metadata": {},
   "source": [
    "Next let's walk through a `guvectorize` example, calculating a moving average. In this case we are illustrating Eager/decoration-time compilation where we pass the type signatures. In general, with `guvectorize` we can work on an arbitrary number of elements of input arrays (as opposed to one element at a time with `vectorize`).\n",
    "\n",
    "The layout in symbolic form `(n),()->(n)`is how we tell NumPy that we take in an input array of shape `n` and a constant which is represented by the `()`. Then the output of this function will also be of the same shape as our original input array `(n)`.\n",
    "\n",
    "Also note in guvectorize we don't have a return, but write our results to one of the arguments, in this case `out`. This is particularly useful because with guvectorize we can have return arrays of differing dimensions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "16d9fd71-2d1e-4a99-a57a-6c1ba03e863b",
   "metadata": {},
   "outputs": [],
   "source": [
    "@guvectorize(['void(float64[:], intp[:], float64[:])'], '(n),()->(n)')\n",
    "def move_mean(a, window_arr, out):\n",
    "    window_width = window_arr[0]\n",
    "    asum = 0.0\n",
    "    count = 0\n",
    "    for i in range(window_width):\n",
    "        asum += a[i]\n",
    "        count += 1\n",
    "        out[i] = asum / count\n",
    "    for i in range(window_width, len(a)):\n",
    "        asum += a[i] - a[i - window_width]\n",
    "        out[i] = asum / count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d553c621-e900-46f6-901d-26246e040fdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input array:\n",
      "[[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9.]\n",
      " [10. 11. 12. 13. 14. 15. 16. 17. 18. 19.]\n",
      " [20. 21. 22. 23. 24. 25. 26. 27. 28. 29.]]\n",
      "Output array:\n",
      "[[ 0.   0.5  1.   2.   3.   4.   5.   6.   7.   8. ]\n",
      " [10.  10.5 11.  12.  13.  14.  15.  16.  17.  18. ]\n",
      " [20.  20.5 21.  22.  23.  24.  25.  26.  27.  28. ]]\n"
     ]
    }
   ],
   "source": [
    "arr = np.arange(30, dtype=np.float64).reshape(3,10)\n",
    "print(\"Input array:\")\n",
    "print(arr)\n",
    "print(\"Output array:\")\n",
    "print(move_mean(arr, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2549ecff-6194-423d-86c6-3577a5d92ae4",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Custom Numba Kernels: Interoperability**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "746fa39d-9742-4241-a68e-b615a097ba8c",
   "metadata": {},
   "source": [
    "In general, we recommend using CuPy arrays when writing Numba code targeting the GPU. Below in our L2 norm calculation we will use CuPy arrays with our Numba kernels.\n",
    "\n",
    "Our first `sum_reduce` function uses the `reduce` decorator which takes our binary summation operation and turns it into a reduction kernel, running one the GPU as indicated by the `cuda` in our decorator. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e287d7d9-7209-491e-9fff-aad83cfef206",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.reduce\n",
    "def sum_reduce(a, b):\n",
    "    return a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "938e2eee-3fe1-4193-9330-9cba48d11596",
   "metadata": {},
   "source": [
    "The next kernel is our first just-in-time kernel targeting the GPU. We use:\n",
    "\n",
    "`cuda.grid` The absolute position of the current thread in the entire grid of blocks \\\n",
    "`cuda.gridsize` The size in threads of the entire grid of blocks\n",
    "\n",
    "We can safely loop over the data using one thread per data element, without assuming we have enough threads to cover the whole array size. Instead we loop over our data one grid size at a time (aka one stride)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "42c959f6-484f-4ebc-9604-ba56c21ea427",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def numba_l2_norm(x):\n",
    "    start = cuda.grid(1)\n",
    "    stride = cuda.gridsize(1)\n",
    "\n",
    "    for i in range(start, x.shape[0], stride):\n",
    "        x[i] = x[i] * x[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccaefc23-f17d-4ed8-879b-a885c575d5a1",
   "metadata": {},
   "source": [
    "Go ahead and change our problem size `x` to see how varying the array size affects both NumPy and Numba performance. Also explore how we can adjust our  <code>threads_per_block</code> to get the best utilization of the GPU.\n",
    "\n",
    "Note, if we weren't comparing the NumPy implementation, we could create our d_x directly using <code>d_x = cp.random.rand(2**26)</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bbf2fbcb-e03b-46aa-a67d-e5835b806818",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.rand(2**22)\n",
    "d_x =cp.array(x)\n",
    "threads_per_block = 1024\n",
    "blocks_per_grid = (d_x.size + (threads_per_block - 1)) // threads_per_block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a9ec8a0c-deec-4dd5-b710-aff3d0dbb07a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/rapids/lib/python3.9/site-packages/numba/cuda/dispatcher.py:488: NumbaPerformanceWarning: \u001b[1mGrid size 64 will likely result in GPU under-utilization due to low occupancy.\u001b[0m\n",
      "  warn(NumbaPerformanceWarning(msg))\n",
      "/opt/conda/envs/rapids/lib/python3.9/site-packages/numba/cuda/dispatcher.py:488: NumbaPerformanceWarning: \u001b[1mGrid size 1 will likely result in GPU under-utilization due to low occupancy.\u001b[0m\n",
      "  warn(NumbaPerformanceWarning(msg))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numba time elapsed (fractional seconds):  4.711460737511516\n",
      "Numba output:  1181.995920294273\n",
      "NumPy time elapsed (fractional seconds): 0.016373537480831146\n",
      "NumPy output:  1181.995920294273\n"
     ]
    }
   ],
   "source": [
    "start_time = time.monotonic()\n",
    "numba_l2_norm[blocks_per_grid, threads_per_block](d_x)\n",
    "output = cp.sqrt(sum_reduce(d_x))\n",
    "print(\"Numba time elapsed (fractional seconds): \", time.monotonic() - start_time)\n",
    "print(\"Numba output: \", output)\n",
    "\n",
    "start_time_np = time.monotonic()\n",
    "np.linalg.norm(x, 1)\n",
    "print(\"NumPy time elapsed (fractional seconds):\" , time.monotonic() - start_time_np)\n",
    "print(\"NumPy output: \", output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c0f4a20-496f-452f-8548-b682d0e39649",
   "metadata": {},
   "source": [
    "You will see the first couple of calls to the Numba kernel are slower than the NumPy implementation, but after a couple reruns, the performance levels out and Numba is much faster.\n",
    "\n",
    "This example also shows us how the GPU and specifically custom Numba kernels as we have written will perform well with very large input data (it also does better when the problem is more complicated than we have demonstrated above). To show this point, when you decrease the input size, the NumPy implementation will greatly outperform our solution. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1602d31e-74d0-47e5-84ca-7e61b2243bf6",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Numba custom kernels: SAXPY**\n",
    "\n",
    "Below we have a simple SAXPY kernel that utilizes the `vectorize` decorator. We specify the expected types therefore doing eager compilation, and make sure to set the target to run on the GPU. You can see when we are calling a function with the decorator it is different from how we call our <code>@cuda.jit</code> functions. We do not have to pass our <code>[blocks_per_grid, threads_per_block]</code>. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0aff03c6-31c1-4a08-9458-47a566e4b775",
   "metadata": {},
   "outputs": [],
   "source": [
    "@vectorize(['float32(float32, float32, float32)'], target='cuda')\n",
    "def saxpy(a, x, y):\n",
    "    return a * x + y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0f13a976-1d1c-47f0-95d5-84655d9ebf75",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 172 ms, sys: 68.1 ms, total: 240 ms\n",
      "Wall time: 260 ms\n"
     ]
    }
   ],
   "source": [
    "N = 2**25\n",
    "A=cp.random.randn(N).astype(cp.float32)\n",
    "B=cp.random.randn(N).astype(cp.float32)\n",
    "\n",
    "%time C = saxpy(2.0, A, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52e2981e-6ffe-4670-8f83-e5c371968e90",
   "metadata": {},
   "source": [
    "Experiment with differnt array sizes by changing `N` and seeing how the performance of our Numba timing changes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0aa9223-7aa9-42f8-9844-dc5c5d1d3ce7",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Numba custom kernels: fast matrix multiplication**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fce2277-9862-4ed5-85a8-6ff34d6fd7da",
   "metadata": {},
   "source": [
    "In this final section of our Numba introduction we will use a common example of GPU performance - matrix multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c53a2174-a474-423e-bd31-aae382245213",
   "metadata": {},
   "outputs": [],
   "source": [
    "threads_per_block = 16\n",
    "blocks_per_grid = 128"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f267578-dbbf-4cdd-b08d-8ccd424fcf15",
   "metadata": {},
   "source": [
    "If you are new to CUDA concepts or thread hierarchy, keep the threads per block and BPG blocks per grid the same, but if you are curious, test out how changing that value affects your performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ad3e7d2c-35fe-4aca-a8d6-7abe582d5106",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit\n",
    "def fast_matmul(A, B, C):\n",
    "    sA = cuda.shared.array(shape=(threads_per_block, threads_per_block), dtype=float32)\n",
    "    sB = cuda.shared.array(shape=(threads_per_block, threads_per_block), dtype=float32)\n",
    "\n",
    "    x, y = cuda.grid(2)\n",
    "\n",
    "    tx = cuda.threadIdx.x\n",
    "    ty = cuda.threadIdx.y\n",
    "    blocks_per_grid = cuda.gridDim.x\n",
    "\n",
    "    tmp = float32(0.)\n",
    "    for i in range(blocks_per_grid):\n",
    "        sA[ty, tx] = 0\n",
    "        sB[ty, tx] = 0\n",
    "        if y < A.shape[0] and (tx + i * threads_per_block) < A.shape[1]:\n",
    "            sA[ty, tx] = A[y, tx + i * threads_per_block]\n",
    "        if x < B.shape[1] and (ty + i * threads_per_block) < B.shape[0]:\n",
    "            sB[ty, tx] = B[ty + i * threads_per_block, x]\n",
    "        cuda.syncthreads()\n",
    "\n",
    "        for j in range(threads_per_block):\n",
    "            tmp += sA[ty, j] * sB[j, tx]\n",
    "            \n",
    "        cuda.syncthreads()\n",
    "    if y < C.shape[0] and x < C.shape[1]:\n",
    "        C[y, x] = tmp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4f1a1f6-ccc6-46bc-8894-8b31218f16a7",
   "metadata": {},
   "source": [
    "Now we will show how to call and time this fast matrix multiplication kernel we have written."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "fdd5eeaf-6e1b-4430-9c35-aabd00044171",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 432 ms, sys: 19.7 ms, total: 452 ms\n",
      "Wall time: 475 ms\n"
     ]
    }
   ],
   "source": [
    "x_d = cp.arange(1000).reshape([50,20])\n",
    "y_d = cp.ones([50,20])\n",
    "z_d = cp.zeros([50,20])\n",
    "\n",
    "%time fast_matmul[blocks_per_grid, threads_per_block](x_d, y_d, z_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "240dbdd3-224b-4980-bc40-d347077a1a81",
   "metadata": {},
   "source": [
    "Below we call our Numba kernel on NumPy arrays, see on the first pass we should get a warning about the conversion to a device array happening. The warning will go go away once you have run the cells again. Go ahead and run the cell above and below again, notice the changes in performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "56ec5bf3-71ea-443a-a6be-30327d033dec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.35 ms, sys: 0 ns, total: 2.35 ms\n",
      "Wall time: 2.35 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/rapids/lib/python3.9/site-packages/numba/cuda/cudadrv/devicearray.py:885: NumbaPerformanceWarning: \u001b[1mHost array used in CUDA kernel will incur copy overhead to/from device.\u001b[0m\n",
      "  warn(NumbaPerformanceWarning(msg))\n"
     ]
    }
   ],
   "source": [
    "x_d = np.arange(1000).reshape([50,20])\n",
    "y_d = np.ones([50,20])\n",
    "z_d = np.zeros([50,20])\n",
    "\n",
    "%time fast_matmul[blocks_per_grid, threads_per_block](x_d, y_d, z_d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63dd0b27-a10b-4c80-a0ae-883d9dd6749c",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Please restart the kernel**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44058450-b452-40eb-9791-c5c2096a6c40",
   "metadata": {},
   "outputs": [],
   "source": [
    "import IPython\n",
    "app = IPython.Application.instance()\n",
    "app.kernel.do_shutdown(True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}