{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import time\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['text.usetex'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MA934\n",
    "\n",
    "## Divide and conquer algorithms\n",
    "\n",
    "### Estimating the computational complexity of numerical computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Big O notation\n",
    "\n",
    "Given functions, $f(x)$ and $g(x)$, we say that $f(x) = \\mathcal{O}(g(x))$ if there exist  $x_0$ and $K > 0$ such that\n",
    "\n",
    "$$ \\left| f(x) \\right| \\leq K\\,g(x) $$\n",
    "\n",
    "for all $x> x_0$.  \n",
    "Used to describe an upper bound on the growth of a function, $f(x$), as its argument, $x$, tends to $\\infty$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### FLOPS and computational complexity \n",
    "\n",
    "FLOPS : Floating point operations per second\n",
    "\n",
    "The *Computational complexity* of an algorithm is the rate at which the total number of operations (additions, multiplications, comparisons) required to solve a problem of size $n$ grows with $n$. \n",
    "\n",
    "We denote this number of operations by $F(n)$. \n",
    "\n",
    "The more efficient the algorithm, the slower the $F(n)$ grows with $n$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example : standard matrix multiplication\n",
    "\n",
    "Consider calculating $C = A \\, B$  where $A$ and $B$ are $n\\times n$ matrices.  \n",
    "Using the standard algorithm, the $i\\,j$ entry of C is:\n",
    "\n",
    "$$c_{i\\,j} = \\sum_{k=1}^n a_{i\\,k}b_{k\\,j}.$$\n",
    "\n",
    "Needs $n$ multiplications and $n-1$ additions,  giving $2n-1$ ops in total. Since there are $n^2$ entries, the total number of ops required is $n^2\\,(2 n -1)$.\n",
    "\n",
    "Thus matrix multiplication has computational complexity $F(n) = \\mathcal{O}(n^3)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Aside : Benchmarking\n",
    "\n",
    "We are ofted interested in accurately recording execution times and benchmarking in general. Let's check matrix multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import timeit\n",
    "\n",
    "# function that computes the standard matrix-matrix multiplication\n",
    "def MMStd(A, B):\n",
    "    # Detect matrix size\n",
    "    n = len(A)\n",
    "    nSize = (n, n)\n",
    "    # Initialise result matrix\n",
    "    C = np.zeros(nSize)\n",
    "    \n",
    "    # Loop through the matrix elements\n",
    "    for i in range(n):\n",
    "        for j in range(n):     \n",
    "            for k in range(n): \n",
    "                # Use standard multiplication formula\n",
    "                C[i][j] += A[i][k]*B[k][j]\n",
    "    return C\n",
    "\n",
    "# This gives us the list of powers of two up to a specified amount\n",
    "sizes = 2 ** np.arange(9)\n",
    "\n",
    "# Initialise the list of runtimes with zeros\n",
    "runtimes = [np.float32(0.0)] * 9\n",
    "\n",
    "# Initialise an iterator to zero as well\n",
    "j=0\n",
    "\n",
    "for i in sizes:\n",
    "    A = np.random.rand(i, i)\n",
    "    B = np.random.rand(i, i)\n",
    "    \n",
    "    starttime = timeit.default_timer()\n",
    "    CStd = MMStd(A, B)\n",
    "    runtimes[j] = timeit.default_timer() - starttime\n",
    "    j=j+1\n",
    "    \n",
    "# Exercise: what could go wrong with the above method? How would you improve it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAHzCAYAAABPMYkzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABmQElEQVR4nO3deXhMZ//H8c8RS1EaqkUtJX59hKJEWlVaVIJuVGuprtaEWGoPqlprJPatJLSqK5Wiq1aCtFWPNai9CFWK2mINIbl/f2SSR0hkkWRmkvfrunI1mXOfM9+JHj5zz/fcxzLGCAAAAIBzymfvAgAAAABkHoEeAAAAcGIEegAAAMCJEegBAAAAJ0agBwAAAJwYgR4AAABwYvntXYCzK1WqlKlUqVKa4y5duqSiRYtmf0EAsgznLeBcOGeRm23evPmUMea+lLYR6O9QpUqVtGnTpjTHRUREqHHjxtlfEIAsw3kLOBfOWeRmlmX9ldo2Wm4AAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAnRqAHAAAAnBiBHgAAAHBiBHoAAADAiRHoAQAAACfGnWJz2NWrV3XmzBlduHBBcXFx9i4HwG0UL15cx44dU8mSJVWoUCF7lwMAQIoI9Dno6tWrOnz4sEqUKKFKlSqpQIECsizL3mUBSIExRtHR0YqLi9Phw4dVsWJFQj0AwCER6HPQmTNnVKJECZUqVcrepQBIg2VZyp8/v0qUKCEp4fwtW7asnasCAOBW9NDnoAsXLqh48eL2LgNABhUvXlwXLlywdxkAAKSIQJ+D4uLiVKBAAXuXASCDChQowDUvAACHRaDPYfTMA86H8xYA4MgI9AAAAIATI9ADAAAAToxADwAAADgxAj0cTkhIiOrWrasSJUrIsixVqVJFbdu2VXh4uL1Lc3je3t6qUqWKvcvItKCgIFmWxZ81AAAZwDr0cCje3t4KDw+Xl5eX2rdvL0k6cOCAvvrqK0mSl5eXPctDNnN1dZWbm5u9ywAAwKkQ6OEwqlSpoqioKAUHB8vHxyfZtuDgYEVFRdmpsuwREhIiNze3PPkmJbXX7uPjc8ufPQAAuD1abnKZMuUryrKsNL/KlK9o71KTCQoKSjXMJ8ptM7f+/v5avHixvcuwi7z82gEAyGrM0OcyJ47+LfVfmfa4yU1zoJr0iY6Olr+/v9zc3JidBQAAyCBm6GF3mzZtkiT5+vpmaL/o6Gj5+vqqSpUqKlGihNq2bXtLW05oaKhKlCihyMhI+fr6qkSJEipRokTSc0VFRcnb2zvp4tuQkJBk+4eHh6tu3bqKjIxUUFCQ6tatK8uy5O3trejo6GRjQ0JCZFmWIiMjkz2eeIFvorZt28qyLEVHRyftY1mWgoKCkr22tm3bqkSJEqpSpYr8/f1vef2JtZUoUUJ169ZVaGiozpw5k6HfYeLvJyoqSuHh4fL29k763aT39dx4nMjISPn7+yf7M7lRWq89pee80z/D9P4+AQBwVgR62F1ieMtIS01UVJQqV66s8PBw+fr6KjAwUFFRUapSpUqyFVLOnDmj6OhoNW2a8IlEYGCgPD09FRISorZt28rb21ve3t4KDAzUmTNn5OvrmyxMRkdHKzIyUnXr1tXGjRvVvn17eXh4JIXpzAgMDFRYWJikhIt8w8LCFBYWlvTpROJri4yMVGBgoHx9fRUUFJTsDU9kZKS8vb0VGRkpHx8f+fr6KiAg4JbwnR7R0dEKDg6Wt7e3zpw5I29v7wwfI/H33LZt26RPXDw9PRUaGposPKf12m937Mz+Gabn9wkAgDOj5cZB9O3bV1u3bs3R52zcuPEd7V+7dm1NnTr1jus4cOCApIwF+sQwlrivlHBBZd26deXr65vscSkhPAYHByeNK1GihEJDQ7V48WK1adNGkuTh4ZG0yo6Hh0ey/W/s7R88eHDSuMQVeTLCzc0t6bWmdGGor6+voqOjdfDgQbm6ukpKWP0l8Y2Lq6urunXrlvT6E4/l4+OT6SUrg4KCkv0uMsvNzS3Z79myLIWGhiowMDBp++1e++1k9s8wPb9PAACcGTP0cBjpbRdJbA8ZOnToLdsSZ+pDQ0OTPZ64BGaixFB5Y4D19PSUJJ0+ffqW4978ZiMxWGb1hZ3R0dEKDw9XmzZtdObMGUVFRSkqKirp+cPDw5M+NfDx8bmlrsyG0zZt2txxmJd0y+y+h4dHhtuAUpOZP8P0/D4BAHB2zNA7iKyY6ZYky7LSPTYiIiJLnvNOJbauREZGpmvGNrGd4uZZdOl/gW7jxo3Jgt7NQbdkyZJ3NDObGAizeinNxOsJQkNDb3lTkvh8ic+ZlTeQujksZ1ZKfyZZJTN/hun5fQIA4OyYoYfdtWvXTtL/Zr3T6+aLUqXUZ/lLliyZ4brsKSwsTMaYW74GDx6cNCYrW0WyaknQ7Pw938mx0/P7BABkjLMulZ0bMUMPu3N1ddXgwYMVFBSk0NDQNFs/EmeBb56Fl/43e//oo49mT7E2iTO7dzojffMbkMRPGG73aUVi+N68efMdPbe9ZVUrzu2k5/cJAMgcZ1wqO7dihh4OITAwUG5ubmrbtm2KrRHS//qd3dzc5OHhoZCQkFtm6f39/eXq6pol/eA3urk1I/Gi3JRaVW4Mqon97ilxdXW9pX5XV1d5eXkpICDglm3R0dGKjo6Wq6triq8/Ojo6W1pI0vt6MiKl154d0vP7BADA2TFDn8uULlchXe+ES5erkAPVZMzmzZvVtGlTtW3bVm3atJG3t7dKliyZdAfZqKgoGWMkJVyMWrduXVWuXDnp4tjEMYnLImalxJVSXF1dFRwcrMjISLVp0ybZDH3izHniii6JyzemxtPTU+Hh4QoKCkpalSc4OFjBwcHJXpurq6s2b96skJCQpBVdAgMD5e3trcqVKyc9X2BgoKKjo7Os7SWjrycjUnvt2SE9v08AAJwZgT6XOX7ksL1LyLTEoBUUFKTg4OBkM/Vt2rRJtqKMm5ubDh48qG7duikgIEBSQkgMCwvLsn7wGy1evFhhYWH66quvJCUsXZkYdBN5eXmpTZs2Cg0N1aZNm+Tp6ang4GAtXrw4xZlzf39/bdq0SQEBAfL09EwKy4mvzd/fP+lNioeHh4KDg5PCZ+Ia7v7+/vL3909a0jE4ODjLZp0z+noyIrXXnh3S8/sEAMCZWYkznsgcT09Pk7iSxu1ERESodOnSqlatWg5UhawSGhqqtm3bKiwsjB7sPOjChQsqVqyYJGn37t2cv4CDi4iIuON7rCD9LMtKVw+9JjcVefPOWZa12RjjmdI2eugBAABwZ+KuS9di7F1FnkWgt7Esy82yrGDLsvgMHgAAIBUxMTEaPXr0/x64Hit93kNa86H9isrj6KGXZFlWYi+FmyTnWrAcAAAgBxhjtHjxYg0aNEiHD99wzV7+glLVxlKpynarLa9jhl6SMSbcGBMuKdretQAAADiaLVu2qFGjRmrfvr1cXV3Vv3//hA1HdyT8t95rUpUn7FdgHscMPXAbbdq04UIeAECedeLECQ0fPlwffvihSpUqpXfeeUcrVqzQ5MmTVaBgIV1b9Haax3DEpbJzG6cP9Lae9/bGmLapbHeVNFTSAdtDVYwx2bdGHgAAgJOLjY3V9OnTNWrUKMXExKh///4yxiggIEClSpXSZ599pldffTVhpRvYndMGesuyEu9Ck1bf+2JJvsaYKNt+bpZlhRljvLO7RgAAAGdijNH333+vAQMGaN++fXr22Wc1adIkubu7a+rUqerevbvGjh0rV1dXe5eKGzhtD70xxtcY46uEwJ4i2+x9VGKYt+0XdcM2AAAASNq1a5datGihli1bysXFRXPnzlVsbKwiIyMlSX379tWsWbMI8w7IaQN9OrWXtDmFx8Mk+eZwLQAAAA7nzJkz6tOnj2rVqqUNGzZowoQJat26tXr27KmNGzfq+vXr9i4RaXDalpt08pIUnMLjUZJSvNMWAABAXnD9+nUFBwdrxIgRio6Olq+vr7y8vDRw4EAdPHhQr7/+uiZOnKjSpUvbu1SkIdfO0NsuhnWVdCaFzdG2bQAAAHlOeHi4ateurV69eql27draunWrPvjgA0nSXXfdpdWrV+vTTz8lzDuJXBvolY4bRNlCvyzL8rAsa7ASZvR9LcsanLgNAAAgt9i/f79efPFFeXt7KyYmRosXL1aLFi0UFhYmSWrdurW2bdumxo0b27dQZEhub7lJF2NMpKRISUHpGW9Zlo8kH0kqXbq0IiIi0tzn4sWLuueee3ThwoU7qBQ3io6O1tKlSyVJW7duVZMmTfTiiy/atyjkKnFxcUnn7JUrV9J1rgOwn4sXL3KepuLSpUv67LPP9PXXXyt//vzq1q2bqlatqkGDBunQoUNq3Lix6tSpwzKUTopAnwnGmBBJIZLk6elp0vMuNiIiQnfddZeKFSuWzdXlHQMHDlRw8P8ukbAsS5s3b5aHh4cdq0JucuHChaRz9q677lKdOnXsXBGA24mIiGBm+Sbx8fH6+OOPNWzYMJ04cUIdO3bUgAEDNHHiRA0cOFAVK1bUN998o5YtW9q7VNyB3Nxyk1LvfDLGmOgcqAPZZNOmTQoPD0/62dXVVVFRUbfZAwCAvOP333/XY489pi5dusjNzU0bNmzQ/PnzdeXKFS1cuFBDhgzRrl27CPO5QK6doTfGRFuWFa2EG09F3rTZTQkXxsKJbd78vxVJo6OjFR0dzew8ACDPO3z4sPz9/bVw4UKVK1dOn3/+udzd3bVy5Uo9+uij8vT01OHDh3X//ffbu1Rkkdw8Qy9Jm5TyxbFVJIWn8DiclL+/vxYvXiw3Nzd7lwIAgF1cvnxZI0eOlLu7u5YtW6YRI0Zo48aN2rBhgx599FFNnjxZ0dHRkkSYz2Vye6BfLMk7hce9JC3K4VqQDaKjoxUSEsJd6wAAeZYxRgsXLpS7u7vef/99vfDCC9q9e7eqV6+uunXravr06erevbt2797Nv5e5VG5ouXFVKmvKG2NCLMvytSzLzRgTJSUsUSnpjDEmNOdKRHZxdXWVj4+PJKlEiRKSpDZt2tizJAAAcszmzZv19ttv6/fff1edOnX0+eef68knn9SJEyfUpUsXubu765tvvtGjjz5q71KRjZx2ht6yrEDLsoIlDZXkZlnWYsuygm2B/UZNlbC2vI9tucn2xpiUZu3hRBJn5m/k5eWVbNUbAAByq+PHj6tLly569NFHtW/fPs2bN0+//vqr/vzzTxljVLp0aa1Zs0br168nzOcBTjtDb4zxt33rm8a4aEn+txsD57Np0yb5+/snzc4DAJAXXL16VdOmTdOYMWN05coVDRgwQMOHD9fatWv1yCOPKCoqSu7u7mrQoIFq165t73KRQ5w20CN3iYyMVHh4uMLCwhQcHCw3NzcFBSXc5+vAgQM6c+aMFi9enDTe09NTQ4cOTXaM8PDwZGMAAMgtjDH69ttvNWDAAB04cEAvvPCCJk2apLvuukudO3fWkiVL5O7urlWrVqlBgwb2Lhc5jEDvQFK6GUa7du3k5+eny5cv69lnn71le8eOHdWxY0edOnUqxd7xHj16qH379vr777/1xhtv3LJ9wIABeuGFF7R37175+t76Ycfw4cPl5eWlrVu3qm/fvsm2ZeXd+BYtWqTAwEAdOHBAvr6+8vDwUGBgYNL2KlWqKCQkJGlG3tXVVV5eXgoKCpKrq6s2b96suXPnysvLK8tqAgDAEezYsUP9+vVTeHi4qlevrp9//lnNmjVTfHy8Hn74YR06dEhjx47VwIEDVbBgQXuXCzsg0MPuIiMjk/r7oqKitGnTpltm2l1dXXXgwIFkj3l4eLDuPAAg1zp9+rTee+89zZkzR8WLF09arWbz5s2KjY1VwYIFNXfuXJUrV06VK1e2d7mwIwK9A7ndjHeRIkVuu71UqVK33V6hQoXbbq9ateptt9euXTtLZ+RvlvjpwqZNmzR06NBbltWKjIxM8RMEAABym2vXrmnOnDl67733dP78eXXv3l0jR45UfHy8fH19NX/+fE2ePFn9+vVTw4YN7V0uHIDTrnKD3CNxlj0qKkrR0dG3tM1ERibc6Jd2GgBAbhcWFqbatWurT58+8vDw0NatWzV9+vSkHvlPP/1UQ4YMYVEIJEOgh8MID0+4ee/NbTSLFi2Sm5sbd4EFAORa+/btU8uWLdWsWTNdvXpVy5YtU1hYmGrUqKEePXrIx8dHNWvW1LZt2xQQEKCiRYvau2Q4EAI9HEZYWFiKs/AhISFJLTlRUVE5XRYAANnm/PnzGjx4sB5++GGtXr1agYGB2rlzp5o0aaJz585Jkrp166ZPPvlEq1evVvXq1e1cMRwRgR4OIzw8XN7eye/5FRkZqejo6KT++dBQbvALAHB+cXFx+vDDD/XQQw9p4sSJeuONN7Rv3z4NGjRIy5Ytk7u7uwYPHiwpYanmN954Q5Zl2blqOCoCPRxCav3zZ86ckaurq9zc3BQVFUXbDQDA6f3222969NFH1bVrVz300EPauHGjPvzwQ124cEHNmzfXK6+8orJly6pr1672LhVOgkAPh5AY1m/un/fy8pKnp6eCgoIUHh6e4lr7AAA4g8OHD+uVV17RU089pZMnT+rLL7/Ub7/9prp16+qrr75SjRo1tH79es2YMUMbNmzQY489Zu+S4SRYthIOwcvL65Z15hOFhYXlcDUAAGSdS5cuKSgoSEFBQbIsS++//74GDRqkIkWKKCYmRoULF9bjjz+uV199VePGjVPZsmXtXTKcDIEeAAAgGxhj9OWXX2rw4ME6evSoOnTooMDAQFWoUEFHjhzRm2++qXPnzmnFihWqWLGi5s+fb++S4aRouQEAAMhiGzduVMOGDfXaa6+pTJky+u233/TFF1+oTJkymjRpktzd3fXDDz+oSZMmio+Pt3e5cHLM0AMAAGSRY8eOadiwYfr4449VunRpffTRR3rrrbeUL18+7du3Ty+//LK2b9+u5557TjNmzFDlypXtXTJyAQI9AADAHbpy5YqmTp2qsWPHKjY2Vv7+/ho2bJiKFy8uY4wkqUyZMrr77ru1dOlStWrVimUokWUI9AAAACkoU76iThz9O81xrvfer5L33K2oqCi1atVKEydO1P/93/8pPj5e8+bN04IFC7Ry5UoVK1ZMv//+O0EeWY5ADwAAkIITR/+W+q9Mc1z05KYqV+a+ZHc837Ztm3r06KH//ve/evLJJ3XmzBmVKVOGMI9sQaAHAAC4Q1u3blX+/PkVExOjYcOGafr06br33nu1YMEC7vKKbMcqNwAAAHcof/6EOdKCBQtqzZo18vHx0Z49e/Tmm28S5pHtmKHPYcYYTmzAySRe0AYAqenQoYNmzZqlkiVLas2aNSpUqJC9S0Iewgx9DnJxcdG1a9fsXQaADLp27ZpcXFzsXQYAR3M9Vlq7QJL0ww8/aOvWrZJEmEeOY4Y+BxUrVkznz59XqVKl7F0KgAw4f/68ihUrZu8yAOQQY4y+//772w86tFFaOV06948kae/evSpbtmwOVAfcihn6HFSyZEmdPXtWp06dUmxsLB/jAw7MGKPr16/r1KlTOnv2rEqWLGnvkgBkM2OMwsLC9Pjjj6tly5a3H7ztWylfPunlIEkizMOumKHPQYUKFVLFihV15swZHTp0SHFxcfYuCcBtxMTEqGzZsqpYsSIfoQO53K+//qp3331Xv/76qypWrKh58+apa9eu/xsQHydtWSq5PS6VKC81GygVKCzlL2i/ogEbAn0OK1SokMqWLcs7ecAJREREyMPDw95lAMhG69ev17vvvquwsDCVLVtWM2fOVNeuXVWoUKH/Bfp/dkrhU6VTUdLVS9ITb0mF77Fr3cCNCPQAACDP2bJli0aMGKHvv/9epUqV0qRJk9SjRw8VLlw4acx9ZR7QyclNk++47pOErxuULlchJ0oGUkUPPQAAyDN27dqltm3bysPDQ2vWrNHYsWN18OBB9e/fP1mYl6TOb72h/Pnza9CgQbpw4YKMMSl+HT9y2E6vBkjADD0AAMj19u3bp5EjR+qLL77Q3XffrREjRqhfv35ydXVNNu6PP/5QbGysPD09NWzYML3++uuqUaOGfYoG0okZegAAkGv99ddf6tq1q6pVq6YlS5Zo8ODBOnjwoEaOHJkszF+4cEEDBgyQh4eHBg4cKEkqXrw4YR5OgRl6AACQ6/zzzz8aO3as5s6dK8uy1KtXLw0ZMkRlypRJNs4Yo6+//lp9+/bV0aNH5ePjo3HjxtmpaiBzCPQAACDX+PfffzV+/HjNnj1b169fV9euXfXOO++ofPnyKY5funSp2rZtq9q1ays0NFSPP/54DlcM3DkCPQAAcHpnzpzR3Llz9c033ygmJkZvvvmmRowYocqVK98y9sqVK9qzZ49q166tli1b6uOPP9Zrr72m/PmJRXBO9NADAACndf78eY0cOVKVK1fWl19+qZYtW2rXrl2aP39+imF+xYoVqlmzppo1a6ZLly4pf/78euuttwjzcGoEegAA4HQuXbqkwMBAVa5cWe+//76aNm2qefPm6YsvvlDVqlVvGX/06FG1b99ezZs3l2VZ+vzzz1W0aFE7VA5kPQI9AABwGleuXNHUqVPl5uamIUOG6PHHH9emTZu0ZMkSubm5pbjPoUOHVK1aNX377bcaPXq0tm/fLm9v7xyuHMg+fL4EAAAcXmxsrD766CONGTNGR48eVdOmTTV69GjVr18/1X2OHTumsmXLqlKlSvL391eHDh1SDf2AM2OGHgAAOKzr169r/vz5qlq1qnr06KFKlSpp1apVCg8PTzXMnz59Wt26dZObm5v2798vSXrnnXcI88i1CPQAAMDhxMXF6YsvvlD16tXVuXNnlSpVSsuXL9dvv/2mJk2apLhPfHy8PvroI1WtWlXz589Xz549b1l3HsiNaLkBAAAOwxijpUuXasSIEdq5c6dq1qypZcuWqWXLlrIsK9X9rl+/rn79+umPP/5QgwYNNHv2bNWsWTMHKwfshxl6AABgd8YY/fDDD/L09NTLL7+s69eva+HChdq6datatWqVapiPjY2VJOXPn1916tTRRx99pF9//ZUwjzyFQA8AAOzGGKOVK1fqiSee0PPPP6+zZ89qwYIF2rFjh9q3b698+VKOKsYYff311/q///s//fbbb5Kkjh07qlOnTqnuA+RW/B8PAADsYs2aNXr66afl5eWlI0eOKCQkRHv37tWbb7552xs9HThwQM8++6zatGmjkiVLqkiRIjlYNeB4CPQAACBHbdy4US1atNCTTz6p3bt3a/r06dq3b5+6deumAgUK3HbfSZMm6eGHH9bvv/+uqVOnatOmTapbt24OVQ44Ji6KBQAAOWLbtm0aMWKEvv32W917772aMGGC/Pz8MjTD7uLiohdffFGTJ0/WAw88kI3VAs6DGXoAAJCtdu/erfbt26t27dr65ZdfNHr0aB08eFADBw5MM8z/888/euWVV/TZZ59Jkt5++20tXLiQMA/cgEAPAACyxYEDB/Tmm2+qRo0a+vHHHzV8+HAdPHhQw4cPV7FixW677/Xr1zV16lS5u7tr2bJlOnPmjCTddulKIK+i5QYAAGSpw4cPa8yYMfroo49UoEAB9e/fX4MHD9Z9992Xrv03btwoHx8fbd26VS1atNDMmTNVpUqVbK4acF4EegAAkCWOHTumcePGKSQkRJLk5+enoUOHqmzZshk6ztGjR3Xy5EmFhobqpZdeYlYeSAOBHgAA3JGTJ08qKChIM2fO1PXr19W5c2e98847qlixYrr2j4+P1yeffKJz587p7bffVqtWrdSsWTOWowTSiR56AACQKWfPntXw4cPl5uamyZMnq127dtqzZ4+Cg4PTHea3b9+uRo0aqVOnTvruu+9kjJFlWYR5IAMI9AAAIEPOnz+v0aNHq3Llyho7dqyee+457dy5UwsWLEh3r/vFixc1aNAg1alTR7t379ZHH32kFStW0F4DZAKBHgCAPKxM+YqyLCvNrzLlK+ry5cuaMGGC3NzcNGLECDVu3Fjbtm3TwoUL5e7unqHn/fPPPzVlyhR17txZe/fuVadOnZQvH7EEyAx66AEAyMNOHP1b6r8y7XGTm8rNzU0nTpxQixYtNGrUKD366KMZeq4DBw7oxx9/VO/eveXh4aEDBw7owQcfzGzpAGwI9AAAIF3c3d0VGhqqhg0bZmi/q1evKigoSOPGjVOBAgXUrl07lS5dmjAPZBECPQAASJfVq1dnuMc9PDxcfn5+2rdvn9q2baspU6aodOnS2VQhkDcR6AEAQLpkNMyfPXtWrVu3VpkyZfTTTz+pefPm2VQZkLdx9QkAAMgy169f11dffSVjjEqUKKEVK1Zo+/bthHkgGxHoAQBAlli/fr0ee+wxtW/fXuHh4ZKk+vXr66677rJzZUDuRqAHACCP2r17d5Yc58yZM+revbvq16+vf//9V4sXL5aXl1eWHBtA2gj0AADkMbGxsRo9erRq1659x8cyxqhZs2aaN2+e+vXrp927d6tNmzbcIArIQVwUCwBAHrJx40Z16dJF27dvV/v27bVo0aJMHWfXrl2qUqWKChUqpIkTJ6pkyZKqVatWFlcLID0I9AAA5AGXL1/WiBEjNGXKFJUpU0bffPONWrZsqYg1a3VictM09y9droIk6eLFixo5cqSmTJmisWPHyt/fX40bN87m6gHcDoEeAIBcbuXKlfLx8VFUVJR8fX0VGBioe+65R5J0/MjhdB3DGKMlS5bo7bff1pEjR9SlSxd16dIlO8sGkE700AMAkEudPXtWXbp0kZeXl/Lly6eIiAjNmTMnKcxnhL+/v15++WWVLFlSv//+u+bNm6dSpUplQ9UAMooZegAAcqElS5aoZ8+eOnnypPz9/fXee++pcOHCGTrG1atXdfXqVRUvXlzt2rVT2bJl1bt3b+XPT3wAHAkz9AAA5CLHjh3Tyy+/rJdfflllypTRhg0bNH78+AyH+ZUrV6pWrVoaMGCAJMnT01P9+vUjzAMOiEAPAEAuYIzRRx99pOrVq+uHH37QuHHjtGHDBnl4eGToOMeOHdOrr74qLy8vxcXF6eWXX86migFkFd5mAwDg5KKiouTj46OVK1fqySef1Ny5c1W1atUMHcMYo+XLl6tDhw66cuWK3nvvPfn7+2d4Zh9AzmOGHgAAJxUXF6fJkyerRo0a2rBhg2bPnq2IiIgMh/nr169LkmrWrKmnn35aO3bs0Pvvv0+YB5wEM/QAADih7du3q0uXLtq4caOef/55zZ49W+XLl8/QMc6ePauhQ4fqwIEDWrFihSpUqKClS5dmU8UAsgsz9AAAOJGrV6/q3XfflYeHhw4dOqQvv/xS3377bYbCvDFGCxYsUNWqVTVv3jzVrFlT165dy8aqAWQnZugBAHASa9euVdeuXbV79269/vrrmjJlSobXgv/777/1+uuv69dff1X9+vU1e/ZsPfLII9lUMYCcwAw9AAAO7uLFi+rTp48aNmyoS5cu6ccff9Snn36aqRs7ubq66ty5c5o7d67WrFlDmAdyAQI9AAAO7KefftLDDz+smTNnqmfPntqxY4eeeeaZdO9vjNGyZcvUvHlzxcbGqlixYtqyZYu6du2qfPmIAUBuwJkMAIADOn36tN58800988wzKlKkiNasWaMZM2aoWLFi6T7GwYMH9cILL6h169b6559/9M8//0iSLMvKrrIB2AGBHgAAB2KM0cKFC1WtWjV9+eWXevfdd7V161Y98cQT6T7GtWvXNHbsWFWvXl0RERGaOHGiIiMjValSpewrHIDdcFEsAAAO4siRI+rRo4e+//57PfroowoPD1etWrUyfJx8+fJpyZIlev755zVlypQML2cJwLkwQ29jWZabZVnBlmW1sXctAIC8JT4+XnPmzFH16tW1cuVKTZo0Sf/9738zFOaPHz+u7t276/Tp03JxcVFERIQWL15MmAfyAAK9JMuyvCS52b5K2rkcAEAe8ueff6pJkybq0aOHHn30UW3fvl39+/eXi4tLuvaPi4vTzJkzVbVqVc2fP1+///67JGWo1x6AcyPQSzLGhBtjwiVF27sWAEDecO3aNY0fP161atXStm3b9OGHHyo8PFxVqlRJ9zE2btyoxx57TL1791a9evW0Y8cOtWzZMhurBuCI6KEHACCHRUZGqkuXLtq6dateeuklzZw5U2XLls3wccaPH69jx45p0aJFatu2LavXAHmUwwV6Ww97e2NM21S2u0oaKumA7aEqxhj/HCoPAIBMi4mJ0ciRIzVx4kTdd999+vrrr/XSSy+le39jjD799FM9/vjj+s9//qPZs2frrrvuUvHixbOxagCOzmECvWVZwbZv0+pjXyzJ1xgTZdvPzbKsMGOMd3bXCABAZv3yyy/q1q2b9u3bp86dO2vixIkqUaJEuvffuXOn/Pz89Ouvv6pv376aMmWK7r///mysGICzcJgeemOMrzHGVwmBPUW22fuoxDBv2y/qhm0AADiUc+fOqXv37mrcuLHi4uIUHh6uDz/8MN1h/tKlSxoyZIhq166t7du3KyQkRJMmTcrmqgE4E4eZoU+n9pLCUng8TJKvpFBJsizLR1LdNI4VZowJzdryAAD4n2+//VZ+fn46duyY+vfvr1GjRqlo0aIZOsaECRMUGBioTp06KTAwUPfdd182VQvAWTlboPeSFJzC41GSPBN/MMaE5FhFAADc5N9//1WfPn20aNEi1axZU0uWLNFjjz2W7v0PHTqks2fPqk6dOhowYIC8vLzUsGHDbKwYgDNzmJabtNguhnWVdCaFzdG2bQAA2I0xRp988omqVaumpUuXatSoUdq0aVO6w3xsbKwCAgJUvXp1+fr6yhijYsWKEeYB3JYzzdCnecMny7JcjTHRGT2wZVkeSpj995LkZnvzEJKZYwEA8qa//vpLvr6++vnnn/XEE09o7ty5ql69err3X716tfz8/LRnzx69/PLLmjJlCstQAkgXyxhj7xqSsfW/+xpj6t70uJsSlqqsa4yJvGmblxL66EvkRAi31egjSaVLl667cOHCNPe5ePGi7r777uwuDUAW4rxFesTFxWnZsmWaN2+eJMnHx0etWrVSvnzp/xB8w4YN8vf31wMPPKA+ffqoXr162VVursY5i9ysSZMmm40xniltc6YZeodh69EPkSRPT0/TuHHjNPeJiIhQesYBcByct0jLrl271KVLF61bt04tWrTQnDlz9OCDD6Zr37i4OO3bt0/u7u568sknVaRIEXXp0kWFCxfO5qpzL85Z5FVO00OvlHvnk6FFBgCQE2JjYzVq1CjVqVNHf/75pz799FP9+OOP6Q7zGzdu1GOPPaYnn3xS586dk4uLi3r16kWYB5ApThPobWE9Wgk3nrqZm20bAADZasOGDapbt67ee+89vfTSS9q9e7def/31dPW7nz17Vn5+fqpXr57++ecfzZgxg7u8ArhjztZys0kpXxxbRVJ4DtcCAMhDLl26pHfffVfTpk1T2bJl9e233+qFF15I9/7Hjh1T7dq1derUKfXu3VujRo3SPffck40VA8grnGaG3maxJO8UHveStCiHawEA5BHh4eGqWbOmpkyZIh8fH+3cuTPdYT46OlqSVLZsWXXr1k2bNm3StGnTCPMAsowjBnpXpbKmvO1iVDfbijeSkpacPMNdXwEAWe3s2bPq3LmzvL29lT9/fv3yyy+aPXt2usL45cuXNXToUFWoUEH79u2TJI0ZM0Z16tTJ7rIB5DEO03JjWVagEoJ8O0mulmUtVsKFsME3LVPZVNJQy7IO2H6uYoxJadYeAIBM+/rrr9WzZ0+dOnVKQ4YM0YgRI9J90eq3336rPn366K+//lLHjh3l6uqavcUCyNMcJtAbY/xt3/qmMS5akv/txgAAkFnHjh1Tz549tXTpUtWpU0fLly9P96x6fHy82rRpo6VLl6pGjRr69ddf9eSTT2ZzxQDyOkdsuQEAIMuUKV9RlmWl+VWmfEXNmzdP1apV0/LlyzV+/Hht2LAhXWE+Li5OkpQvXz5Vr15dEyZMUGRkJGEeQI5wmBl6AACyw4mjf0v9V6Y9bnJTdevWTY0aNdLcuXP10EMPpev4q1evVs+ePTV79mw1atRIY8aMudOSASBDmKEHAMBmzpw5WrVqVbrC/IkTJ/TGG2/o6aef1pUrV2SMyYEKAeBWBHoAAGx8fX2VL1/a/zR+9NFHqlq1qhYtWqThw4drx44daty4cfYXCAApoOUGAIAMunDhgjw9PTVr1ixVrVrV3uUAyOOYoQcAIA3R0dHq1auXPvvsM0lS7969FRYWRpgH4BCYoQcAIBXGGH3++ecaMGCATp06pfvuu0+S0tWWAwA5hUAPAEAK9uzZIz8/P61evVr16tXTTz/9xF1eATgkphgAALnSlStXNHTo0Ezv/+eff2rr1q0KDg7W2rVrCfMAHBYz9ACAXGfDhg3q1KmTdu3alaH9vvvuOx0+fFg9e/ZUy5YtFRUVJVdX1+wpEgCyCIEeAJBrXLlyRe+//74mTJigBx54QMuXL1fHrj46MblpmvsWuusutWzZUh4eHurevbtcXFwI8wCcAi03AIBcYf369fLw8FBgYKA6d+6sHTt2qEWLFjp+5LCMMSl+Xb16VQEBASpcuLBc8uVTYGCg1q1bJxcXF3u/HABINwI9AMCpXblyRf7+/nriiSd04cIF/fTTT5o7d67uueeeNPfds2eP3nnnHbVo0UK7d+/W4MGDVaBAgRyoGgCyDi03AACntX79enXs2FF79uxR165dNXHixDSD/L///qvvvvtOXbp0Ua1atbR9+3ZVr149hyoGgKzHDD0AwOncOCt/6dKldM3Kx8XFafbs2apatar8/Pz0999/SxJhHoDTI9ADAJzKunXrVKdOHQUFBalz587avn27mjdvftt9Nm/erPr168vPz08eHh7atm2bKlSokEMVA0D2ouUGAOAUrly5ohEjRmjSpEkqV66cfv75ZzVr1izN/S5evKimTZuqcOHC+uKLL/TKK6/IsqwcqBgAcgaBHgDg8NatW6dOnTppz5496tatmyZOnKjixYunOt4Yo59//lnNmzfX3XffraVLl8rDwyNdF8oCgLOh5QYA4LBiYmI0aNAgNWjQQJcuXdLPP/+skJCQ24b5PXv2yMvLS88884y+/fZbSVKTJk0I8wByLQI9AMAh/fe//1WdOnU0ceJEde3aVTt27Lhti83ly5f1zjvvqFatWoqMjNTs2bP1/PPP52DFAGAftNwAABxKTEyMRowYocmTJ6t8+fJasWKFvL2909zv2Wef1S+//KI333xTQUFBKl26dA5UCwD2xww9AMBh3Dwrv3379tuG+cOHD+vq1auSpOHDhysiIkILFiwgzAPIUwj0AAC7i4mJ0cCBA9WgQQPFxMQoLCxMwcHBqfbKx8bGKjAwUNWqVdOECRMkSV5eXmrUqFFOlg0ADoGWGwCAXa1du1adOnXSn3/+KV9fXwUFBd32otdffvlFfn5+2rVrl1q1aqU33ngjB6sFAMfDDD0AwC4SZ+UbNmyoK1euKCwsTHPmzLltmA8ICFDjxo116dIlffvtt1q2bJkefPDBHKwaABwPM/QAgBx346x89+7dFRQUpGLFiqU4Nj4+XjExMSpatKiaN2+uCxcuaPjw4SpSpEgOVw0AjokZegBAjrl8+bIGDBighg0b6urVqwoPD9fs2bNTDfORkZGqX7+++vTpI0ny8PDQuHHjCPMAcAMCPQAgR/z++++qXbu2Jk+eLF9fX23fvl1NmzZNcey5c+fUp08fPfroo/rrr79SHQcAINADALLZ5cuX1b9/fz355JOKjY1Nc1b+t99+k7u7u2bOnCk/Pz/t2bNHr776ag5XDQDOgx56AEC2+f3339WpUyft27dPPXr0UGBg4G175fPly6fKlSurWrVq+v7771W3bt0crhgAnE+mAr1lWbUleUmqIqmkJFdJUZKiJZ2WFGmMWZUlFQIAnM7ly5c1fPhwTZ06VRUrVtTKlSv19NNPpzp23Lhx2rx5s3788UeVL19eq1bxTwgApFe6W24sy3rasqyvLMvaJ8lXkiUpXFKIpCBJoZI22R5vZ1nWJsuyFlmWlfLf4ACAXGnNmjWqXbu2pkyZou7du2v79u2phvkffvhBDz/8sMaOHav77rtPMTExOVwtADi/NGfoLcuqLClQ0gFJ/saYg+k9uGVZ90jysSzL17bvocwWCgBwbJcvX9Y777yjadOm6cEHH7ztrPy///4rX19fLVu2TNWqVdPq1avVuHHjnC0YAHKJ2wZ6y7KaKqG1ppsx5lxGD27bZ4It2A+1LGsFrTgAkPusWbNGnTp10v79++Xn56fAwEDdfffdqY4vUqSIdu/erYCAAPXv318FCxbMwWoBIHdJteXGNjPvaowZmpkwfyNjzDljzJCEw1qV7uRYAADHcfnyZfXr109PPfWUrl+/rlWrVmnWrFkphvnffvtNL730kq5evaq7775bO3bs0JAhQwjzAHCHUg30xpiDxpivs/LJjDErabsBgNzht99+0yOPPKKpU6eqR48e2r59u5o0aXLLuJMnT6pjx4566qmnFBkZqYMHEzo38+dnoTUAyAqsQw8AyJDLly+rb9++atSokeLi4lKdlY+Pj1dwcLCqVq2qL774QkOHDtWuXbvk7u5up8oBIHfK7LKVTZVwoexpY0xz22MvKWEJy3Bm4QEgd/rtt9/UqVMnHThwQD179tT48eNT7ZU3xmjevHmqXbu2Zs2apWrVquVwtQCQN2R2ht7LGOOphCUrZVnWHEmLJbWTFGpZ1sAsqg8A4AAuXbqkt99+W40aNVJ8fLxWr16tmTNn3hLmz507p8GDB+vUqVNycXHRTz/9pJUrVxLmASAbZTbQn5GkG3rsfSSFGGOa2YL+OduMPQDAyf3666965JFHNH36dPXs2VN//PHHLUtMGmO0cOFCubu7a+LEifr5558lSffee68sy7JD1QCQd2Q20JvEb2ztN0YJM/QJG42ZK8ntzkoDANhTSrPyM2bMuGVW/s8//1SzZs3UoUMHlStXThs2bNBrr71mp6oBIO/JbKC/cQbeX5JSWF/+jpa6BADYz42z8r169UpxVj7RiBEjtHHjRs2aNUvr16+Xp6dnzhYLAHlcpgK9bQb+Mcuy4pVw46khKQ27k8IAADnv0qVL6tOnjxo1aiRjjCIiIlKclf/xxx+1d+9eSdKUKVO0Z88e+fn5ycXFxR5lA0CelullK40xQ4wx+WxfE6SEi2Mty1poG1IiSyoEAOSIX375RbVq1dKMGTPUu3dv/fHHH2rUqFGyMX///bdeeuklPffcc5o4caIkqWzZsipTpow9SgYAKOvXoQ+X9H+WZQXItgIOAMCxXbp0Sb17905qqYmIiND06dNVtGjRpDHXrl3TxIkTVa1aNf30008KCAjQrFmz7FQxAOBGWXqbPmNMqKTQrDwmACD7REREqEuXLoqKilKfPn00bty4ZEE+0eTJkzVkyBC1bNlS06ZNU6VKlXK+WABAirjvNgDkQRcvXtTQoUM1c+ZMubm56ZdfftFTTz2VbMzJkyd1/Phx1axZU35+fnr44Yf1/PPP26liAEBqUm25sSyrsmVZxbPyySzLKm5ZVqWsPCYAIGMiIiJUq1YtzZw5U3369NEff/yRLMzHx8crJCREVatW1WuvvSZjjIoVK0aYBwAHlWqgN8YclORrWVbtrHgiy7LqSPIxxhzKiuMBAP6nTPmKsiwrza8ixV3VpEkT5cuXT7/88oumTZuWrMVmy5YteuKJJ+Tr66tatWpp4cKF3BgKABzcbVtujDETLMvqZlmWr6TAzIRx24z8EEmbjDETM1UlAOC2Thz9W+q/Ms1xMZObptor/9tvv6lx48YqVaqUPvnkE73++uuEeQBwAmn20Btj5lqWdY+kQMuyKksKkxSphIB+/ubxtjYdT0l1JXlLOiDJ3xjDjaYAwAFMmzYt6XtjjA4fPqwHH3xQTzzxhEaPHq0ePXqoRAlWHgYAZ5GuZSuNMeeMMd2NMc0lHZTUTtIqy7L2W5Z12rKsM5ZlxVmWdVrSZkltJUUZY5oZY3oQ5gHA8fz5559q3ry56tatq9OnT8vFxUXDhg0jzAOAk8nwKjfGmK8lfZ0NtQAAckBMTIwCAgIUGBiou+66S+PGjZOrq6u9ywIAZBLLVgJAHlOjRg1FRUXptdde08SJE7nLKwA4OQI9AOQF165IBe6SJLVu3VrPPvusnn76aTsXBQDICgR6AHBi0dHRGjFiROoD4q5LW5ZIG76QOsyQJE2cyIJjAJCbEOgBwAkZY/Tpp59q0KBBOnnyZMqDjm6XVk6TTh2U3B6XXArmbJEAgByRrlVuAACOY/v27Xrqqaf01ltvqXLlytq4cWPyAcZIYZOlRX2lq5eklqOkF8dKxUvbpV4AQPZihh4AnMT58+f1/vvva/r06XJ1ddXcuXPVuXNn5cuXT6XLVdCJyU1v3enCv9K3/2vJKV2uQg5WDADICczQA4CDM8boiy++UNWqVTV16lR16dJFe/fuVdeuXZUvX8Jf4z//8J3q16+v1atXyxiT6tfxI4ft/GoAAFmNQA8ADmzXrl1q2rSpXnvtNZUrV07r1q1TcHCw7r33XkkJs/b9+vWTh4eH9u/fr/Pnb7mBNwAglyPQA4ADunjxogYPHqxHHnlEW7du1ezZs7V+/Xo99thjSWOWLl2qatWqadq0afLx8dGePXvUsmVLO1YNALCHO+qhtyyruCQvSW7GmIm2x+6R1NQYsyQL6gOAPMUYo9DQUPXr109Hjx5V586dNX78eN133323jD1y5IhKly6tJUuWqF69enaoFgDgCDI9Q29Z1mxJZyWFSgq8aXOoZVld7qQwAMhr9u7dq2bNmqldu3a67777tHbtWn344YdJYT4mJkbvvfeePvvsM0mSn5+fNmzYQJgHgDwuU4HesqzxknwldZfU7MZtxphzkubatgEA0nDp0iUNGzZMNWvW1MaNGzVjxgxt3LhR9evXTxrz008/qWbNmho1apQ2bNggSXJxcVH+/CxWBgB5XWZn6NtIGmyMmSvpQArbN0nyyHRVAJAHGGO0dOlSVa9eXQEBAerQoYP27t2rXr16JQX1o0ePqm3btnrmmWeUP39+hYeHa/r06XauHADgSDIb6Esqod0mNa6SojJ5bADI9fbv36/nnntOL730kooXL65ff/1VCxYsUOnSyW/+tGXLFn3//fcaM2aMtm3bpqZNU1hrHgCQp2X2s9qVkoZYlvXVzRtsF8UOlbToTgoDgNwoJiZGAQEBCgwMVKFChTRlyhT17NlTBQoUSBrz+++/a9euXerWrZuef/55HTx4UGXKlLFj1QAAR5bZGfqukkpJipY0XpIsy2piWdZAJczMG0n+WVEgAOQW3333napXr67Ro0erTZs22rNnj/r27ZsU5k+dOqUuXbqoYcOGmjBhgmJjYyWJMA8AuK1MBXrbha91Ja2S1FaSpYRZ+yBJmyV5GmO4uwkASDp48KBatmypli1bqnDhwlq1apU+//xzPfDAA5Kk+Ph4zZs3T1WrVtUnn3yiQYMGKTIyUgULFrRz5QAAZ5Dp5RGMMVGSvG0tNp62hzfZwj4A5HlXrlxRUFCQAgIC5OLiogkTJqhPnz63BPW9e/fK19dXDRo00AcffKAaNWrYqWIAgDO64/XObAF+ZRbUAgC5xvLly9W7d28dOHBA7dq106RJk1S+fPmk7RcuXND333+vDh06qFq1alq3bp08PT1lWZYdqwYAOKNMB3rLsmor4S6x96YyxBhjhmX2+ADgjP766y/17dtXy5Yt03/+8x+tWLFC3t7eSdsT7wTbt29fHT9+XJ6ennrooYf06KOP2rFqAIAzy1Sgtyyrm6Q5SuidT42RRKAHkCdcvXpVkyZN0pgxYyRJ48aNU//+/VWoUKGkMfv371evXr30888/q06dOlqyZIkeeughe5UMAMglMjtD7y/pnBIuiN2UdeUAgPMJCwtTr1699Oeff+qll17SlClTVLFixWRjrly5oieeeEJXr17V9OnT1aNHD+7yCgDIEpn918RNCXeKpXceQJ515MgR9e/fX4sXL1aVKlW0fPlytWjRItmYtWvXqn79+rrrrrv0ySef6JFHHlHZsmXtVDEAIDfK7Dr04VlahZ1ZluVqWZaP7SvYsqw29q4JgOOKjY3VhAkT5O7uru+++06jRo3Sjh07koX5o0ePqm3btmrQoIEWL14sSWrRogVhHgCQ5TI7Q99d0grLsg4YY5ZmZUF2EmiM8U38wbIsY1lWXWNMpD2LAuB4Vq9erZ49e2r37t1q2bKlpk6dqsqVKydtv379umbMmKERI0bo+vXrGj16tFq1amXHigEAuV2mAr0xJsqyrEhJobYl1qJTHmZK3UFtOcnTsiwvY0ziJw/RSmgrItADkCT9888/GjhwoL788ktVrlxZ3333nZ5//vlbxr388sv69ttv9eyzz2rGjBlyc3OzQ7UAgLwks6vczJHURgnBN0rSmSysKccZY+omfm9ZlqskVxHmAUi6du2aZs6cqffee0+xsbEaMWKEhgwZosKFCyeNOX36tIoWLaq77rpLPXv21FtvvaXWrVuzpjwAIEdktuWmnaQwY0zzrCxGkmz96+2NMW1T2e4qaaikA7aHqhhj/LOwhEBJbW13wgWQh23btk29e/fWjh079Mwzz2j69On6v//7v6Tt8fHx+vjjjzV48GD16tVL77//vpo1a2bHigEAedGdrJm2OMuqkGRZVrDtWzdJJdN4Xt/EwG1ZlptlWWHGGO/b7JOe53dVwhuV6Ds5DgDnd+LECQ0aNEiffvqpKlasqKVLl6pVq1bJZtz/+OMP9ejRQ2vXrlXDhg3Vpg3X0gMA7COzq9wsllQ3zVEZYIzxtV2YmuobBdvsfdSNs+c3BPs7+tfUGBNtjAmxzfbPZaUbIO9JvKD1P//5jxYuXKjXXntNu3bt0osvvpgszH/wwQfy8PDQn3/+qfnz5+vXX39VjRo17Fg5ACAvy+wM/RxJKy3LKiHpKyXMat/SR2+M2ZrpylLWXlJYCo+HSfKVFCpJlmX5KO03HGHGmNDEmXljTMgN28JvPB6A3G/t2rXq2bOntm7dKm9vb82cOVP//POPihYtKkkyxig2NlaFChVS/fr11bVrV40bN04lS97uA0UAALJfZgP9SiVcONpOCXeLvZklyUhyyeTxU+MlKTiFx6MkeSb+cFM4T4unEvrmM7IPgFzi33//1ZAhQzR//nyVL19eixcv1ssvvyzLsvTPP/9Ikvbv369evXqpbNmymj9/vurUqaM5c+bYuXIAABJkNtB30+373LPcDavPpLSiTrRtW2ZskhRw02NeSvmNCoBcIi4uTiEhIRo2bJguXryowYMH691339Xdd9+dNCY2NlYjR45UQECAChYsqLFjx8oYw+o1AACHktl16L/O6kLSIc03EJZluRpjojNyUGNMtGVZ4ZZlDVbCG4O6krrdsCY9gFxmw4YN8vPz0+bNm9WkSRPNmjVL1apVSzYmMjJSnTt31tGjR/XKK69o0qRJeuCBB+xUMQAAqbuTVW5yDdsdYVl3HsjlTp8+raFDh2revHkqU6aMvvzyS7Vv3z7ZjHviDHzp0qVVtGhRrVixQt7ed7SIFgAA2SrNQG9Z1tOSZIxZdcNjtdNz8Gy4KNYh2C669ZGk0qVLKyIiIs19Ll68mK5xALJefHy8fvzxR82dO1cXL15UmzZt9NZbb6lo0aL65ZdfJCW04CxdulRbtmzRmDFjZFmWJk6cqAIFCnDuAk6Cf2uRV6Vnhj5cCTdxeuiGxyKVcNFrarLjotg070ab0XabzLJddBsiSZ6enqZx48Zp7hMREaH0jAOQtTZv3iw/Pz9t2LBBTz31lGbNmnXLEpPr1q1Tjx49tHXrVj3zzDOqW7euihcvznkLOBnOWeRV6Qn0QySdvemx7rp9oM9ytl73aCXceOrm9hg3cUMoADc4e/as3nnnHc2ZM0f333+/Pv30U7322mvJ2mvOnTunQYMGae7cuSpfvry+/vprtW7dmoteAQBOJc1Ab4wJSuExey3xuEkpXxxbRQmfJADIxcqUr6gTR/9Oc1zxkqVUMJ905swZ9e7dWyNHjpSrq+st4/Lnz69Vq1Zp4MCBeu+995KtcAMAgLPI1EWxlmUNlBRijDmfyvbZkhbf2HefRRZL8tata8Z76dalJwHkMieO/i31X5nmuPOTm+qJJ57QrFmzVLt27WTb/vjjD02YMEFz585V0aJFtWPHDt11113ZVDEAANkvXyb3C9QNN3JKgSXJP5PHdlUqa8rbPhlwsyzLLemJLMtD0hljDHd1BZDkt99+SxbmL1y4oAEDBsjDw0PLly/X7t27JYkwDwBwepldtjKtBtMDyuCNmSzLCtT/7j7ralnWYiVcCBtsW1YyUVNJQy3LOmD7uYoxhjXlACSTL1/CfIUxRl9//bX69u2ro0ePqlu3bgoICNC9995r5woBAMga6Q70lmU1ldTmhof8Lcu6ObSXVEIo91IGL1I1xiTO6PumMS5amZ/9B5DHGGM0YcIElSpVSosXL1b9+vXtXRIAAFkqIzP0bpLa2743SuhlT0m0pJVKaMsBALsICAhQ165ddd9992nZsmW67777lD8/99IDAOQ+6e6hN8bMNcaUNMaUVELLjZcxJl8KXyWNMc2MMWlfuQYAWe2vTZKkYcOGacmSJZKksmXLEuYBALlWZv+FC5UUlZWFAEBq/v77b7399tu3H3TxlPTLbGlvhCRpxYoV8vbm8hoAQO6XqVVujDHtjDGHsrgWAEjm+vXrmjx5sqpVq6affvrp9oPXzJP2/y7Vf0uSCPMAgDwjs8tWAkC2WrdunTw9PTVgwAA1atRIO3fuvHXQP7uks7YbTTXsKr31kVT/zZwtFAAAO6OpFIBDOXv2rIYOHaqQkBA98MAD+vrrr9W6dWtZlqXS5SroxOSmaR6jdLkKOVApAACOgRl6AA7BGKPPPvtMVatW1bx589SvXz/t3r1bL730kiwr4dYX/xw+pI8++kilSpWSi4uLBgwYoPPnz8sYk+zr+JHDdn41AADkHGboAdjdnj175Ofnp9WrV6tevXpasWJFsru8Jvrggw/Uu3dvPfHEE5o9e7Zq1aqV88UCAOBgCPQA7CYmJkbjxo1TYGCgihYtqjlz5qhbt25Jd3mVpIsXL+rIkSNyd3dXx44ddc899+i1115LNgYAgLyMfxEB2MXPP/+sGjVqaMyYMWrfvr327NkjX1/fpKBujNGSJUtUrVo1tW7dWnFxcbr77rv1xhtvEOYBALgB/yoCyFH//POP2rdvrxYtWqhAgQJauXKlPv30U5UuXTppTFRUlJ5//nm9/PLLKlmypD788EO5uLjYsWoAABwXLTcAckRcXJxmzZql4cOHKzY2VqNHj9agQYNUqFChZOM2b96shg0bKn/+/Jo8ebJ69+7NXV4BALgN/pUEkO02btyo7t27KzIyUs2aNdOsWbP0f//3f8nGnDhxQqVLl1bt2rXVr18/9ezZU+XKlbNTxQAAOA9abgBkm3PnzqlXr16qV6+ejh07pkWLFumnn35KFub/+ecfdejQQdWrV9epU6fk4uKicePGEeYBAEgnAj2ALGeM0cKFC+Xu7q7Zs2erV69e2r17t9q1a5e0pvz169c1ffp0ubu7a+nSperTp4/uvvtuO1cOAIDzoeUGQJbat2+fevbsqbCwMNWtW1ffffedPD09k425cOGCGjVqpC1btqh58+aaOXPmLS04AAAgfZihB5Alrl69qlGjRqlmzZpav369Zs6cqfXr1ycL89euXZMkFStWTE888YS++uorLV++nDAPAMAdINADuGMrV65UrVq19N5776l169bas2ePevbsmbTUpDFGCxYsUOXKlbV3715J0syZM9W2bdukFhwAAJA5BHoAmXb8+HG99tpr8vLyUlxcnH7++Wd9+eWXKlu2bNKYHTt2qFGjRurYsaMqVqwoY4wdKwYAIPch0APIsLi4OM2ePVvu7u4KDQ3Ve++9px07dqhZs2bJxr3zzjuqU6eOdu7cqXnz5mnNmjVyd3e3U9UAAOROXBQLIEO2bNmi7t27a8OGDWratKk++OAD/ec//0lxbGxsrN566y2NHz9epUqVyuFKAQDIG5ihB5AuFy5cUL9+/eTp6alDhw7p888/V1hYWLIwf/DgQb3wwgtatWqVJCkoKEjz5s0jzAMAkI0I9ABuyxij0NBQubu7a9q0afL19dWePXv06quvJl3QevXqVY0dO1bVq1dXRESEjh49Kklc8AoAQA6g5QZAqqKiotSrVy8tX75ctWvX1pIlS1SvXr1kYyIiItS9e3ft3btXbdu21eTJk1W+fHk7VQwAQN5DoAdwi9jYWE2cOFGjR49W/vz5NWXKFPXq1Uv589/6V8bOnTt1/fp1LV++XC1atLBDtQAA5G0EegDJ/PLLL+rRo4d2796tl19+WVOnTk024379+nXNnj1b9957r1599VV1795dnTt3VuHChe1YNQAAeRc99AAkSSdPnlTHjh3VuHFjxcTE6IcfflBoaGiyML9hwwY99thj6tOnj77//ntJkouLC2EeAAA7ItADeVx8fLzmzZunqlWr6osvvtCwYcO0c+dOPfvss0ljzp49q+7du+vxxx/X8ePHtWjRIn3++ed2rBoAACSi5QbIw/744w91795d//3vf/XUU09p9uzZql69+i3j/vvf/2ru3Ll6++23NXLkSBUvXtwO1QIAgJQwQw/kQRcvXtSgQYPk4eGhffv26eOPP1ZERESyML9z504tWLBAkvTss89q3759mjJlCmEeAAAHQ6AH8phvvvlG1atX18SJE9W5c2ft2bNHb731VtKa8ZcuXdKQIUNUu3ZtDR06VDExMZIkNzc3e5YNAABSQaAH8oi//vpLrVq10osvvihXV1etWbNGISEhuvfeeyUl3EBq2bJlql69ugIDA/Xmm2/qjz/+4IJXAAAcHD30QC537do1TZkyRSNHjpQkTZgwQW+//bYKFCiQbNzBgwfVpk0bVa9eXWvWrFGDBg3sUS4AAMggAj2Qi61Zs0Y9evTQjh071KpVK02fPl0VK1ZM2n716lX99NNPatWqldzc3BQeHq4GDRrcEvYBAIDjouUGyIVOnz6trl276sknn9S5c+e0bNkyLVu2LFmYX7VqlR555BG9+OKL2r59uySpcePGhHkAAJwMgR7IRYwx+vjjj1W1alV9/PHHGjRokHbt2qVWrVoljTl+/Lhee+01NW3aVNeuXdOPP/6omjVr2rFqAABwJ2i5AXKJXbt2qUePHvr111/VoEEDzZ49+5agfu3aNdWrV0/Hjx/XiBEjNGTIEC56BQDAyRHoASd3+fJljRkzRhMmTFDx4sU1b948derUSfny/e8DuB07dujhhx9WgQIFNGPGDLm7u+s///mPHasGAABZhZYbwIn98MMPevjhhxUQEKDXX39de/bsUZcuXZLC/NmzZ+Xn56datWrp888/lyS1bNmSMA8AQC7CDD3gQMqUr6gTR/9Oc9x9ZcrpySfqacmSJapWrZoiIiLUqFGjpO3GGH322WcaOHCgTp06pbffflstW7bMztIBAICdEOgBB3Li6N9S/5Vpjjs5uamWL1+ugIAA9e/fXwULFky2vWPHjvrkk0/0+OOP6+eff1bt2rWzqWIAAGBvBHrASe3cuVOVK1dO+vnSpUtycXHRXXfdpQ4dOqhhw4bJ2m8AAEDuxL/0gJO6Mcx/8803ql69ugICAiRJLVq0ULdu3QjzAADkAfxrDzixQ4cOqWXLlnrxxRdVvHhxeXt727skAACQwwj0gJP6/PPPVb16da1atUoTJkxQZGSkGjZsaO+yAABADqOHHnAQV65cSd/A+DhJ0sMPP6znnntOkydPVoUKFbKxMgAA4MgI9ICdxcTEKCQkRIGBgbcfeOmM9MscydYXX7t2bS1evDgHKgQAAI6MlhvATi5fvqzJkyercuXK6tu3r6pWrZrywPg4aesy6eOO0r5fpWKlc7JMAADg4JihB3LYpUuXNHv2bE2YMEH//vuvnn76aX311Vd66qmnZFlW8sFnDkvLx0sn9koVPaSn+0glK0jrP7NP8QAAwOEQ6IEccvHiRc2aNUsTJ07UqVOn5O3trREjRiS7kLV0uQo6Mblpygc4HJkwS28bBwAAINFyA2S7CxcuKCAgQJUqVdKQIUPk6emptWvXasWKFbesSnPs77/06aefqk2bNoqPj5cxRnFxcTLGJPs6fuSwnV4NAABwNAR6IJucO3dOY8aMUaVKlTRs2DDVq1dP69at0/Lly1W/fv1bxu/atUtNmjTRG2+8ob///ltnzpyRJG4OBQAAboukAGSx6OhojRw5UpUqVdK7776rBg0aaOPGjfrhhx9Ur169W8bHxMRo6NCheuSRR/THH38oODhYa9eu1b333muH6gEAgLOhhx7IImfPntXUqVM1bdo0nTt3Tq1atdKIESPk4eFx2/2MMVq4cKFef/11BQYG6v7778+higEAQG5AoAfu0OnTpzVlyhRNnz5dFy5c0EsvvaR3331XtWvXTnWfQ4cOKSgoSJMnT1aRIkW0bds2FS9ePOeKBgAAuQYtN0AmnTp1SsOGDVOlSpU0duxYNW/eXNu2bdPXX3+dapiPjY3V+PHjVb16dS1YsECbN2+WJMI8AADINGbogQw6efKkJk6cqFmzZuny5ctq166dhg8frho1atx2v9WrV8vPz0979uxR69atNXXqVFWsWDGHqgYAALkVgR5IpxMnTmjixIn64IMPFBMTo1deeUXDhw9X9erV09zXGKOhQ4fq6tWr+v777/Xcc8/lQMUAACAvINADaTh+/LiCgoI0Z84cXb16Va+++qreeecdubu733a/uLg4zZ07Vy+99JLuv/9+ffXVVypVqpSKFCmSQ5UDAIC8gEAPpOKff/5RUFCQgoODde3aNb3++usaNmyY/vOf/6S576ZNm9SjRw9t2rRJ58+f1+DBg2mvAQAA2YJAD9zkyJEjCgwM1Ny5c3X9+nW9+eabGjZsmP7v//4vzX2jo6M1fPhwffDBBypdurS+/PJLtW/fPgeqBgAAeRWBHrD5+++/NX78eM2bN0/x8fHq2LGjhg4dKjc3t3Qfw9/fX/PmzVOvXr00evRo3XPPPdlYMQAAAIEe0F9//aWAgAB99NFHkqROnTpp6NChqlSpUrr23717t/Lnz6+HHnpI7733nnx9fdO8mRQAAEBWYR165FkHDx6Uj4+PHnroIc2fP19du3bV/v37FRwcnK4wf/nyZQ0bNkyPPPKIBg4cKEl64IEHCPMAACBHMUOPPCcqKkpjx47VJ598onz58snHx0f+/v6qUKFCuo/x3XffqXfv3vrrr7/05ptvasKECdlYMQAAQOoI9Mgz9u/fr7Fjx+rTTz9V/vz55efnp8GDB6tcuXIZOs7HH3+sTp06qXr16oqIiFCjRo2yqWIAAIC0EeiR6/35558aM2aMPv/8cxUsWFC9e/fW4MGDVbZs2XQfIzY2VkeOHJGbm5vatm2r8+fPq3v37ipYsGA2Vg4AAJA2Aj1yrd27d2vMmDFauHChChUqpH79+mngwIEqU6ZMho7zyy+/yM/PT9euXdPOnTtVtGhR9enTJ5uqBgAAyBguikWus2vXLnXo0EEPP/ywli1bpgEDBujQoUOaOHFihsL8v//+q7feekuNGzfW5cuXNXnyZBUoUCAbKwcAAMg4ZuiRa2zfvl2jR49WaGioihYtKn9/f/Xv31/33Xdfho+1a9cuNWjQQJcuXdI777yjYcOGqUiRItlQNQAAwJ0h0MPp/fHHHxo1apS+/vprFStWTEOHDlX//v117733ZvhY586d0z333CN3d3e9+eab8vPzU9WqVbOhagAAgKxByw2c1pYtW9S6dWs98sgjCgsL07vvvqtDhw5p7NixGQ7z0dHR6t27tx566CGdPHlS+fLl07Rp0wjzAADA4TFDD6ezefNmjRo1St9++63uuecevffee3r77bdVokSJDB/LGKMvvvhCAwYM0MmTJ+Xn58fKNQAAwKkQ6OE0NmzYoFGjRumHH35QiRIlNGrUKPXp00f33HNPpo4XExOj559/XqtWrZKnp6d++OEH1a1bN4urBgAAyF603CDHlSlfUZZlpflVpnxFSdK6dev0zDPPqF69evrvf/+rMWPG6NChQ3r33XczFebj4+MlSYULF5abm5s++OADrVu3jjAPAACcEjP0yHEnjv4t9V+Z9rjJTdW8eXOtWLFC9957rwICAtSzZ08VK1Ys08/9/fffa8CAAVq2bJmqVaumuXPnZvpYAAAAjoAZeji0LVu2KCgoSIcOHdKQIUMyHeYPHz6s1q1b64UXXpCLi4suXbqUxZUCAADYBzP0kizL8pFURdIiSSUl+Rpj2tq3KkjSwYMHVbRo0Ts6xpQpUzR8+HBJUmBgoPr27cuFrwAAINcg0CdwleQjabCkcEm+dq0GSe40zEvSsWPH1KxZM02dOlUPPvhgFlQFAADgOAj0CaKNMRlf8xBJLl26pJMnT+rff/9N+rr558THstu///6rwYMH6/XXX5eXl5cCAgLk4uKS7c8LAABgDwT6bFCmfMWECz/TULpcBR0/cjgHKsq42NjYZIE8rbB++fLlFI9TtGhR3X///br//vtVoUIFeXh46KOPPsqWmuPj4xUSEqKhQ4fq0qVLqlu3rry8vAjzAAAgV3O4QG9ZVhtJ7VPrYbcsy1XSUEkHbA9VMcb4Z9HzSpKbpEhjTHhmj5WRVVxySlxcnE6fPp3qrPnNP0dHR6d4nIIFC+r+++/Xfffdp/vvv19Vq1ZNCuyJj934c5EiRW45RnYE+i1btqh79+7asGGDmjRpolmzZqlatWpZ/jwAAACOxmECvWVZwbZv3ZRwYWpqFivhotUo235ulmWFGWO87+DpoyRtMsZE2455wLKsuok/OyJjjM6dO5dmME/8/tSpUzLG3HKcfPnyqVSpUkkhvG7durcE8xvDevHixWVZlh1e8e2tW7dOhw4d0meffaZXX33VIWsEAADIDg4T6I0xvlLSijMpXpRqm0WPSgzztv2ibDciamOMCc3kc988Gx+phE8B7njmPyMuXbqUrvaWxMeuXbuW4nFKlCiRFMKrVq2qJ598MtVZ9BIlSjhlS4oxRgsXLpRlWXrllVfk4+OjDh06yNXV1d6lAQAA5CiHCfTp1F5SWAqPhynhTUColPSmIK3bfoYlvgGwLMvDGBN5w7YzSvikINvVq1cvzT70u+++OymIV6hQ4baz6KVKlXL4JRlLl6uQrnaj0uUqpPj43r171bNnT61cuVLNmzfXK6+8IhcXF8I8AADIk5wt0HtJCk7h8ShJnok/GGNC0ntAy7I8JK2UdPMqN2cyU2BGlSxZMlkf+s2z6Kn1oTuzzF4IHBMTo3HjxikoKEiFCxfWBx98IB8fnyyuDgAAwLk4TaC3XQzrqpSDdrRtW4YZYyItywq46WEvSXfSk59uy5cvz4mnyRXWrFmjMWPG6I033tCECRNUunRpe5cEAABgd04T6HX7C2UlJYT+TF7IGm5Z1mAlvDGoohsuuoV9HT58WGvXrtUrr7wib29vbd++XTVq1LB3WQAAAA7DmQJ9trH1z0emOdDG1qPvI0mlS5dWREREpp/7TvbNza5fv67Q0FAtWLBA+fPnV7FixZLuGsvvDDnl4sWL/P8GOBHOWeRVBPpMsPXoh0iSp6enady4caaPdSf75la//vqr+vbtq507d6ply5aaNm2aKlWqZO+ykAdFRERwjgJOhHMWeZUzBfo0L1J1lHXj73QVl7zs6NGjatq0qcqVK6dvvvlGLVu2tHdJAAAADi2fvQtIL1tYj1bKy0m62bY5hONHDssYk+xr9erVtzyW2dVecpv4+HiFhyfcCiAxyCfOzgMAAOD2nCbQ22xSyhfHVpF0882h4AS2bNmiJ554Qt7e3tqwYYMk6dlnn03qlwcAAMDtOVugX6yUl5P0krQoh2vBHTh37pzefvtteXp66uDBg/rkk0/06KOP2rssAAAAp+OIPfSuSmVNeWNMiGVZvpZluSUuK2m7MdSZxLu+wvHFxcWpfv362rNnj3r06KGxY8dyl1cAAIBMcphAb1lWoBKCfDtJrpZlLVbChbDBtmUlEzWVNNSyrAO2n6sYY3LkJlC4MwcPHtSDDz4oFxcXjR49Wg8++KA8PT3T3hEAAACpcphAb4zxt33rm8a4aEn+txsDxxITE6Nx48YpKChIc+bMUadOnfTyyy/buywAAIBcwWECPXKnH3/8Ub169dLBgwf1+uuv65lnnrF3SQAAALmKs10UCyfSt29fPffccypUqJBWrVqlTz/9VGXKlLF3WQAAALkKM/TIUteuXVN8fLwKFSqkFi1aqHTp0howYIAKFixo79IAAAByJWbokWXWrFkjDw8PjR07VpLUokULDR06lDAPAACQjQj0uGMnT55U586d9eSTT+r8+fOsXAMAAJCDaLnBHfnuu+/01ltv6cKFC/L399e7777LXV4BAAByEIEemWKMkWVZevDBB1WnTh1Nnz5dDz/8sL3LAgAAyHMI9MiQ8+fP691339X58+c1f/581apVSytXrrR3WQAAAHkWPfRIF2OMFi1aJHd3d82YMUOFCxdWXFycvcsCAADI8wj0SNNff/2l5s2b65VXXtEDDzyg9evX64MPPpCLi4u9SwMAAMjzaLlBmgoWLKi9e/dq5syZ6t69O0EeAADAgRDokaKffvpJX3zxhRYsWKCyZctq//79KlCggL3LAgAAwE1ouUEyR44cUZs2bfTMM89ow4YNOn78uCQR5gEAABwUgR6SpGvXrmnSpElyd3fXDz/8oLFjx2rbtm0qW7asvUsDAADAbdByA0kJgX7mzJlq3LixZsyYocqVK9u7JAAAAKQDM/R52KlTp+Tv76+YmBgVKVJE69ev13fffUeYBwAAcCIE+jwoPj5e8+bNU9WqVTV58mStWbNGknT//ffLsiw7VwcAAICMINDnMdu2bVPDhg3VrVs31ahRQ1u3bpW3t7e9ywIAAEAm0UOfx/Tu3Vv79+/XggUL9MYbbzAjDwAA4OQI9LmcMUaLFy9Wo0aNVLp0aX388cdydXVVyZIl7V0aAAAAsgAtN7nYvn371KJFC7Vv314zZ86UJLm5uRHmAQAAchFm6HOhK1euaPz48Ro/frwKFSqkGTNmqEePHvYuCwAAANmAQJ8LDRkyRNOmTVOHDh00adIkbg4FAACQixHoc4mjR4/qypUrqlKligYPHqznn39eXl5e9i4LAAAA2Yweeid3/fp1TZ48We7u7vLz85MkPfDAA4R5AACAPIJA78R+//131a1bVwMGDNBTTz2l2bNn27skAAAA5DACvZP6+uuv1bBhQ509e1ZLly7V999/Lzc3N3uXBQAAgBxGoHci8fHxOnLkiCTpmWee0ejRo7Vr1y69+OKL3CAKAAAgjyLQO4lt27apYcOGaty4sa5cuaIiRYpo+PDhuvvuu+1dGgAAAOyIQO/gLly4oP79+6tu3brat2+fhg8frkKFCtm7LAAAADgIlq10YFFRUXryySd17Ngx+fj4aNy4cdzlFQAAAMkQ6B1QTEyMChcurEqVKumZZ55R165d9fjjj9u7LAAAADggWm4cyJUrVzRy5EhVqVJF//77r/Lly6d58+YR5gEAAJAqZugdxIoVK9SzZ0/t379fHTp0sHc5AAAAcBLM0NvZtWvX1L59ezVv3lz58uVTWFiYvvjiC91///32Lg0AAABOgEBvJ8YYSVKBAgVUuHBhjR49Wn/88Ye8vLzsXBkAAACcCYHeDtauXavHHntMu3btkiTNnz+f5SgBAACQKQT6HHT69Gl17dpVDRo00PHjx/Xvv/9KEnd5BQAAQKZxUWwO+emnn9SmTRudO3dOgwYN0ogRI7jLKwAAAO4YgT6HHD58WNWqVdPs2bNVo0YNe5cDAACAXIJAn0M6deokLy8v2msAAACQpeihzyEFChQgzAMAACDLEegBAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAnRqAHAAAAnBiBHgAAAHBiBHoAAADAiRHoAQAAACdGoAcAAACcGIEeAAAAcGIEegAAAMCJEegBAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAnRqAHAAAAnBiBHgAAAHBiBHoAAADAiRHoAQAAACdGoAcAAACcGIEeAAAAcGIEegAAAMCJEegBAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAnRqAHAAAAnBiBHgAAAHBiBHoAAADAiRHoAQAAACdGoAcAAACcGIEeAAAAcGIEekmWZS22LMvDsizXG7/sXRcAAACQlvz2LsBBeElqc/ODlmWVMMZE53w5AAAAQPowQ5/A3xhjJX5JqiKpLWEeAAAAji7PB3pba81XNz3cxhgTaodyAAAAgAxxuJYby7LaSGpvjGmbynZXSUMlHbA9VMUY45/Z57t5Ft6yLB9JIZk9HgAAAJCTHCbQW5YVbPvWTVLJ2wxdLMnXGBNl28/NsqwwY4x3FtTgKsmVVhsAAAA4C4dpuTHG+BpjfJUQ2FNkm72PSgzztv2ibth2p4ZKCs+C4wAAAAA5wmFm6NOpvaSwFB4Pk+QrKVRKapupm8axwlLok28jKTilwQAAAIAjcrZA76WUA3eUJM/EH4wxGe6Bt7XbuN04+w8AAAA4OodpuUlLYn+7pDMpbI62bbsTbne4PwAAAJDjnCbQ6/YXykpKCv13IvIO9wcAAABylLO13GQbY0yk0u67l5TUo+9j+/GKZVk707FbKUmnMlleXnCPpHP2LiId7FFndj5nVh77To+V2f0zul9GxnPepo5z1j7PyTl7e5yzqeOctc9zZuWxH0p1izHGob6UEJQ3p/C4myQjySOFbV62ba52qDckneM22ft368hf6f092vvLHnVm53Nm5bHv9FiZ3T+j+2VkPOdt9v155+Y6OWezdj/OWcf4887NdeaGc9aZWm5S6p1Pxthn/fjv7PCcuZGz/B7tUWd2PmdWHvtOj5XZ/TO6n7P8v+bonOX3yDmbfcfinHUuzvJ75JzNxLEsW+J3GLZ2Fl9jzC3tL5ZlnZXUzdy03KRtn0BjTIkcKjPDLMvaZIzxvM32wUpYraekEtbaZz18wM5ud95aluUmyV8pL4ELwA5SO2dt19i1U8ICGo9KCubfWeQmztZDv0kpXxxbRY5/Q6hUl9K03SV3ceJfLpZlhdn+UorOqeIApCjF89ayLC/bt2nd2RpAzkrt39qhxhh/KSncn7Usq4phqWrkEs7UciMl3EXWO4XHvSQtyuFaMsTcfm38djfNFIQpYSYBgB2ldt4aY8Jt52x0zlYE4HZSOmdtn6a53jAmWgk3ovTPscKAbOaIgd5VqawpbztR3WwnpyTJsiwPSWec9SPvxPpvejhaKb9xAQAAGdcuhaWt+XQNuYbDtNxYlhWohCDfTpKrZVmLlRB0g03CkpKJmkoaalnWAdvPVYwxdgu/lmW1kdTeGNM2le2ukoZKurHeG2cF3HTrLN8ZcaMrIFtkwTkLIAfd6Tlra6u5+Ro7D6V853nAKTlMoL/h5PNNY1y0HOBjMlvfu5R2D+1iJVzkG2Xbz82yrLB0vAlxvfMqASTKgXMWQBbKrnPW9sm4jDFBWVkvYE+O2HLjFIwxvsYYXyX8RZIi26xC1I0X3dzwF06bG4a63rRrSdGbC2SpLD5nAWSzbDxn54q2VuQyBPrs1V7S5hQeD9P/PomI1q0zD65KWMISQM5KzzkLwHFk6Jy1tfd2Y3Ub5DYE+uzlpZSDeZQkTylhtYwUtt+rhL+MAOSsNM9ZAA4l3ees7Z41ixKvy0tsvQFyAwJ9NrFdpOOqlO9wG63kbTbhN/3F4iHpq2wqDUAKMnjOArCzjJyztntHnJEUZVmWq21f3qQj13CYi2JzoTSXw7Isy9V2kW83JazcU1IJYT6Qm0oBOS7d56ztDbiX7cvNFg5COG+BHJWuc9Y2LqVPvWmjQ65BoHcAN63c4+h3vAXyPNtH9pGSWCUDcHC2fnnL3nUA2YmWGwAAAMCJEeizT0o9fcnw8TzgUDhnAefCOQvYEOizie0vkWilfMfXlO4OC8COOGcB58I5C/wPgT57bVLKF+1UEb3ygCPinAWcC+csIAJ9dluslO9G5yVpUQ7XAiBtnLOAc+GcBUSgzwquSmV9amNMiBKWtEv6ONC23N0ZY0xojlQH4Gau4pwFnImrOGeB27KMMfauwSnZbh/tKqmd7b+hSrhAJzjxLnS2ca6Shko6YHuoijHGXwByFOcs4Fw4Z4H0I9ADAAAAToyWGwAAAMCJEegBAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAnRqAHAAAAnBiBHgAAAHBiBHoAAADAiRHoAQAAACdGoAcAAACcGIEeAAAAcGIEegAAAMCJ5bd3AQCA3MuyLA9JXpK8JfkaY6Isyxps21xFUkljTFu7FQgAuQCBHgCQndobY/wty6oiKdiyrEhjjH/iRsuyDliW5WOMCbFjjQDg1Aj0AIBsYZud32j70U2Sp6SbZ+OjlTBTDwDIJHroAQDZxhgTavvWU1KAMSb6piEekg7kaFEAkMsQ6AEA2cIYEylJlmW5SXKVFH7jdtsMvm5+HACQMQR6AEB285L+F/Bv0F5SlDEmKudLAoDcg0APAMhu3kp5Ft5HUqiUNIsPAMgEAj0AILt5SQq78QFbu42rpGDbQ21yuCYAyDUI9ACAbJNa/7ykkpKibevSu0mi7QYAMskyxti7BgBALmVZlpekYGPMLUtTWpYVpoSZ+2jWoQeAzCPQAwAAAE6MlhsAAADAiRHoAQAAACdGoAcAAACcGIEeAAAAcGIEegAAAMCJEegBAAAAJ0agBwAAAJwYgR4AAABwYgR6AAAAwIkR6AEAAAAn9v/fM7AT1QrEGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams.update({'font.size': 22})\n",
    "\n",
    "plt.figure(figsize=(12, 8))\n",
    "\n",
    "plt.loglog(sizes, [i for i in runtimes], \"-ks\", label=r\"Computed runtime\", markersize=10, markerfacecolor=(0, 0.447, 0.741, 1))\n",
    "plt.loglog(sizes, 1e-6*np.power(sizes, 3.0), \"--k\", label=r\"$n^3$\")\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "\n",
    "plt.xlabel(r\"$n$\")\n",
    "plt.ylabel(r\"time $(s)$\")\n",
    "plt.xlim([0.9, 600])\n",
    "plt.ylim([1e-7, 1e2])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### A recursive approach to matrix multiplication\n",
    "\n",
    "Assuming that $n=2^m$, divide $\\mathbf{A}$, $\\mathbf{B}$ and $\\mathbf{C}$ into smaller matrices of size $\\frac{n}{2} \\times \\frac{n}{2}$,\n",
    "$$\n",
    "\\mathbf{A} = \\left(\n",
    "\\begin{array}{cc}\\mathbf{A}_{11}&\\mathbf{A}_{12}\\\\\n",
    "\\mathbf{A}_{21}&\\mathbf{A}_{22}\n",
    "\\end{array} \\right)\n",
    "\\hspace{0.5cm}\n",
    "\\mathbf{B} = \\left(\n",
    "\\begin{array}{cc}\\mathbf{B}_{11}&\\mathbf{B}_{12}\\\\\n",
    "\\mathbf{B}_{21}&\\mathbf{B}_{22}\n",
    "\\end{array} \\right)\n",
    "\\hspace{0.5cm}\n",
    "\\mathbf{C} = \\left(\n",
    "\\begin{array}{cc}\\mathbf{C}_{11}&\\mathbf{C}_{12}\\\\\n",
    "\\mathbf{C}_{21}&\\mathbf{C}_{22}\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "The matrices $C_{i\\,j}$ are\n",
    "$$\n",
    "\\nonumber \\mathbf{C}_{11} = \\mathbf{A}_{11}\\mathbf{B}_{11}+\\mathbf{A}_{12}\\mathbf{B}_{21}\\\\\n",
    "\\nonumber \\mathbf{C}_{12} = \\mathbf{A}_{11}\\mathbf{B}_{21}+\\mathbf{A}_{12}\\mathbf{B}_{22}\\\\\n",
    "\\label{eq-DC1} \\mathbf{C}_{21} = \\mathbf{A}_{21}\\mathbf{B}_{11}+\\mathbf{A}_{22}\\mathbf{B}_{21}\\\\\n",
    "\\nonumber \\mathbf{C}_{22} = \\mathbf{A}_{21}\\mathbf{B}_{21}+\\mathbf{A}_{22}\\mathbf{B}_{22}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Divide and conquer strategy\n",
    "\n",
    "Multiplication of $n \\times n$ matrices can be expressed as multiplication (and addition) of $\\frac{n}{2} \\times \\frac{n}{2}$ matrices. Recursion leads to trivial $1 \\times 1$ matrices.  \n",
    "\n",
    "This is an example of what is called a \"divide and conquer\" strategy. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split"
   },
   "source": [
    "![Float32](files/images/recursionTree.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Divide-and-conquer**: break a problem into smaller subproblems of the same type, recursively solve the subproblems, and finally combine the solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computational complexity of recursive matrix multiplication\n",
    "\n",
    "Is the recursive strategy better or worse than the iterative one?\n",
    "\n",
    "Let $F(n)$ be the ops required to multiply 2 matrices of size $n \\times n$.  \n",
    "We must have\n",
    "\n",
    "$$ F(n) =  8\\, F\\left(\\frac{n}{2}\\right) + 4\\,\\left(\\frac{n}{2}\\right)^2, $$\n",
    "\n",
    "since we have to do 8 multiplications and 4 additions of size $\\frac{n}{2}$.\n",
    "\n",
    "When the recursion reaches $n=1$, we have $F(1) = 1$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving for $F(n)$\n",
    "\n",
    "How do we solve such equations?  \n",
    "Trick: let  \n",
    "> $n = 2^p$  \n",
    "> $a_p = F(2^p).$ \n",
    "\n",
    "This gives a linear first order recursion relation\n",
    "$$ a_p = 8\\, a_{p-1} + 2^{2p}$$ \n",
    "\n",
    "with $a_0 = 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The solution is\n",
    "$$a_p = 2^{2p}\\left(2^{p+1}-1 \\right).$$\n",
    "Returning to the original variables: $p = \\log_2 n$.\n",
    "\n",
    "$$ F(n) = n^2 \\left(2n - 1 \\right).$$\n",
    "\n",
    "In this case, the divide and conquer strategy gives the same answer as the iterative strategy ... but this is not always the case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The \"Master Theorem\"\n",
    "\n",
    "Let $a$ and $b$ be integers greater than or equal to 1 (usually $b=2$). Let $c$ and $d$ be positive and real. \n",
    "\n",
    "Given a recurrence relation of the form\n",
    "$$\n",
    "F(n) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "a\\, F(n/b) + n^c&\\mbox{if $n>1$}\\\\\n",
    "d&\\mbox{if $n=1$,}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "then if $n$ is a power of $b$, there are three cases:\n",
    "1. if $\\log_b(a) < c$ then $F(n) = O(n^c)$\n",
    "2. if $\\log_b(a) = c$ then $F(n) = O(n^c \\log(n))$\n",
    "3. if $\\log_b(a) > c$ then $F(n) = O(n^{\\log_b(a)})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Matrix multiplication again: Strassen's algorithm\n",
    "\n",
    "An alternative (very non-obvious) divide and conquer strategy is:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mathbf{C}_{11} = \\mathbf{M}_{1}+\\mathbf{M}_{4}-\\mathbf{M}_{5}+\\mathbf{M}_{7} & \\mathbf{C}_{12} =\\mathbf{M}_{3}+\\mathbf{M}_{5}\\\\\n",
    "\\mathbf{C}_{21} = \\mathbf{M}_{2}+\\mathbf{M}_{4} & \\mathbf{C}_{22} = \\mathbf{M}_{1}-\\mathbf{M}_{2}+\\mathbf{M}_{3}+\\mathbf{M}_{6},\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{M}_1$ to $\\mathbf{M}_7$ are\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mathbf{M}_{1} = \\left(\\mathbf{A}_{11}+\\mathbf{A}_{22}\\right)\\left(\\mathbf{B}_{11}+\\mathbf{B}_{22}\\right) &\n",
    "\\mathbf{M}_{2} = \\left(\\mathbf{A}_{21}+\\mathbf{A}_{22}\\right) \\mathbf{B}_{11}\\\\\n",
    "\\mathbf{M}_{3} = \\mathbf{A}_{11} \\left(\\mathbf{B}_{12}-\\mathbf{B}_{22}\\right)&\n",
    "\\mathbf{M}_{4} =  \\mathbf{A}_{22} \\left(\\mathbf{B}_{21}-\\mathbf{B}_{11}\\right)\\\\\n",
    "\\mathbf{M}_{5} =  \\left(\\mathbf{A}_{11}+\\mathbf{A}_{12}\\right) \\mathbf{B}_{22}&\n",
    "\\mathbf{M}_{6} = \\left(\\mathbf{A}_{21}-\\mathbf{A}_{11}\\right)\\left(\\mathbf{B}_{11}+\\mathbf{B}_{12}\\right)\\\\\n",
    "\\mathbf{M}_{7} =  \\left(\\mathbf{A}_{12}-\\mathbf{A}_{22}\\right)\\left(\\mathbf{B}_{21}+\\mathbf{B}_{22}\\right)&\n",
    " .\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computational complexity of Strassen multiplication\n",
    "\n",
    "These crazy formulae were dreamt up by Strassen (1969). The key point is that there are only 7 multiplications rather than 8. We have:\n",
    "$$ F(n) = 7\\,F\\left(\\frac{n}{2}\\right) + 18\\,\\left(\\frac{n}{2}\\right)^2$$,\n",
    "with $F(1)=1$. The solution is\n",
    "\n",
    "$$ F(n) = 7 n^{\\log_2(7)} - 6 n^2 = \\mathcal{O}(n^{\\log_2(7)}).$$\n",
    "\n",
    "Check that this is consistent with the master theorem.\n",
    "\n",
    "Note that $\\log_2(7) \\approx 2.807 < 3$. Amazing!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving linear recursion relations\n",
    "\n",
    "We don't need to memorise the master theorem: linear recursion relations are very similar to linear differential equations.\n",
    "\n",
    "1st order case:\n",
    "$$a_{n} = b_1\\, a_{n-1} + f(n)\\ \\ \\ \\ \\  \\text{with initial condition}\\ \\  a_1 = A_1.$$ \n",
    "2nd order case:\n",
    "$$a_{n} = b_1\\, a_{n-1}+ b_2\\, a_{n-2} + f(n)\\ \\ \\ \\ \\  \\text{with initial conditions}\\ \\  a_1 = A_1,\\ a_2=A_2 $$ \n",
    "\n",
    "Solution procedure is the same in both cases:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving linear recursion relations\n",
    "1. Find the *general solution* of the homogeneous equation,\n",
    "$$ a_{n} = b_1\\, a_{n-1}\\ \\ \\  \\text{or}\\ \\ \\ a_{n} = b_1\\, a_{n-1}+ b_2\\, a_{n-2}$$\n",
    "with the ansatz $a_n = x^n$.  \n",
    "The solution will usually be of the form\n",
    "$$ a_n = C_1\\,x^n\\ \\ \\text{or}\\ \\ a_{n} = C_1\\,x_1^n + C_2\\,x_2^n,$$\n",
    "where the $C$'s are arbitrary constants.\n",
    "2. Find a *particular solution*, $\\alpha_n$, of the inhomogeneous equation.\n",
    "3. The full solution is\n",
    "$$ a_n = C_1\\,x^n +\\alpha_n\\ \\ \\text{or}\\ \\ a_{n} = C_1\\,x_1^n + C_2\\,x_2^n + \\alpha_n,$$\n",
    "Use the initial conditions to evaluate the constants, $C$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solving linear recursion relations\n",
    "The tricky part is step 2 - need to \"guess\" a particular solution.\n",
    "\n",
    "A good guess is often a term, $\\text{const} \\times f(n)$, proportional to $f(n)$ (or sometimes a polynomial in $f(n)$).\n",
    "\n",
    "As for resonant ODEs, it can happen that the solution of the homogeneous equation is already proportional to $f(n)$\n",
    "\n",
    "In this case a good guess for the particular solution is $\\text{const}\\times n\\,f(n)$.\n",
    "\n",
    "See the excellent online notes available __[here](https://www.tutorialspoint.com/discrete_mathematics/discrete_mathematics_recurrence_relation.htm)__."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}