Practice 06 - SLIM MSE with Gradient Descent and Cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2021/20\n",
    "\n",
    "### Practice - SLIM MSE implemented with Python and Cython\n",
    "\n",
    "\n",
    "### Cython is a superset of Python, allowing you to use C-like operations and import C code. Cython files (.pyx) are compiled and support static typing.\n",
    "\n",
    "### Why do we use it (or any other compiled language)? If the code is written properly it is fast... I mean, FAST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's implement something simple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isPrime(n):\n",
    "    \n",
    "    i = 2\n",
    "    \n",
    "    # Usually you loop up to sqrt(n)\n",
    "    while i < n:\n",
    "        if n % i == 0:\n",
    "            return False\n",
    "        \n",
    "        i += 1\n",
    "        \n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is 2 prime? True\n",
      "Is 3 prime? True\n",
      "Is 5 prime? True\n",
      "Is 15 prime? False\n",
      "Is 20 prime? False\n"
     ]
    }
   ],
   "source": [
    "print(\"Is 2 prime? {}\".format(isPrime(2)))\n",
    "print(\"Is 3 prime? {}\".format(isPrime(3)))\n",
    "print(\"Is 5 prime? {}\".format(isPrime(5)))\n",
    "print(\"Is 15 prime? {}\".format(isPrime(15)))\n",
    "print(\"Is 20 prime? {}\".format(isPrime(20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is 80000023 prime? True, time required 8.87 sec\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "result = isPrime(80000023)\n",
    "\n",
    "print(\"Is 80000023 prime? {}, time required {:.2f} sec\".format(result, time.time()-start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load Cython magic command, this takes care of the compilation step. If you are writing code outside Jupyter you'll have to compile using other tools. See at the end of the notebook for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Declare Cython function, paste the same code as before. The function will be compiled and then executed with a Python interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def isPrime(n):\n",
    "    \n",
    "    i = 2\n",
    "    \n",
    "    # Usually you loop up to sqrt(n)\n",
    "    while i < n:\n",
    "        if n % i == 0:\n",
    "            return False\n",
    "        \n",
    "        i += 1\n",
    "        \n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is Prime 80000023? True, time required 4.73 sec\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "result = isPrime(80000023)\n",
    "\n",
    "print(\"Is Prime 80000023? {}, time required {:.2f} sec\".format(result, time.time()-start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As you can see by just compiling the same code we got some improvement.\n",
    "#### To go seriously higher, we have to use some static typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "# Declare the tipe of the arguments\n",
    "def isPrime(long n):\n",
    "    \n",
    "    # Declare index of for loop\n",
    "    cdef long i\n",
    "    \n",
    "    i = 2\n",
    "    \n",
    "    # Usually you loop up to sqrt(n)\n",
    "    while i < n:\n",
    "        if n % i == 0:\n",
    "            return False\n",
    "        \n",
    "        i += 1\n",
    "        \n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is 80000023 prime? True, time required 1.20 sec\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "result = isPrime(80000023)\n",
    "\n",
    "print(\"Is 80000023 prime? {}, time required {:.2f} sec\".format(result, time.time()-start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cython code with two type declarations, for n and i, runs 50x faster than Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Main benefits of Cython:\n",
    "* Compiled, no interpreter\n",
    "* Static typing, no overhead\n",
    "* Fast loops, no need to vectorize. Vectorization sometimes performes lots of useless operations\n",
    "* Numpy, which is fast in python, when operations are not vectorizable often becomes slooooow compared to a carefully written Cython code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A few info about gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from scipy import stats\n",
    "from scipy.optimize import fmin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Gradient descent</b>, also known as <b>steepest descent</b>, is an optimization algorithm for finding the local minimum of a function. To find a local minimum, the function \"steps\" in the  direction of the negative of the gradient. <b>Gradient ascent</b> is the same as gradient descent, except that it steps in the direction of the positive of the gradient and therefore finds local maximums instead of minimums. The algorithm of gradient descent can be outlined as follows:\n",
    "\n",
    "&nbsp;&nbsp;&nbsp; 1: &nbsp; Choose initial guess $x_0$ <br>\n",
    "&nbsp;&nbsp;&nbsp;    2: &nbsp; <b>for</b> k = 0, 1, 2, ... <b>do</b> <br>\n",
    "&nbsp;&nbsp;&nbsp;    3:   &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $s_k$ = -$\\nabla f(x_k)$ <br>\n",
    "&nbsp;&nbsp;&nbsp;    4:   &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; choose $\\alpha_k$ to minimize $f(x_k+\\alpha_k s_k)$ <br>\n",
    "&nbsp;&nbsp;&nbsp;    5:   &nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{k+1} = x_k + \\alpha_k s_k$ <br>\n",
    "&nbsp;&nbsp;&nbsp;    6: &nbsp;  <b>end for</b>\n",
    "\n",
    "As a simple example, let's find a local minimum for the function $f(x) = x^3-2x^2+2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x**3-2*x**2+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGiCAYAAADA0E3hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABMmUlEQVR4nO3deVxU9f4/8NcwMwwgMMi+Iwii4obggootFi43vy3WbbumbTf7udwi7y3te+vWvfdre9a1NG9q19t6Cy27mkUmuKGC4o64sYkgsg37DDNzfn8MTJGgDMzMmeX1fDzm8YjTOfKew9F58VklgiAIICIiIrIDLmIXQERERNRbDC5ERERkNxhciIiIyG4wuBAREZHdYHAhIiIiu8HgQkRERHaDwYWIiIjsBoMLERER2Q0GFyIiIrIbDC5ERERkN0wKLqtXr8aoUaPg7e0Nb29vpKSk4LvvvrvmNdnZ2UhKSoKbmxtiYmKwZs2afhVMREREzsuk4BIeHo5XXnkFeXl5yMvLw80334zbb78dJ0+e7Pb8oqIizJo1C6mpqcjPz8fy5cuxZMkSZGRkmKV4IiIici6S/m6y6Ovri9dffx2PPvroVf/v2WefxZYtW1BQUGA8tmDBAhw9ehQ5OTn9+bZERETkhGR9vVCn0+HLL79Ec3MzUlJSuj0nJycHaWlpXY5Nnz4d69atQ3t7O+RyebfXqdVqqNVq49d6vR61tbXw8/ODRCLpa8lERERkRYIgoLGxEaGhoXBxMc+wWpODy/Hjx5GSkoK2tjZ4enpi8+bNGD58eLfnVlZWIigoqMuxoKAgaLVaVFdXIyQkpNvrVqxYgZdeesnU0oiIiMgGlZWVITw83Cx/lsnBJT4+HkeOHEF9fT0yMjIwb948ZGdn9xheft1C0tkzda2Wk2XLliE9Pd34tUqlQmRkJMrKyuDt7W1qyURERA7jh5OVSP/PUQwP8cZ/FnTf42ErGhoaEBERAS8vL7P9mSYHF1dXV8TGxgIAkpOTkZubi3feeQcffPDBVecGBwejsrKyy7GqqirIZDL4+fn1+D0UCgUUCsVVxztnMxERETmrC6pLcFF4IDE21G4+E805zKPfHU6CIHQZj/JLKSkpyMzM7HLshx9+QHJyco/jW4iIiKhnx8tVAIARYfYRWszNpOCyfPly7N69G8XFxTh+/Dief/55ZGVl4cEHHwRg6OJ56KGHjOcvWLAAJSUlSE9PR0FBAdavX49169Zh6dKl5n0XRERETkAQBJy81AAAGBmmFLkacZjUVXT58mXMnTsXFRUVUCqVGDVqFLZv345bb70VAFBRUYHS0lLj+dHR0di2bRuefvppvPfeewgNDcW7776LOXPmmPddEBEROYFLqjbUNmsgc5FgSJD5xo3Yk36v42INDQ0NUCqVUKlUdtOfR0REZG7fn6zEE/8+hOEh3tj2h1Sxy7kuS3x+c68iIiIiO3HsYj0A5+0mAhhciIiI7Maxi4aBuaMiGFyIiIjIhgmCgKNl9QCA0eE+otYiJgYXIiIiO1Bc04KGNi0UMhfEBzvnwFyAwYWIiMgudI5vGR7qDbnUeT++nfedExER2ZGjZYbxLc7cTQQwuBAREdmFzhaXUeHOOzAXYHAhIiKyeVqdHicudcwoYosLERER2bIzl5vQ1q6Hl0KGGP8BYpcjKgYXIiIiG2dceC5cCRcX8+20bI8YXIiIiGzcUeP4Fh9R67AFDC5EREQ27ucZRc49MBdgcCEiIrJpbe06FF5uBACMivARtxgbwOBCRERkw05eaoBOL8DfU4FQpZvY5YiOwYWIiMiG/bw/kRISiXMPzAUYXIiIiGzaMQ7M7YLBhYiIyIYdu9ix8FwEB+YCDC5EREQ2S9XajgvVzQC4R1EnBhciIiIbdbyjtSXC1x2+A1xFrsY2MLgQERHZKC48dzUGFyIiIhvVOTCXC8/9jMGFiIjIRnWumMsWl58xuBAREdmgClUrKhva4CIBRrHFxYjBhYiIyAbll9YDAIYGe8PDVSZuMTaEwYWIiMgG5ZfWAQASI33ELcTGMLgQERHZoM4Wl7GRA8UtxMYwuBAREdkYjVaP4+WGgblscemKwYWIiMjGnK5sgFqrh9Jdjmj/AWKXY1MYXIiIiGxMZzdRYqQPd4T+FQYXIiIiG2McmBvB8S2/xuBCRERkY/LL6gFwfEt3GFyIiIhsSHWTGiU1LQCA0RE+4hZjgxhciIiIbMiRjvEtsYGeULrLxS3GBjG4EBER2ZD8ss7xLT7iFmKjGFyIiIhsiHHhuSgOzO0OgwsREZGN0OkFHOXA3GticCEiIrIRZ6sa0azRYYCrFHGBXmKXY5MYXIiIiGxEZzfR6AgfSF248Fx3GFyIiIhsBHeEvj4GFyIiIhthXOqfK+b2iMGFiIjIBqha23G2qgkAMIYtLj1icCEiIrIBnbOJIn094O+pELcYG8bgQkREZAPySgzjW8ayteWaGFyIiIhswKGSWgBA8iBfkSuxbQwuREREItPq9MaBucmDODD3WhhciIiIRHa6shEtGh283GQYwoXnronBhYiISGR5xYZuorGRA+HCheeuicGFiIhIZLkdA3OTubHidTG4EBERiUgQBBwq7gguHJh7XQwuREREIiqvb0VlQxtkLhKMifARuxybx+BCREQkokMd3UQJod5wd5WKXI3tMym4rFixAuPGjYOXlxcCAwNxxx13oLCw8JrXZGVlQSKRXPU6ffp0vwonIiJyBHkd3URJUewm6g2Tgkt2djYWLlyI/fv3IzMzE1qtFmlpaWhubr7utYWFhaioqDC+4uLi+lw0ERGRo8gt7lx4jgNze0Nmysnbt2/v8vWGDRsQGBiIQ4cOYerUqde8NjAwED4+Pr36Pmq1Gmq12vh1Q0ODKWUSERHZhYa2dhRebgTAGUW91a8xLiqVCgDg63v95q3ExESEhIRg2rRp2Llz5zXPXbFiBZRKpfEVERHRnzKJiIhsUn5pPQTBsLFioLeb2OXYhT4HF0EQkJ6ejilTpmDEiBE9nhcSEoK1a9ciIyMDmzZtQnx8PKZNm4Zdu3b1eM2yZcugUqmMr7Kysr6WSUREZLMOdXYTsbWl10zqKvqlRYsW4dixY9izZ881z4uPj0d8fLzx65SUFJSVleGNN97osXtJoVBAoeCW3kRE5Ng6d4RO4viWXutTi8vixYuxZcsW7Ny5E+Hh4SZfP3HiRJw9e7Yv35qIiMghtP9yY0XOKOo1k1pcBEHA4sWLsXnzZmRlZSE6OrpP3zQ/Px8hISF9upaIiMgRFFQ0oLVdB283GeICPcUux26YFFwWLlyITz/9FN988w28vLxQWVkJAFAqlXB3dwdgGJ9SXl6OjRs3AgBWrlyJQYMGISEhARqNBh9//DEyMjKQkZFh5rdCRERkP35ev4UbK5rCpOCyevVqAMCNN97Y5fiGDRswf/58AEBFRQVKS0uN/0+j0WDp0qUoLy+Hu7s7EhISsHXrVsyaNat/lRMREdmxzhVzuT+RaSSCIAhiF3E9DQ0NUCqVUKlU8Pb2FrscIiKifhEEAeP/bweuNKrx+e8nYmKMn9glWYQlPr+5VxEREZGVFVU340qjGq4yF26saCIGFyIiIis7WGRYv2VMhA/c5NxY0RQMLkRERFbWGVwmRHN8i6kYXIiIiKzsgDG4OObYFkticCEiIrKii3UtKK9vhcxFgrFRPmKXY3cYXIiIiKzowAVDa8vIcCU8XPu8847TYnAhIiKyos7xLeM5vqVPGFyIiIis6EBRDQBgIse39AmDCxERkZVcbmhDcU0LJBLuCN1XDC5ERERW0jmbaHiIN7zd5CJXY58YXIiIiKzkYEc3EadB9x2DCxERkZVwYG7/MbgQERFZQW2zBmcuNwFgcOkPBhciIiIr6GxtGRLkCd8BriJXY78YXIiIiKzgAMe3mAWDCxERkRV0rpjLbqL+YXAhIiKysLpmDU5VNAAAJsawxaU/GFyIiIgsbP8FQzfRkCBPBHgpRK7GvjG4EBERWVhOR3CZNNhf5ErsH4MLERGRhe07bwguKYPZTdRfDC5EREQWVNXQhnNVTZBIuLGiOTC4EBERWVBnN1FCqDeUHtyfqL8YXIiIiCwo5zzHt5gTgwsREZEFcXyLeTG4EBERWUhZbQtKa1sgc5Fg3CAuPGcODC5EREQW0jm+ZVS4Ep4KmcjVOAYGFyIiIgvZz/EtZsfgQkREZAGCIBjHt0zi+BazYXAhIiKygKLqZlQ2tMFV5oKxUQPFLsdhMLgQERFZQGdry9hIH7jJpSJX4zgYXIiIiCyA+xNZBoMLERGRmen0AvaeqwbA8S3mxuBCRERkZicvqVDf0g4vhQxjInzELsehMLgQERGZ2e6zhtaWiYP9IJPyo9aceDeJiIjMbE9HcJkax/Et5sbgQkREZEYtGi3ySmoBAFPiAkSuxvEwuBAREZnRgaJatOsEhPm4Y5Cfh9jlOBwGFyIiIjPq7CZKjfOHRCIRuRrHw+BCRERkRrvPXgEApLKbyCIYXIiIiMzkckMbzlxugkTC9VsshcGFiIjITDq7iUaGKTFwgKvI1TgmBhciIiIz+bmbiNOgLYXBhYiIyAwEQcCec4b9iabEcnyLpTC4EBERmcHpykZUN6nhLpdibJSP2OU4LAYXIiIiM+jsJpoQ4wuFTCpyNY6LwYWIiMgMOvcnmhLL8S2WxOBCRETUTy0aLQ5cMCzzf2M8x7dYEoMLERFRP+Wcr4FGp0f4QHcMDvAUuxyHxuBCRETUTzsLqwAYWlu4zL9lMbgQERH1gyAIyCo0DMy9cUigyNU4PpOCy4oVKzBu3Dh4eXkhMDAQd9xxBwoLC697XXZ2NpKSkuDm5oaYmBisWbOmzwUTERHZkvNXmnGxrhWuUhdMiuUy/5ZmUnDJzs7GwoULsX//fmRmZkKr1SItLQ3Nzc09XlNUVIRZs2YhNTUV+fn5WL58OZYsWYKMjIx+F09ERCS2rI5uogkxvvBwlYlcjeMz6Q5v3769y9cbNmxAYGAgDh06hKlTp3Z7zZo1axAZGYmVK1cCAIYNG4a8vDy88cYbmDNnTt+qJupBW7sO9S3taNfpodULAIABCim8FHK4yV3Y90xEZtfZTXTDEM4msoZ+RUOVSgUA8PX17fGcnJwcpKWldTk2ffp0rFu3Du3t7ZDL5Vddo1aroVarjV83NDT0p0xyMIIgoKy2FYdKa3HmchPOVDaiqKYZVxrVaGzT9nidq9QFoT5uCB/ogQhfDwwP9UZCqDeGBXvD3ZWLRRGR6ZrVWhws6pwGzfEt1tDn4CIIAtLT0zFlyhSMGDGix/MqKysRFBTU5VhQUBC0Wi2qq6sREhJy1TUrVqzASy+91NfSyAG1aLT46XQVdhRUYf+FGlSo2no8V+YigVzqAplUAghAs0YLvQBodHoU17SguKblqvMTI30wabA/psT5Y0yED+RSjlsnouvrOg16gNjlOIU+B5dFixbh2LFj2LNnz3XP/XXzvCAI3R7vtGzZMqSnpxu/bmhoQERERF9LJTul0wvIPlOFL/MuYmdhFdra9cb/J5dKMDJMieGh3hgS5IXYAE8EKd0Q4KWAl0LW5dkSBAGt7TrUNmtQXteKi3WtuFDdhJOXGnCivAHVTWrkFtcht7gO7+w4i4EecswYEYzbRoViQrQvZAwxRNSDrDOG8S03xQeyK9pK+hRcFi9ejC1btmDXrl0IDw+/5rnBwcGorKzscqyqqgoymQx+ft2PvlYoFFAoFH0pjRyAqqUdnx4sxScHSnCxrtV4PNLXAzNHBmNqXADGRg7sdfeORCKBh6sMHq4yhA/0wIRf/f/SmhbsPV+NvecMr7qWdnx2sAyfHSyDv6cCv00Ox/3jIxHh62HGd0lE9k4QBOw83TENmqvlWo1JwUUQBCxevBibN29GVlYWoqOjr3tNSkoKvv322y7HfvjhByQnJ3c7voWcl6qlHev2FmHDniI0qg1jVZTuctyTFI47EsOQEOptkd9oIv08EOkXifvHR0Kr0+NAUS3+e+wStp+oRHWTGu9nncfq7PO4YUgAHp4cjalx/vzNiohw/koTyusN06BTBnMatLWYFFwWLlyITz/9FN988w28vLyMLSlKpRLu7u4ADN085eXl2LhxIwBgwYIFWLVqFdLT0/H4448jJycH69atw2effWbmt0L2SqPVY/3eIrz30zljYIkP8sJjqdGYPToUbnLrDZyVSV0wOdYfk2P98fLtI7Cj4DI+OVCK3WerkVV4BVmFVzAizBv/78ZYTE8IhtSFAYbIWXXOJuI0aOuSCJ0DTnpzcg+/ZW7YsAHz588HAMyfPx/FxcXIysoy/v/s7Gw8/fTTOHnyJEJDQ/Hss89iwYIFvS6yoaEBSqUSKpUK3t7evb6ObF9WYRVe/vYULlQb1gKKD/LCU7fEYXpCMFxsKBQUVTdjY04xPj9YhtZ2HQBgcMAA/HF6PKYnBLMFhsgJPfDP/dh3vgb/+5theCw1RuxybJIlPr9NCi5iYXBxPPUtGrzwzUlsOXoJAODv6YpnZwzFnLHhNhVYfq22WYOP9hbho33FaOiYej020gfLZg3DuEE9LwtARI5F1dqOpL9mQqsXsHPpjYj254yi7lji85ttW2R1O09X4dmMY6hqVEPqIsH8SYPwh1vi4O1m+2OefAe4Ij0tHo9PjcHaXRfw4e4iHC6txz1rcnDr8CC8cNtwDuIlcgLZZ65AqxcQG+jJ0GJlDC5kNRqtHv+3rQAf7SsGYOhqefO3YzAmwkfUuvrCy02OZ9LiMXdiFN7+8Sz+k1eGzFOXsfvsFSy+OQ6Pp8bAVcZp1ESOakfBZQDAtGFcdM7a+C8rWUWFqhX3rc0xhpZHp0Rj65JUuwwtvxTo7YYVd43E90+lIiXGD23terz+fSFmvrML+85Xi10eEVlAu06PnacN67fcOizoOmeTuTG4kMUdLKrFbe/uweHSeni5yfDhQ8n4823DrTpbyNJiA73w6eMTsPLeMfD3dMX5K8144J8H8Pzm42hW97wNARHZn9ziWjS0aeE7wBWJkQPFLsfpMLiQRX1zpBy/+/AAapo1GB7ija2LU3HLcMf8DUUikeCOxDDseOZGPDghEgDwyYFSzHhnF/ZfqBG5OiIylx0FP6+WyyURrI/BhSxCEAS8t/Mc/vD5EWh0esxICEbGk5MQ6ef4A1eV7nL8/c6R+OSxCQjzcUdZbSvu/+d+/PW/p6DW6sQuj4j6QRAE/NgxvuUWjm8RBYMLmZ0gCHjp21N4/ftCAMBjU6Lx3oNjnW4H5smx/tj+VCruGxcBQQDW7SnCnNX7UNyxZg0R2Z/zV5pQUtMCV6kLUodwmX8xMLiQWen1Ap7/+oRxEO6Ls4fjf28b7rTNqV5ucrwyZxTWz0/GQA85TpQ34LZ/7ME3R8rFLo2I+iDzlKGbaOJgP3gqODFXDAwuZDY6vYA/fnUMnx4ohYsEeP3uUXh48vX3s3IGNw8NwrY/pGL8IF80qbX4w+dH8FzGMbRq2HVEZE86p0Hfym4i0TC4kFno9QKeyziGjMMXIXWR4O17x+Ce5Aixy7IpIUp3fPr4BCy5ORYSCfB5bhnmrN6Hi3UtYpdGRL1Q06TGodI6AMDNnAYtGgYX6jdBEPC3rQX48pAhtPzj/kTcPiZM7LJskkzqgvS0eHz86AT4DXDFqYoG/M+qvVzzhcgO/HS6CoIADA/xRpiPu9jlOC0GF+q3d3ecw/q9RQCA1+aMwqyRISJXZPsmx/pjy+IpGBmmRG2zBnPXHcSGvUWwg63DiJwWZxPZBgYX6peP95fg7R/PAAD+Mns45iSFi1yR/QjzcceXC1JwZ2IYdHrDTKylXx7jlGkiG9Si0SL7zBUAQFpCsMjVODcGF+qzrMIqvLjlJADgqVviMJ8DcU3mJpfird+Oxp87Zl5lHL6Ih9YdhKqlXezSiOgXdp25grZ2PcIHuiMh1Dy7HFPfMLhQn5yubMCiT/Oh0wuYMzYcf5gWJ3ZJdksikeDRKdHYMH8cPBUyHCiqxZ2r96K0hoN2iWzFdycqAQAzRwRDInHO5R1sBYMLmayqsQ2PbMhFk1qLiTG+WHHXSP5FNoOpQwLw1ZMpCFW64cKVZtz5/l7kd8xgICLxqLU6/NSxzP+MEewmEhuDC5mkXafHok/zcUnVhhj/AVjzuyS4yvgYmcvQYG9sXjgZCaHeqGnW4L61+7H9RIXYZRE5tX3na9Co1iLAS4HECG6qKDZ+4pBJXv3uNA4W1cJTIcM/5yXDx8NV7JIcTpC3G/7zRApuHhoItVaPJz85jE8OlIhdFpHT+r6jm2h6QhBcnHQVcFvC4EK9tvVYBT7cY5j2/MY9ozE4wFPkihzXAIUMa+cm4YEJkRAE4PnNJ/DeznOcLk1kZTq9gB9OGaZBzxzBpR5sAYML9cq5qkb86aujAIAnbohhP68VyKQu+PsdI7DoplgAwOvfF+LvWwug1zO8EFnLwaJa1DZr4OMhx/hoX7HLITC4UC+0tevw5MeH0azRISXGD39Mixe7JKchkUiwdHo8/vc3wwAAH+4pwh+/OgatTi9yZUTO4fuThm6iW4YFQS7lR6Yt4E+Bruv/thXgbFUTArwUePf+RMj4l9fqHkuNwZv3jDau9fLkJ4e5UB2Rhen1Arb/Yho02QZ+AtE17TxdhY05hoGhb94zGgFeCpErcl5zksLxwe+SoJC5IPPUZTzx70Noa2d4IbKUY+UqVDa0YYCrFJNj/cUuhzowuFCPrjSq8ceOcS2PTI7G1CEBIldEtwwPwvr54+Amd0FW4RU8vjEPrRqGFyJL+K5jKYKbhgbCTS4VuRrqxOBC3RIEAX/66iiqmzQYGuyFP83guBZbMTnWHx89PB4erlLsPluNhz86iGa1VuyyiByKIAjYeswQXLhxrG1hcKFufXawDDsLr8BV5oJ37kvkbxs2ZmKMHzY+Mh6eChn2X6jF/A0H0cTwQmQ2R8rqcbGuFR6uUtwUz92gbQmDC12lQtWKFdsKAAB/mh6P+GAvkSui7iQP8sW/Hx0PLzcZcovrMHfdATS0cXNGInP4b0dryy3DguDuyl/cbAmDC3UhCAL+d/MJNKq1SIz0wcPc8dmmJUYOxKePTYTSXY780nrMX8+WF6L+0ut/7iaaPTpU5Gro1xhcqIstRy9hx+kquEpd8NqcUZByeWubNzJciU8fnwCluxyHS+vxyEe5aNEwvBD1VV5JHSob2uDlJsPUIZxNZGsYXMiopkmNl749BQBYdHMs4oLYRWQvEkKVhm4jhQwHi2rx+MY8TpUm6qP/HrsEAJieEAyFjN1EtobBhYxe/u8p1DYbZhEtuGGw2OWQiUaF++CjR8ZjgKsUe8/VYMHHh7hIHZGJtDo9th03dBPdNoqziWwRgwsBAPadr8Y3Ry7BRQK8dvcouMr4aNijpKiBWD9/HNzlUmQVXsGiT/PRzu0BiHrtQFEtqps0GOgh56JzNoqfToR2nR4vfnMSAPDghCiMCvcRtyDqlwkxfvhwXrJxhd0/fJ7PvY2Ieunbo4ZuohkjQrg3kY3iT4Xw0d5inK1qgu8AVyzlBooOYXKsPz6YmwRXqQu2Ha/Ec5uOc1dpouvQaPXY3rGp4mx2E9ksBhcnd7mhDSt/PAMAeG7GUCg95CJXROZyY3wgVj2QCKmLBF8duoi/bS2AIDC8EPVk77lq1Le0w99TgQkxfmKXQz1gcHFyf99agGaNDomRPrg7KVzscsjM0hKC8frdowAA6/cWYdVP50SuiMh2benoJvrNyGAuBWHDGFyc2MGiWmw5egkSCfDX20fAhX9RHdJdY8Pxwm3DAQBvZp7BxpxicQsiskHNai22nzB0E92eGCZyNXQtDC5OSq8X8PethjVb7h8fiRFhSpErIkt6ZEo0/jAtDgDwwjcn8XV+ucgVEdmW709WorVdh2j/AUiM8BG7HLoGBhcn9e2xSzh6UQVPhQzptw4RuxyygqduicP8SYMAAM98eRQ7Ci6LWxCRDdncEebvGBMGiYStz7aMwcUJtbXr8Nr2QgDAkzcOhr+nQuSKyBokEgleuG047kwMg04v4P99chgHLtSIXRaR6CpVbdh7rhoAcCe7iWweg4sT+te+YpTXtyJE6YZHuImiU3FxkeC1u0fhlmGBUGv1eOxfeThRrhK7LCJRfXOkHHoBSI4aiEg/D7HLoetgcHEytc0arNppmFmyNC2e27U7IbnUBaseGIsJ0b5oVGvx8Ee5KKttEbssItF0dhPdNZYzK+0Bg4uTeXfHWTS2aTE8xJtNok7MTS7FP+clY2iwF640qjFv/UHUNmvELovI6k5dasDpyka4Sl3wm5FcdM4eMLg4kYt1LfjkQAkAYPmsYZz+7OS83eT46OHxCPNxx4XqZjz6r1y0argpIzmXzfkXAQDThgVyAU47weDiRFb9dA7tOgGTBvthShw3DyMgWOmGfz0yDkp3OfJL67H4s8Pc14ichlanx9dHDIvOsQXafjC4OIni6mZ8ecjwm8UzaZz+TD+LDfTCuo5NGX8sqMKfvznBrQHIKew9X4MrjWoM9JDjxvhAscuhXmJwcRLv7jgLnV7ATfEBSIryFbscsjHJg3zxzn2JcJEAnx0swzs7zopdEpHFZXT8Mjd7dChcZfw4tBf8STmBs5cbsfmIYdR8+q3c/Zm6N2NEMF66fQQAYOWPZ/HZwVKRKyKynPoWjXEn6HuSIkSuhkzB4OIEVv54FoIATE8IwshwLu1PPZs7MQqLbooFADy/+Th+PMXVdckxfZ1fDo1Wj2Eh3hgR5i12OWQCBhcHV1DRgK3HKyCRAE9zaX/qhWfShuCepHDoBWDRZ4dxuLRO7JKIzEoQBHyRZ+gmujc5nEv82xkGFwfXudjcrJEhGBrM3yro+iQSCf7vrpG4MT4Abe2G1XWLq5vFLovIbE6UN6CgogGuMhfcwdlEdsfk4LJr1y7Mnj0boaGhkEgk+Prrr695flZWFiQSyVWv06dP97Vm6qULV5qw7XgFABib/4l6Qy51wfsPjsXIMCVqmzV4+KNcLlBHDuOLPMP4rRkJwfDxcBW5GjKVycGlubkZo0ePxqpVq0y6rrCwEBUVFcZXXFycqd+aTLQ66zwEAbhlWCCGhbC1hUzj4SrDuvnJCPNxR1F1M36/MQ9t7Vygjuxbq0aHb/INa7fcO46Dcu2RzNQLZs6ciZkzZ5r8jQIDA+Hj49Orc9VqNdRqtfHrhoYGk7+fs7tY12Lcf2MhW1uojwK93LDh4XGYs3of8krq8MyXR/GP+xK56jLZre9OVKBRrUWErztSYvzELof6wGpjXBITExESEoJp06Zh586d1zx3xYoVUCqVxldEBFOxqdbuugCtXsDkWD8kRg4UuxyyY0OCvPDB75Igl0qw9VgFXvu+UOySiPrsi9wyAMBvkyIYwO2UxYNLSEgI1q5di4yMDGzatAnx8fGYNm0adu3a1eM1y5Ytg0qlMr7KysosXaZDqWpsw+cdfznZ2kLmMCnWH6/cNQoAsCb7vHHPKyJ7cv5KEw4U1cJFAtydzJ2g7ZXJXUWmio+PR3z8z4uepaSkoKysDG+88QamTp3a7TUKhQIKhcLSpTmsdXuKoNHqMTbSh02hZDZzksJxsa4Vb/94Bn/++gRCle64aSiXSSf78cl+w6Dcm4cGIkTpLnI11FeiTIeeOHEizp7lkuKW0NjWjk87/nIuvCmW6xOQWS2ZFou7O9Z4WfjpYZwoV4ldElGvtGp0+OqQoSX6wYlRIldD/SFKcMnPz0dISIgY39rhfZFbhka1FoMDBuAmbhpGZiaRSPB/d47E5Fg/tGh0eOSjXFyqbxW7LKLr+vboJTS0GQbl3hAXIHY51A8mB5empiYcOXIER44cAQAUFRXhyJEjKC01/Ja/bNkyPPTQQ8bzV65cia+//hpnz57FyZMnsWzZMmRkZGDRokXmeQdkpNXpsWFvMQDgsdQYDjwji3CVueD9B5MwJMgTVY1qPLwhFw1t7WKXRXRNH3eMy3pwQhT/bbRzJgeXvLw8JCYmIjExEQCQnp6OxMREvPDCCwCAiooKY4gBAI1Gg6VLl2LUqFFITU3Fnj17sHXrVtx1111megvU6bsTlSivb4XfAFfcydUgyYKU7nJseHg8ArwUKLzciP/38WG06/Ril0XUraNl9Th2UQVXqQvuSeKgXHsnEQRBELuI62loaIBSqYRKpYK3NxdS644gCLjjvb04elGFp26Jw1O3cF8isrwT5Sr89oMctGh0uCcpHK/dPYrjqsjm/PHLo/jy0EXcMSYUK+9LFLscp2KJz2/uVeQgcovrcPSiCgqZC+Zy4BlZyYgwJf5xfyJcJMCXhy5i1U/nxC6JqAtVSzu+PWZYKXduCv9tdAQMLg7in7svAADuGhsOP09OJSfrmTYsCC/9TwIA4M3MM9icf1Hkioh+9uWhMrS16zE02AtjuRinQ2BwcQDF1c34seAyAODRKdEiV0POaG7KIPx+agwA4E9fHUPO+RqRKyICdHoB/8opBmBobWE3pmNgcHEA/95fAkEAbooPQGygp9jlkJN6bsZQ/GZkCNp1Ap74dx7OVTWKXRI5ucxTl1FW2wofDznuSuSgXEfB4GLnWjRa/CfPsKjSQ5MGiVsMOTUXFwne/O1oJEUNREObFvPW56KqsU3sssiJrd9bBAB4YHwk3F2lIldD5sLgYue+zr+ExjYtovw8uKgSic5NLsU/H0rGID8PlNe34rF/5aFFoxW7LHJCJ8pVOFhUC5mLBA+lDBK7HDIjBhc7JggCNnb2307kokpkG3wHuOKjh8djoIccxy6qsOSzfOj0Nr/qAjmY9XsMrS2/GRWCYKWbyNWQOTG42LGDRbU4XdkId7kU9yRFiF0OkdEg/wH4cF4yXGUu+LGgCi9/exJ2sGQUOYiqhjbjFOhHJnPCgqNhcLFjG3MMS1jfkRgGpYdc5GqIukqK8sXKe8cAAP6VU4J1Hb8BE1nav/eXoF0nIDlqIEZH+IhdDpkZg4udqlS1YfvJSgDAQ1xUiWzUrJEhWD5rKADg79sK8N3xCpErIkfXqtHhkwOGbWce4fIQDonBxU59eqAEOr2A8dG+GBbCbRDIdj2eGoO5E6MgCMBTXxzBoZI6sUsiB/afvDLUNmsQPtAdacODxC6HLIDBxQ5pdXp80TkFmq0tZOMkEglenD0c04YGQq3V4/GNeSiubha7LHJA7To91u4yrCL+xNQYyKT8iHNE/KnaoZ2FV3C5QQ2/Aa5IGx4sdjlE1yWTuuAfDyRiZJgStc0aPPxRLuqaNWKXRQ5m67EKlNe3wm+AK+5J5oQFR8XgYoc+P2jov52TFA5XGX+EZB88XGVYNy8ZYT7uKKpuxuMb89DWrhO7LHIQgiBgddZ5AMDDkwfBTc4F5xwVP/XsTIWqFTsLqwAA947jbxRkXwK93bDh4XHwcpMhr6QOz3x5FHqu8UJmsLOwCoWXGzHAVYq5EweJXQ5ZEIOLnfky7yL0AjB+kC8GB3BfIrI/Q4K88MHcJMilEmw9VoHXvi8UuyRyAGuyDGNbHpwYxeUhHByDix3R6wV8kWsYlHvfeLa2kP2aNNgfr84ZBQBYk30enxwoEbkismeHSmpxsLgWrlIXPMop0A6PwcWO7D5XjfL6Vni7yTBrZIjY5RD1y11jw5F+6xAAwJ+/PoGdp6tErojs1Ts7zgEA7kwMQ5A3l/d3dAwudqRzUO6diWEceEYOYfHNsbg7KRx6AVj46WGcKFeJXRLZmUMlddh15gpkLhIsvClW7HLIChhc7MSVRjUyT10GANw3PlLkaojMQyKRYMVdIzEl1h8tGh0e+SgX5fWtYpdFdmTlj2cAAHPGhiPSz0PkasgaGFzsRMbhi9DqBYyO8OFKueRQ5FIXvP+7sRga7IWqRjUe2ZCLhrZ2scsiO3CopBa7z1aztcXJMLjYAUEQ8NWhiwCA+zgFmhyQt5sc6+ePQ6CXAoWXG/H7jXlQa7nGC13b25lnAbC1xdkwuNiB4+UqnKtqgkLmgt+M4qBcckyhPu7Y8PA4eCpk2H+hFulfHIWOa7xQD3KLa7HnnKG1ZdHNbG1xJgwudiCjo7UlLSEY3m5cn4AcV0KoEms713g5XoGXvz0JQWB4oat1jm25OykcEb5sbXEmDC42TqPVY8vRSwCAOWPDRK6GyPImxfrjrd+OgUQC/CunBO93LONO1GnP2WrsPVcDuZRjW5wRg4uNyyqsQl1LOwK8FJgS6y92OURWMXt0KF64bTgA4PXvC/GfjoUXifR6Aa9uPw0AeHBCFFtbnBCDi43LOGzoJrozMYxbtJNTeXhyNJ68cTAAYNnm49hRcFnkisgWbD1egePlKngqZFjMsS1OiZ+ENqyuWYOfOlYTvYvdROSE/jQ9HncnhUOnF7Dw08M4VFIndkkkIo1Wjzd+MOxt9fupMfDzVIhcEYmBwcWG/ffYJbTrBAwP8cbQYK7dQs6nc4G6m+ID0NauxyMf5eLs5UaxyyKRfJ5bipKaFvh7KrgnkRNjcLFhXx0uBwDMSQoXuRIi8cilLnjvwbEYE+EDVWs7Hlp/EBUqrq7rbJrUWry7w7Buyx9uicMAhUzkikgsDC426lxVE46W1UPqIsH/jA4VuxwiUXm4yrB+/jjEBAxAhaoND607iLpmjdhlkRWtyTqP6iYNov0HcCFOJ8fgYqM25xsG5d44JAABXuzHJfId4IqNj4xHkLcCZ6uaMH/DQTRyawCnUFrTgrW7LwAAnp0xFHJOVHBq/OnbIEEQjGu33JHIQblEncIHeuDjRydgoIccRy+q8Oi/8tDWzq0BHN3ft52CRqvH5Fg/TE8IErscEhmDiw3KL6tHWW0rPFyluGUY/5IS/VJckBc2PjIBXgoZDhbV4smPD0Gj1YtdFlnI3nPV+P7kZUhdJHhxdgIkEonYJZHIGFxs0LcdrS23Dg+Cu6tU5GqIbM/IcCXWzR8HN7kLdhZewdP/OcJ9jRxQu06Pl749CQCYOzEKQ4K8RK6IbAGDi43R6QX891gFAHBQLtE1jI/2xQdzkw37Gh2rwPJNx7mvkYP5eH8JzlxuwkAPOZ6+ZYjY5ZCNYHCxMQeKanClUQ2luxypcQFil0Nk024YEoB370uEiwT4Iq8Mf9tawPDiIKoa2vBWpmEjxWfS4qH04AazZMDgYmM6u4lmjgiGq4w/HqLrmTkyBK/OGQUAWLenCCt/PCtyRWQOL/33FBrbtBgZpsT94yPFLodsCD8ZbYhGq8d3JyoBGDaZI6LeuSc5An+ZbdiU8Z0dZ/HeznMiV0T98dPpy9h6rAJSF8PKyVIXDsilnzG42JA9566gvqUd/p4KTIzxE7scIrsyf3I0np0xFIBhR+n3sxhe7FGzWos/f20YkPvI5EEYEaYUuSKyNQwuNmTLEUM30W2jQvgbBlEfPHnjYPxxejwA4LXthfgg+7zIFZGp3so8g/L6VoT5uOPpWzkgl67G4GIjWjU6ZJ66DIDdRET9sfCmWDzT8YG34rvT+OeuCyJXRL117GI9NuwtAgD87c4R8HDlfkR0NQYXG/HT6So0a3QI83HH2EgfscshsmuLp8UZp8/+fVsBPtzN8GLr2tp1SP/PUegFwy9vN8UHil0S2SgGFxux9XhHN9HoEK4MSWQGf7glDkumxQEA/ra1AOv3FIlcEV3L698X4lxVEwK8FHjpfxLELodsGIOLDWjV6LDz9BUAwG0j2U1EZC5P3xKHxTfHAgBe/u8prM7imBdbtO98NdZ1BMtX54yE7wBXkSsiW8bgYgOyz1ShtV2H8IHuGBHmLXY5RA5DIpEg/dYhxpaXV7efxps/FHKROhvS2NaOP355DABw//gI3DyU+7PRtTG42IDOtVtmJASzm4jIzDrDy3MzDVOl//HTOa6wa0Ne+vYUyutbEeHrjud/M1zscsgOMLiITK3V4aeCKgCGFUCJyDIW3DAYL99uGDuxbk8Rlm8+AT03ZhTVpsMX8dWhi5BIgDfvGQNPBWcR0fUxuIhs77lqNKq1CPJWIDHCR+xyiBzaQymD8Prdo+AiAT47WIpnvjwKrU4vdllO6VxVE/736xMAgD9Mi8P4aF+RKyJ7weAism3Hf+4mcuGic0QWd09yBN65LxEyFwk255fjiX8fQqtGJ3ZZTqVVo8PCTw6jRaPDpMF+WHxznNglkR0xObjs2rULs2fPRmhoKCQSCb7++uvrXpOdnY2kpCS4ubkhJiYGa9as6UutDqddpzcuOjdjBLuJiKxl9uhQfDA3CQqZC3acrsIDH+5HbbNG7LKcgiAIeHHLCRReboS/pwIr7xvDlcLJJCYHl+bmZowePRqrVq3q1flFRUWYNWsWUlNTkZ+fj+XLl2PJkiXIyMgwuVhHs/9CDVSt7fAb4MpmUiIrmzYsCJ8+PgE+HnLkl9bj7jX7UFbbInZZDu/j/SX4T55hXMvKe8cg0MtN7JLIzpg8EmrmzJmYOXNmr89fs2YNIiMjsXLlSgDAsGHDkJeXhzfeeANz5swx9ds7lM7ZRGkJwfyNg0gESVG++GpBCuatz8WFK82Ys3ofPnp4PIaHclkCS8g5X4OXvj0FAHh2xlBMifMXuSKyRxYf45KTk4O0tLQux6ZPn468vDy0t7d3e41arUZDQ0OXl6PR6QX8cNIQXGaOCBa5GiLnFRvohYwnJ2FosBeqGtW494Mc7D57ReyyHE5ZbQsWfnoYWr2A/xkdiiemxohdEtkpiweXyspKBAV1XVAoKCgIWq0W1dXV3V6zYsUKKJVK4ysiIsLSZVpdXnEtqps0ULrLkTLYT+xyiJxasNINXzyRggnRvmhUazF/Qy425hSLXZbDaGhrx+Mb81DbrMGIMG+8OmcU16yiPrPKrKJfP6CdCz/19OAuW7YMKpXK+CorK7N4jdbW2U106/AgyKWc3EUkNqW7HBsfHY+7xoZBpxfwwjcn8eevT3C6dD9ptHo8sfEQTlc2IsBLgQ/mJsPdVSp2WWTHLL7aT3BwMCorK7scq6qqgkwmg59f9y0NCoUCCoXC0qWJRhAE/FhgmE00PYHdRES2QiGT4s17RmNIkBde3X4a/95fgqLqZrz34Fgo3eVil2d39HoBf/zqKHIu1GCAqxQb5o9DmI+72GWRnbP4r/opKSnIzMzscuyHH35AcnIy5HLn/IfgdGUjLta1wk3ugimxHJxGZEskEgkW3DAYa36XBA9XKfacq8ad7+/FuaomsUuzK4IgYMV3BfjmyCXIXCRY/bskjAhTil0WOQCTg0tTUxOOHDmCI0eOADBMdz5y5AhKS0sBGLp5HnroIeP5CxYsQElJCdLT01FQUID169dj3bp1WLp0qXnegR36sWPtlimxAWwyJbJR0xOC8eWCFIQq3XDhSjNuX7UHW49ViF2W3Xgr8wz+uduw4/Mrc0Zh6pAAkSsiR2FycMnLy0NiYiISExMBAOnp6UhMTMQLL7wAAKioqDCGGACIjo7Gtm3bkJWVhTFjxuCvf/0r3n33XaeeCp3Z0U106/BAkSshomtJCFXim0VTkBLjh2aNDgs/PYyXvz2Fdo57uaZ3d5zFP346BwB4cfZw3J0ULnJF5Egkgh1skdrQ0AClUgmVSgVvb/teX6FS1YaJK3ZAIgEOLr8FAV6OO5aHyFFodXq8mXkGq7POAwCSowZi1QNjEazk4mm/9t7Oc3j9+0IAwPJZQ/H7qYNFrojEZInPb05nsbIdpw2tLYkRPgwtRHZCJnXBszOGYu3cJHi5yZBXUoeZ7+zC9ycrr3+xk+gc09IZWpamDWFoIYtgcLGyzr2JbhkedJ0zicjWpCUE49tFU5AQ6o26lnY88e9DeC7jGJrVWrFLE5VOL2DZpuP4IPsCAENLyyJunEgWwuBiRc1qLfadqwEA3DqMwYXIHg3yH4DN/28ynrghBhIJ8HluGX7z7m4cKasXuzRRNKu1WPDxIXyeWwYXCfDqnJFsaSGLYnCxol1nrkCj02OQnwdiAz3FLoeI+shV5oJlM4fhk8cmIETphuKaFsxZvQ+vbj+Ntnad2OVZzcU6w/vOPHUZrlIXvPfAWNw7LlLsssjBMbhYUedsoluGBXG5ayIHMGmwP7b/YSpmjw6FTi9gddZ5zFi5C/vOd7+diSM5WFSL21ftxenKRvh7KvD5ExMxc2SI2GWRE2BwsRKtTo+fTlcBMCzzT0SOQekhxz/uT8TauUkI8laguKYFD/zzAJ796hhqmtRil2d2Or2Af+w4i/vW5qCmWYPhId74ZtFkjI0cKHZp5CQYXKzkUEkd6lva4eMhR1IU/4ITOZq0hGBkpt+A3000dJV8kVeGG9/Iwoe7L0CjdYx1Xy43tGHuugN4M/MM9AJwZ2IYvnoyhcv4k1UxuFhJ595EN8cHQsZNFYkckrebHH+7YyS+XJCC4SHeaGzT4m9bCzBj5S78eOoy7GDZrG4JgoDPD5bilreyse98DdzlUrxxz2i8fe8YeLhafMs7oi74xFnJzsIrAIBpnE1E5PDGDfLFt4un4KtDZXj9+0JcqG7GYxvzkBjpg6dvGYLUOH+7Ged2/koT/nfzCeRcMMyIHBWuxNv3jsHgAE4wIHEwuFhBWW0LzlU1QeoiwZQ4bqpI5AykLhLcOy4Ss0aG4P2s89iwtwj5pfV4aP1BJEcNxMKbY3FDXABcXGwzwNQ0qfHOjrP45EApdHoBbnIXLE2Lx8OToyG10ZrJOTC4WEFWoWFQblLUQCjdnXNHbCJn5eUmx7MzhuLhyYOwJusCPjlQgrySOjy8IReDAwbg4cnRmDM23GY2XK1pUuOjfcXYsLcYTR0L600bGogXZycg0s9D5OqIuFeRVTzyUS5+Ol2FZ2cMxZM3cmEmImd2uaEN/9x1AV/klqGxIxh4u8kwe3Qo5iSFIzHCR5RupHNVTdiYU4z/5JWhrd0wmHhkmBLLZw1DymA/q9dDjsESn98MLhbW1q7DmJd/QFu7HtufSsXQYPuqn4gso7GtHV/mXcSGfUUoq201Ho8JGIBZI0IwbVggRof7WLQrqa5Zg+0nK/FlXhkOl9Ybj48KV2LBDYMxIyHYZruyyD4wuNhhcMkqrML8DbkIUbph33M3282APCKyDp1eQM75GmQcvojvTlQYWzsAwN9TgalD/DFukC/GRg5EXKBnv4JEW7sOJy+pcKCoFj8VVOFwaR30HZ8AUhcJbhwSgEenRCNlsB//rSKzsMTnN8e4WFhWx2yiG+MD+Q8BEV2lc9D+lDh//PWOEfjhZCV2FFQh+8wVVDepselwOTYdLgcAeLnJEBfoiZgAT0T7D0CglwL+ngoMHOAKuVQCqYsEEkjQpG5HQ6sWdS0alNa2oKTGMEHgdGUD2nVdf1cdGuyFOxPDcGdiGAK93cS4BUQmYXCxIEEQjKvl3hQfIHI1RGTrPBUy3DU2HHeNDYdGq8eBohocuFCLQyV1OHqxHo1tWhwure/SrWMqf09XjIkYiBviA3Dz0EAuHkd2h8HFgoqqm1Fa2wK5VILJsZwGTUS95ypzQWpcAFLjDL/0aHV6nK1qwvkrTbhwpRnFNc2obtKgpkmN+pZ2tOv00AsCBAHwdJPBy00Gpbsc4T4eiPL3wCC/ARgZpkT4QHe2/pJdY3CxoM5F5yZE+2GAgreaiPpOJnXBsBBvDAuxr3F+RObGtectqHP9lhvZTURERGQWDC4W0qzW4sCFWgDATUMDRa6GiIjIMTC4WMi+8zXQ6PSI8vNAjP8AscshIiJyCAwuFmLsJhoSwIFwREREZsLgYiG7z1YDAKYO4fgWIiIic2FwsYCSmp+nQU+M4R4fRERE5sLgYgG7OlpbxkYO5DRoIiIiM2JwsYDdZwzrt7CbiIiIyLwYXMysXadHzvkaAMDUOAYXIiIic2JwMbOjZfVoVGsx0EOOhFCucElERGRODC5m1jm+ZUpcQL+2nyciIqKrMbiY2e6zhvEtqXHcVJGIiMjcGFzMSNXSjqNl9QAYXIiIiCyBwcWM9p6vhl4A4gI9EaJ0F7scIiIih8PgYkY/dxNxNhEREZElMLiYiSAI2HXGMDA3dQi7iYiIiCyBwcVMiqqbUV7fClepCyZE+4pdDhERkUNicDGTzk0VkwcNhIcrl/knIiKyBAYXM9l7rnP9FnYTERERWQqDixno9AL2XzAs8z95MIMLERGRpTC4mMHJSyo0tGnh5SbjMv9EREQWxOBiBvs6NlWcEO0HmZS3lIiIyFL4KWsGneNbJsf6iVwJERGRY2Nw6SeNVo/c4loAwCSObyEiIrIoBpd+yi+tQ1u7Hv6erhgS5Cl2OURERA6NwaWfOse3pAz2h0QiEbkaIiIix8bg0k/7zhvGt0wazPEtRERElsbg0g8tGi3yS+sBcP0WIiIia2Bw6YeDRbXQ6gWE+bgjwtdd7HKIiIgcHoNLP+R0jG+ZHOvH8S1ERERWwODSD50DczkNmoiIyDoYXPqovkWDE5dUADgwl4iIyFoYXPpo/4VaCAIQG+iJQG83scshIiJyCn0KLu+//z6io6Ph5uaGpKQk7N69u8dzs7KyIJFIrnqdPn26z0Xbgs7doFNi2NpCRERkLSYHly+++AJPPfUUnn/+eeTn5yM1NRUzZ85EaWnpNa8rLCxERUWF8RUXF9fnom1BZ3CZyOBCRERkNSYHl7feeguPPvooHnvsMQwbNgwrV65EREQEVq9efc3rAgMDERwcbHxJpdI+Fy22+hYNCi83AgDGR/uKXA0REZHzMCm4aDQaHDp0CGlpaV2Op6WlYd++fde8NjExESEhIZg2bRp27tx5zXPVajUaGhq6vGzJwSLD+JbBAQMQ4KUQuxwiIiKnYVJwqa6uhk6nQ1BQUJfjQUFBqKys7PaakJAQrF27FhkZGdi0aRPi4+Mxbdo07Nq1q8fvs2LFCiiVSuMrIiLClDIt7kCRYTfoCewmIiIisipZXy769WJrgiD0uABbfHw84uPjjV+npKSgrKwMb7zxBqZOndrtNcuWLUN6errx64aGBpsKLweKDONbJrCbiIiIyKpManHx9/eHVCq9qnWlqqrqqlaYa5k4cSLOnj3b4/9XKBTw9vbu8rIVDW3tOHXJ0HU1IZotLkRERNZkUnBxdXVFUlISMjMzuxzPzMzEpEmTev3n5OfnIyQkxJRvbTMOFddBLwBRfh4IVnL9FiIiImsyuasoPT0dc+fORXJyMlJSUrB27VqUlpZiwYIFAAzdPOXl5di4cSMAYOXKlRg0aBASEhKg0Wjw8ccfIyMjAxkZGeZ9J1ayn91EREREojE5uNx7772oqanByy+/jIqKCowYMQLbtm1DVFQUAKCioqLLmi4ajQZLly5FeXk53N3dkZCQgK1bt2LWrFnmexdWdOBCx8BcdhMRERFZnUQQBEHsIq6noaEBSqUSKpVK1PEuzWotRr30A3R6AXuevQnhAz1Eq4WIiMjWWeLzm3sVmeBQSR10egFhPu4MLURERCJgcDEBp0ETERGJi8HFBMbxLTEMLkRERGJgcOmlVo0ORy/WA+DAXCIiIrEwuPRSflkd2nUCgrwViPLj+BYiIiIxMLj00i+nQfe0vQERERFZFoNLLxkH5nJ8CxERkWgYXHpBo9Ujv7QeAGcUERERiYnBpRdOXFJBrdVjoIccgwM8xS6HiIjIaTG49MKh4joAQFKUL8e3EBERiYjBpRdyiw0Dc8cNGihyJURERM6NweU6BEHAoRJDi0sygwsREZGoGFyuo6i6GTXNGrjKXDAiTCl2OURERE6NweU68jpaW0aHK6GQSUWuhoiIyLkxuFxHXsf4lqQoToMmIiISG4PLdXS2uHBgLhERkfgYXK6hpkmNC1eaAQBJUQwuREREYmNwuYbO2URxgZ7w8XAVuRoiIiJicLkGToMmIiKyLQwu15DLgblEREQ2hcGlB23tOpwobwDAgblERES2gsGlB8fLVdDo9PD3VCDS10PscoiIiAgMLj365f5E3FiRiIjINjC49ODnHaHZTURERGQrGFy6odcLv1h4jgNziYiIbAWDSzfOX2mCqrUd7nIphod6i10OERERdWBw6YZxY8UIJeRS3iIiIiJbwU/lbuSXcnwLERGRLWJw6cbh0noAQGIEgwsREZEtYXD5FVVrO85VNQEAEiN9xC2GiIiIumBw+ZWjZfUAgCg/D/h5KsQthoiIiLpgcPmVfGM3kY+odRAREdHVGFx+5XDHwNzESI5vISIisjUMLr+g1ws40tFVNJbBhYiIyOYwuPxCUU0zVK3tUMhcMDTES+xyiIiI6FcYXH7hcMfCc6PCufAcERGRLeKn8y/ks5uIiIjIpjG4/IJxRhHXbyEiIrJJDC4dmtRaFFY2AOCMIiIiIlvF4NLh2MV66AUgVOmGIG83scshIiKibjC4dDB2E3FjRSIiIpvF4NKBK+YSERHZPgYXAIIgIJ8r5hIREdk8BhcAZbWtqGnWQC6VICHUW+xyiIiIqAcMLgDyywytLQmhSrjJpSJXQ0RERD1hcMHPK+Zy/RYiIiLbxuCCn1fM5fgWIiIi2+b0waWtXYeCio6F5zijiIiIyKY5fXApqGhAu06A7wBXhA90F7scIiIiuganDy7HLqoAAKPDlZBIJCJXQ0RERNfi9MHl6MV6AMCocB9R6yAiIqLr61Nwef/99xEdHQ03NzckJSVh9+7d1zw/OzsbSUlJcHNzQ0xMDNasWdOnYi3haMfA3DEc30JERGTzTA4uX3zxBZ566ik8//zzyM/PR2pqKmbOnInS0tJuzy8qKsKsWbOQmpqK/Px8LF++HEuWLEFGRka/i++vxrZ2XKhuBgCMCleKXA0RERFdj0QQBMGUCyZMmICxY8di9erVxmPDhg3DHXfcgRUrVlx1/rPPPostW7agoKDAeGzBggU4evQocnJyuv0earUaarXa+LVKpUJkZCTKysrg7W2+lW0PXKjBo//KQ4jSDZnpN5jtzyUiIiKgoaEBERERqK+vh1JppgYCwQRqtVqQSqXCpk2buhxfsmSJMHXq1G6vSU1NFZYsWdLl2KZNmwSZTCZoNJpur3nxxRcFAHzxxRdffPHFlwO8zp8/b0rcuCYZTFBdXQ2dToegoKAux4OCglBZWdntNZWVld2er9VqUV1djZCQkKuuWbZsGdLT041f19fXIyoqCqWlpeZLbA6iM82auzXKEfDedI/3pWe8Nz3jvekZ703POntMfH19zfZnmhRcOv162rAgCNecStzd+d0d76RQKKBQKK46rlQq+VD0wNvbm/emB7w33eN96RnvTc94b3rGe9MzFxfzTWI26U/y9/eHVCq9qnWlqqrqqlaVTsHBwd2eL5PJ4OfnZ2K5RERE5MxMCi6urq5ISkpCZmZml+OZmZmYNGlSt9ekpKRcdf4PP/yA5ORkyOVyE8slIiIiZ2Zy2016ejo+/PBDrF+/HgUFBXj66adRWlqKBQsWADCMT3nooYeM5y9YsAAlJSVIT09HQUEB1q9fj3Xr1mHp0qW9/p4KhQIvvvhit91Hzo73pme8N93jfekZ703PeG96xnvTM0vcG5OnQwOGBehee+01VFRUYMSIEXj77bcxdepUAMD8+fNRXFyMrKws4/nZ2dl4+umncfLkSYSGhuLZZ581Bh0iIiKi3upTcCEiIiISg9PvVURERET2g8GFiIiI7AaDCxEREdkNBhciIiKyGzYZXP7+979j0qRJ8PDwgI+PT6+uEQQBf/nLXxAaGgp3d3fceOONOHnypGULFUFdXR3mzp0LpVIJpVKJuXPnor6+/prXzJ8/HxKJpMtr4sSJ1inYgt5//31ER0fDzc0NSUlJ2L179zXPz87ORlJSEtzc3BATE4M1a9ZYqVLrM+XeZGVlXfV8SCQSnD592ooVW8euXbswe/ZshIaGQiKR4Ouvv77uNc7y3Jh6b5zluVmxYgXGjRsHLy8vBAYG4o477kBhYeF1r3OG56Yv98Ycz41NBheNRoN77rkHTz75ZK+vee211/DWW29h1apVyM3NRXBwMG699VY0NjZasFLre+CBB3DkyBFs374d27dvx5EjRzB37tzrXjdjxgxUVFQYX9u2bbNCtZbzxRdf4KmnnsLzzz+P/Px8pKamYubMmSgtLe32/KKiIsyaNQupqanIz8/H8uXLsWTJEmRkZFi5cssz9d50Kiws7PKMxMXFWali62lubsbo0aOxatWqXp3vTM+Nqfemk6M/N9nZ2Vi4cCH279+PzMxMaLVapKWlobm5ucdrnOW56cu96dSv58Zs2zVawIYNGwSlUnnd8/R6vRAcHCy88sorxmNtbW2CUqkU1qxZY8EKrevUqVMCAGH//v3GYzk5OQIA4fTp0z1eN2/ePOH222+3QoXWM378eGHBggVdjg0dOlR47rnnuj3/T3/6kzB06NAux5544glh4sSJFqtRLKbem507dwoAhLq6OitUZzsACJs3b77mOc703PxSb+6Nsz43VVVVAgAhOzu7x3Oc9bnpzb0xx3Njky0upioqKkJlZSXS0tKMxxQKBW644Qbs27dPxMrMKycnB0qlEhMmTDAemzhxIpRK5XXfZ1ZWFgIDAzFkyBA8/vjjqKqqsnS5FqPRaHDo0KEuP28ASEtL6/E+5OTkXHX+9OnTkZeXh/b2dovVam19uTedEhMTERISgmnTpmHnzp2WLNNuOMtz0x/O9tyoVCoAuOZux8763PTm3nTqz3PjEMGlcxPHX2/0GBQUdNUGj/assrISgYGBVx0PDAy85vucOXMmPvnkE/z000948803kZubi5tvvhlqtdqS5VpMdXU1dDqdST/vysrKbs/XarWorq62WK3W1pd7ExISgrVr1yIjIwObNm1CfHw8pk2bhl27dlmjZJvmLM9NXzjjcyMIAtLT0zFlyhSMGDGix/Oc8bnp7b0xx3MjM0fBvfGXv/wFL7300jXPyc3NRXJycp+/h0Qi6fK1IAhXHbNFvb03wNXvEbj++7z33nuN/z1ixAgkJycjKioKW7duxV133dXHqsVn6s+7u/O7O+4ITLk38fHxiI+PN36dkpKCsrIyvPHGG8atPJyZMz03pnDG52bRokU4duwY9uzZc91zne256e29McdzY7XgsmjRItx3333XPGfQoEF9+rODg4MBGFJuSEiI8XhVVdVVqdcW9fbeHDt2DJcvX77q/125csWk9xkSEoKoqCicPXvW5Fptgb+/P6RS6VUtCNf6eQcHB3d7vkwmg5+fn8Vqtba+3JvuTJw4ER9//LG5y7M7zvLcmIsjPzeLFy/Gli1bsGvXLoSHh1/zXGd7bky5N90x9bmxWnDx9/eHv7+/Rf7s6OhoBAcHIzMzE4mJiQAMff3Z2dl49dVXLfI9zam39yYlJQUqlQoHDx7E+PHjAQAHDhyASqXCpEmTev39ampqUFZW1iXk2RNXV1ckJSUhMzMTd955p/F4ZmYmbr/99m6vSUlJwbffftvl2A8//IDk5GTI5XKL1mtNfbk33cnPz7fb58OcnOW5MRdHfG4EQcDixYuxefNmZGVlITo6+rrXOMtz05d70x2Tn5s+D+u1oJKSEiE/P1946aWXBE9PTyE/P1/Iz88XGhsbjefEx8cLmzZtMn79yiuvCEqlUti0aZNw/Phx4f777xdCQkKEhoYGMd6CxcyYMUMYNWqUkJOTI+Tk5AgjR44Ubrvtti7n/PLeNDY2Cs8884ywb98+oaioSNi5c6eQkpIihIWF2fW9+fzzzwW5XC6sW7dOOHXqlPDUU08JAwYMEIqLiwVBEITnnntOmDt3rvH8CxcuCB4eHsLTTz8tnDp1Sli3bp0gl8uFr776Sqy3YDGm3pu3335b2Lx5s3DmzBnhxIkTwnPPPScAEDIyMsR6CxbT2Nho/PcEgPDWW28J+fn5QklJiSAIzv3cmHpvnOW5efLJJwWlUilkZWUJFRUVxldLS4vxHGd9bvpyb8zx3NhkcJk3b54A4KrXzp07jecAEDZs2GD8Wq/XCy+++KIQHBwsKBQKYerUqcLx48etX7yF1dTUCA8++KDg5eUleHl5CQ8++OBV08p+eW9aWlqEtLQ0ISAgQJDL5UJkZKQwb948obS01PrFm9l7770nREVFCa6ursLYsWO7TMGbN2+ecMMNN3Q5PysrS0hMTBRcXV2FQYMGCatXr7ZyxdZjyr159dVXhcGDBwtubm7CwIEDhSlTpghbt24VoWrL65yK+evXvHnzBEFw7ufG1HvjLM9Nd/fk158/zvrc9OXemOO5kXR8cyIiIiKb5xDToYmIiMg5MLgQERGR3WBwISIiIrvB4EJERER2g8GFiIiI7AaDCxEREdkNBhciIiKyGwwuREREZDcYXIiIiMhuMLgQERGR3WBwISIiIrvx/wECZJIio1pptgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-1,2.5,1000)\n",
    "plt.plot(x,f(x))\n",
    "plt.xlim([-1,2.5])\n",
    "plt.ylim([0,3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from plot above that our local minimum is going to be around 1.4 or 1.5 (on the x-axis), but let's pretend that we don't know that, so we set our starting point (arbitrarily, in this case) at $x_0 = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Local minimum occurs at: 1.33\n",
      "Number of steps: 17\n"
     ]
    }
   ],
   "source": [
    "x_old = 0\n",
    "x_new = 2 # The algorithm starts at x=2\n",
    "n_k = 0.1 # step size\n",
    "precision = 0.0001\n",
    "\n",
    "x_list, y_list = [x_new], [f(x_new)]\n",
    "\n",
    "# returns the value of the derivative of our function\n",
    "def f_gradient(x):\n",
    "    return 3*x**2-4*x\n",
    " \n",
    "while abs(x_new - x_old) > precision:\n",
    "    \n",
    "    x_old = x_new\n",
    "    \n",
    "    # Gradient descent step\n",
    "    s_k = -f_gradient(x_old)\n",
    "    \n",
    "    x_new = x_old + n_k * s_k\n",
    "    \n",
    "    x_list.append(x_new)\n",
    "    y_list.append(f(x_new))\n",
    "    \n",
    "print (\"Local minimum occurs at: {:.2f}\".format(x_new))\n",
    "print (\"Number of steps:\", len(x_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figures below show the route that was taken to find the local minimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAEnCAYAAABhZjbYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABspElEQVR4nO3dd1hTZxsG8DvsIaKgLFHAiVvEAS5U3KNWraNaRVvrqNZVraO1am0dre1nh7NVW+uoVdxSq1XAhdaBo3Vb3OAWEJWV9/vjbYIREILAScL9u65cSU7OSZ5zAufNc96lEkIIEBERERERmRAzpQMgIiIiIiLKb0x0iIiIiIjI5DDRISIiIiIik8NEh4iIiIiITA4THSIiIiIiMjlMdIiIiIiIyOQw0SEiIiIiIpPDRIeIiIiIiEwOEx0iIiIiIjI5THQo3506dQrvvPMOKlSoAFtbW9ja2qJSpUoYMmQIjh49WmhxTJs2DSqVSmeZt7c3BgwYUKCfe/DgQUybNg2PHj16pfcZMGAAvL298yUmQ3Hr1i1MmzYNJ06cUDoUIiIiMnFMdChfLV68GP7+/jh8+DBGjRqFbdu2Yfv27Rg9ejT++ecf1K9fH5cvX1Ysvo0bN2LKlCkF+hkHDx7E9OnTXznRMUW3bt3C9OnTmegQGThesOIFq+zk9YLVihUrULp0aSQmJhZMYEbip59+gkqlwpUrV166XlZ/+7m1e/duFCtWDDdv3szT9qbEQukAyHQcOHAA7733Hjp27Ij169fDyspK+1rLli0xfPhwrFu3Dra2ti99nydPnsDOzq5AYvTz8yuQ9yUiMhWLFy/GiBEjUKVKFYwaNQrVq1eHSqXC2bNnsWbNGtSvXx+XLl1ChQoVFIlv48aNKF68eIF+huaC1YABA1CiRIkC/Sxjo7lg5e3tjTp16uRqmydPnmDy5MmYMGECHBwcCjZAEzFo0CC0a9cuT9sGBwejQYMGmDx5Mn7++ed8jsy4sEaH8s3MmTNhbm6OxYsX6yQ5z+vRowc8PDy0zwcMGIBixYrh9OnTaNOmDRwcHBAcHAwA2LVrF7p06QJPT0/Y2NigYsWKGDJkCO7du5fpfbdv3446derA2toaPj4+mDt3bpafn9WVwISEBIwbNw4+Pj6wsrJCmTJlMHr0aCQlJemsp1KpMGLECPzyyy+oWrUq7OzsULt2bWzbtk27zrRp0zB+/HgAgI+PD1QqFVQqFSIiIl567H766SdUqVIF1tbWqFq1KlasWJHleikpKfjss8/g6+sLa2trlC5dGgMHDsTdu3d11tuzZw+aN28OZ2dn2Nraoly5cujevTuePHmiXSc5ORmffvopqlatChsbGzg7O6NFixY4ePCgdh0hBBYsWIA6derA1tYWJUuWxBtvvIF///1X5/OaN2+OGjVq4MiRI2jatCns7OxQvnx5zJ49G2q1GgAQERGB+vXrAwAGDhyoPTbTpk176bEhosKjuWDVvn17HD9+HCNHjkRwcLD2YtX+/fvx22+/5eqCVUHx8/NTLMmivPn5559x//59DBo0SOlQjIanpycCAgLyvP3w4cOxatUqXL9+PR+jMkKCKB+kpaUJW1tbERgYqNd2ISEhwtLSUnh7e4tZs2aJ3bt3iz/++EMIIcTChQvFrFmzxJYtW0RkZKT4+eefRe3atUWVKlVESkqK9j3+/PNPYW5uLpo0aSI2bNgg1q1bJ+rXry/KlSsnXvwT9/LyEiEhIdrnSUlJok6dOqJUqVLi66+/Fn/++af45ptvhKOjo2jZsqVQq9XadQEIb29v0aBBA/Hbb7+JsLAw0bx5c2FhYSEuX74shBDi+vXr4v333xcAxIYNG0RUVJSIiooS8fHx2R6D5cuXCwCiS5cuYuvWrWLlypWiYsWKomzZssLLy0u7Xnp6umjXrp2wt7cX06dPF7t27RI//vijKFOmjKhWrZp48uSJEEKImJgYYWNjI1q3bi02bdokIiIixKpVq0S/fv3Ew4cPhRBCpKamihYtWggLCwsxbtw4ERYWJrZs2SImT54s1qxZo/3Md999V1haWooPPvhA7NixQ6xevVr4+voKV1dXERcXp10vKChIODs7i0qVKolFixaJXbt2iffee08AED///LMQQoj4+Hjtvn788cfaY3P9+vXc/KkQUSHo0KGDsLS0FLdu3cr1NiEhIcLe3l6cOnVKtG7dWhQrVkwEBAQIIYTYuXOneO2110SZMmWEtbW1qFChghg8eLC4e/dupvfZtm2bqF27trCyshLe3t7iyy+/FFOnTs3xPC6EPL988MEHwtvbW1haWgoPDw8xatQo8fjxY531AIjhw4eLFStWCF9fX2Fraytq1aoltm7dql1H85kv3sLDw196HJYvXy4qV64srKyshK+vr/j5559FSEiIznlcCCGSk5PFjBkzRJUqVYSVlZUoVaqUGDBggLhz547Oert37xZBQUHCyclJ2NjYiLJly4pu3bqJpKQk7TrPnj0T06dPF76+vsLa2lo4OTmJ5s2biwMHDmjXUavVYv78+aJ27drCxsZGlChRQnTv3l1bbmkEBQWJ6tWri7/++ks0adJE2NraCh8fHzFr1iyRnp4uhBAiPDw8y2MzderUlx6bmjVrih49eugsy+44A9D5fu/fvy+GDRsmPDw8hKWlpfDx8RGTJ08Wz54903m/p0+fiokTJ+r8Dbz33nvack/Dy8tLdOzYUWzdulXUqVNH2NjYCF9fX+3fwPLly4Wvr6+ws7MT9evXF0eOHMm0P0eOHBGdO3cWJUuWFNbW1qJOnTpi7dq1mdaLiooSjRo1EtbW1sLd3V1MnDhRLFmyRAAQMTExLz1m2f3td+zYUfz+++/Cz89P2NjYiCpVqoilS5dm2j45OVk4OjqKKVOmvPRzTB0THcoXcXFxAoDo3bt3ptfS0tJEamqq9vZ88hASEiIAiGXLlr30/dVqtUhNTRVXr14VAMTmzZu1rzVs2FB4eHiIp0+fapclJCQIJyenHAvIWbNmCTMzs0wnsvXr1wsAIiwsTLsMgHB1dRUJCQk6+21mZiZmzZqlXfbll1/m6iQmhExePDw8RN26dXWOy5UrV4SlpaVOAblmzRoBQISGhuq8x5EjRwQAsWDBAp3YT5w4ke3nrlixQgAQP/zwQ7brREVFCQDiq6++0ll+/fp1YWtrKz788EPtsqCgIAFAHD58WGfdatWqibZt22aKdfny5dl+LhEpgxeseMGqIC5YXb9+XaeMen65ZnvNbfz48QKA+OKLL4QQMnmpVauWsLe3F3PnzhU7d+4UU6ZMERYWFqJDhw7a91Kr1aJt27bCwsJCTJkyRezcuVPMnTtX2NvbCz8/P52kyMvLS3h6eooaNWqINWvWiLCwMNGwYUNhaWkpPvnkE9G4cWOxYcMGsXHjRlG5cmXh6uqq/V6EEGLPnj3CyspKNG3aVKxdu1bs2LFDDBgwIFPZ9s8//wg7OztRrVo1sWbNGrF582bRtm1b7d90XhMdT09PUa1aNbFixQrxxx9/iB49eggAIjIyMtN7tG/fXtStW/eln2PqmOhQvnhZolO7dm2dKzVffvml9jVNopNVAXL79m0xZMgQ4enpKczMzHTeY/bs2UIIIR4/fizMzMzEiBEjMm2vee/nvVhANm7cWNSqVUsnEUtNTRWJiYlCpVLp/JjPbv/c3NzE0KFDtc/1SXTOnDkjAIi5c+dmei0oKEingOzbt68oUaKESElJyRSvm5ub6NmzpxBCiEuXLgkrKyvRoEED8dNPP2W6aieEEG+++aawsbHRXqXLykcffSRUKpW4fft2ps8LCAgQDRo00InVzc0t03v07t1b+Pr6ap8z0SEyXLxgxQtWBXHBau3atQKAOHTo0EvX27dvn7CxsRF9+/bVHsdFixYJAOK3337TWXfOnDkCgNi5c6cQQogdO3boJEgvfvaSJUu0y7y8vIStra24ceOGdtmJEycEAOHu7q5TY7Zp0yYBQGzZskW7zNfXV/j5+YnU1FSdz+rUqZNwd3fXlqu9evUStra2OslkWlqa8PX1faVEx8bGRly9elW77OnTp8LJyUkMGTIk03t89NFHwszMLFPNZlHCPjqUL0qVKgVbW1tcvXo102urV6/GkSNHsGXLliy3tbOzy9SxVK1Wo02bNtiwYQM+/PBD7N69G3/99RcOHToEAHj69CkA4OHDh1Cr1XBzc8v0vlkte9Ht27dx6tQpWFpa6twcHBwghMjUH8jZ2TnTe1hbW2vj0df9+/ezjfXFZbdv38ajR49gZWWVKd64uDhtrBUqVMCff/4JFxcXDB8+HBUqVECFChXwzTffaN/r7t278PDwgJlZ9qeA27dvQwgBV1fXTJ936NChAj82RGQ4/P39dc4BX331VaZ1unfvnmnZnTt3MHToUJQtWxYWFhawtLSEl5cXAODs2bMAgKSkJBw5cgTdunWDjY2NdlsHBwd07tw5x9i2bduGGjVqoE6dOkhLS9Pe2rZtm2UfyRYtWuh0iHd1dYWLi0uW5VdunD9/Hrdu3UKfPn10Rsny8vJCo0aNMsVaokQJdO7cWSfWOnXqwM3NTRtrnTp1YGVlhcGDB+Pnn3/O1C8SAH7//XfY2Njg7bffzja2bdu2QaVS4a233tL5PDc3N9SuXTvTsXFzc0ODBg10ltWqVSvPxwaQgxcAgIuLS7brnD17Fq+99hoaNWqEZcuWaY/jnj17YG9vjzfeeENnfU1f2927d2vXe365Ro8ePWBvb69dT6NOnTooU6aM9nnVqlUByP6mzw+GpFmu2f9Lly7h3Llz6Nu3LwDoHNMOHTogNjYW58+fBwCEh4cjODgYrq6u2vczNzdHr169sj0OuVGnTh2UK1dO+9zGxgaVK1fO8jtycXGBWq1GXFzcK32mMeOoa5QvzM3N0bJlS+zcuROxsbFwd3fXvlatWjUAyHYoxayGT/z7779x8uRJ/PTTTwgJCdEuv3Tpks56JUuWhEqlyvKfODf/2JoEbdmyZdm+XpA0yUFu4i9VqhScnZ2xY8eOLN/r+YK7adOmaNq0KdLT03H06FF89913GD16NFxdXdG7d2+ULl0a+/fvh1qtzjbZKVWqFFQqFfbt2wdra+tMr2e1jIiMV04XrJ48eYLY2Fi89tprmV5/2QWrW7duYcqUKahZsybs7e2hVqsREBCQrxesLl26BEtLyyxfV/qC1fNl3/MXrF4Wq+aC1RdffIHhw4cjKSkJ5cuXx8iRIzFq1CgA+l+wykr58uV1nhfEBSvNts8nsc+7desW2rVrB09PT2zYsEHn2Ny/fx9ubm6Zfie4uLjAwsJCe+zv378PCwsLlC5dWmc9lUoFNzc37XoaTk5OOs81n5nd8mfPngGQxxMAxo0bh3HjxmW5P5rvUBP7i3LzN/0y+nxHmmNelC84MtGhfDNp0iT8/vvvGDp0KNavX59toZMbmpPaiz+mFy9erPPc3t4eDRo0wIYNG/Dll19q/6kTExOxdevWHD+nU6dOmDlzJpydneHj45PneJ+niTk3J5YqVarA3d0da9aswdixY7X7ffXqVRw8eFBnhLpOnTrh119/RXp6Oho2bJirWMzNzdGwYUP4+vpi1apVOH78OHr37o327dtjzZo1+Omnn7K9GtipUyfMnj0bN2/eRM+ePXP1eTnR59gQUeHiBau84QWrl9Mc/wcPHuj8TQFy1NMOHTpArVYjLCwMjo6OOq87Ozvj8OHDEELo/I3duXMHaWlp2vd2dnZGWloa7t69q5PsCCEQFxenHfEzv/Zl0qRJ6NatW5brVKlSRRtTXv+m88uDBw8AFPz/gCFjokP5pnHjxpg/fz7ef/991K1bF4MHD0b16tVhZmaG2NhYhIaGAkCu5j/w9fVFhQoVMHHiRAgh4OTkhK1bt2LXrl2Z1p0xYwbatWuH1q1b44MPPkB6ejrmzJkDe3t77T95dkaPHo3Q0FA0a9YMY8aMQa1ataBWq3Ht2jXs3LkTH3zwQa6TCo2aNWsCAL755huEhITA0tISVapUyXLuADMzM8yYMQODBg1C165d8e677+LRo0eYNm1apqs+vXv3xqpVq9ChQweMGjUKDRo0gKWlJW7cuIHw8HB06dIFXbt2xaJFi7Bnzx507NgR5cqVw7Nnz7Q/AFq1agUAePPNN7F8+XIMHToU58+fR4sWLaBWq3H48GFUrVoVvXv3RuPGjTF48GAMHDgQR48eRbNmzWBvb4/Y2Fjs378fNWvWxLBhw/Q6NprJB1etWoWqVauiWLFi8PDw0EnoiEg5vGAFnZh5wSozfS9Y+fr6AgAuX76M6tWra5enpKSga9euuHLlCvbv3w9PT89M2wYHB+O3337Dpk2b0LVrV+1yzRQMmukogoOD8cUXX2DlypUYM2aMdr3Q0FAkJSVp13tVVapUQaVKlXDy5EnMnDnzpeu2aNECW7Zswe3bt7U1aunp6Vi7dm2+xJIb//77L5ydnbOt0SsKmOhQvho6dCgCAwPxzTff4H//+x9u3boFlUoFT09PNGrUCLt370bLli1zfB9LS0ts3boVo0aNwpAhQ2BhYYFWrVrhzz//1GmbCgCtW7fGpk2b8PHHH6NXr15wc3PDe++9h6dPn2L69Okv/Rx7e3vs27cPs2fPxpIlSxATE6Odd6ZVq1Z5mtG6efPmmDRpEn7++Wf88MMPUKvVCA8PR/PmzbNc/5133gEAzJkzB926dYO3tzcmT56MyMhInfbT5ubm2LJlC7755hv88ssvmDVrFiwsLODp6YmgoCBtglWnTh3s3LkTU6dORVxcHIoVK4YaNWpgy5YtaNOmDQDAwsICYWFhmDVrFtasWYN58+bBwcEBtWvX1pmgbPHixQgICMDixYuxYMECqNVqeHh4oHHjxpnaceeGnZ0dli1bhunTp6NNmzZITU3F1KlTOZcOkYHgBSuJF6yyp+8Fq4YNG8LW1haHDh3SafY4ZswY7NmzBzNnzsTjx4+1fXABoHTp0qhQoQL69++P+fPnIyQkBFeuXEHNmjWxf/9+zJw5Ex06dNAei9atW6Nt27aYMGECEhIS0LhxY5w6dQpTp06Fn58f+vXrp9c+vszixYvRvn17tG3bFgMGDECZMmXw4MEDnD17FsePH8e6desAAB9//DG2bNmCli1b4pNPPoGdnR3mz5+faY6+gnTo0CEEBQVlWeNaZCg4EAIREREZoBMnToiBAwcKHx8fYW1tLWxsbETFihVF//79xe7du3XW1cyjk5UzZ86I1q1bCwcHB1GyZEnRo0cPce3atSznXtmyZYuoVauWsLKyEuXKlROzZ8/O9Tw6jx8/Fh9//LF2bhpHR0dRs2ZNMWbMGJ1Rr/DfPDovyuo9J02aJDw8PLSjfuY0j86PP/4oKlWqJKysrETlypXFsmXLspxHJzU1VcydO1c7r02xYsWEr6+vGDJkiLh48aIQQo6W1rVrV+Hl5SWsra2Fs7OzCAoK0hn9Swg54tYnn3yi/VxnZ2fRsmVLcfDgQZ31li1bJho2bCjs7e2Fra2tqFChgujfv784evSodh3NPDovymof1qxZI3x9fYWlpWWu5tHp16+fqFatms4yzShvWd1enEdn6NChwt3dXVhYWAgvLy8xadKkLOfRmTBhgvDy8hKWlpbC3d1dDBs2LNt5dF6U1d9GTExMptFihRDi5MmTomfPnsLFxUVYWloKNzc30bJlS7Fo0SKd9Q4cOCACAgKEtbW1cHNzE+PHj8+XeXReFBQUJIKCgnSWXbp0KcsR/ooalRBCFH56RURERERFwdGjR1G/fn0cOnRI79o1ypspU6ZgxYoVuHz5Miwsim4DLiY6RERERFSgevXqhaSkJGzbtk3pUEzeo0ePUL58eXz33XfaobCLKs6jQ0REREQF6quvvkL9+vWRmJiodCgmLyYmBpMmTUKfPn2UDkVxeiU6CxcuRK1atVC8eHEUL14cgYGB+P3331+6TWRkJPz9/WFjY4Py5ctj0aJFrxQwERHR81g2ERk+T09PTJ06NcsBHSh/+fn5Yfz48UV7EIL/6JXoeHp6Yvbs2Th69CiOHj2Kli1bokuXLvjnn3+yXD8mJgYdOnRA06ZNER0djcmTJ2PkyJHaUVuIiIheFcsmIiLKyiv30XFycsKXX36pHSL3eRMmTMCWLVtw9uxZ7bKhQ4fi5MmTiIqKepWPJSIiyhbLJiIiyvMwDOnp6Vi3bh2SkpIQGBiY5TpRUVHaeTs02rZti6VLlyI1NTXbiciSk5ORnJysfa5Wq/HgwQM4OzuzGo6IqBAJIZCYmAgPD49sZ183JCybiIhMX27LJr0TndOnTyMwMBDPnj1DsWLFsHHjRlSrVi3LdePi4jLNxurq6oq0tDTcu3cP7u7uWW43a9asHCd6JCKiwnP9+vUsZy43FCybiIiKnpzKJr0TnSpVquDEiRN49OgRQkNDERISgsjIyGwLlBevcmlayr3s6tekSZMwduxY7fP4+HiUK1cO169fz9VszEREpkQIwNcXiIsDwsKAxo0L77MTEhJQtmxZg+9AzLKJiKjoyG3ZpHeiY2VlhYoVKwIA6tWrhyNHjuCbb77B4sWLM63r5uaGuLg4nWV37tyBhYUFnJ2ds/0Ma2trWFtbZ1quGVGHiKgouXpVJjkWFkBQEGBnV/gxGHrTLJZNRERFT05l0ys3uBZC6LRZfl5gYCB27dqls2znzp2oV69etm2giYhI18GD8t7PT5kkxxixbCIiIr0SncmTJ2Pfvn24cuUKTp8+jY8++ggRERHaWVcnTZqE/v37a9cfOnQorl69irFjx+Ls2bNYtmwZli5dinHjxuXvXhARmTBNotOokbJxGCqWTURElBW9mq7dvn0b/fr1Q2xsLBwdHVGrVi3s2LEDrVu3BgDExsbi2rVr2vV9fHwQFhaGMWPGYP78+fDw8MC3336L7t275+9eEBGZMCY6L8eyiYiIsvLK8+gUhoSEBDg6OiI+Pp7toImoSHn8GChRAkhPB65fBwp74DOef7PHY0NEpIzcnn8Nf1IEIqIi7MgRmeSULVv4SQ4REZExY6JDRGTA2GyNiIgob5joEBEZsAMH5D0THSIiIv0w0SEiMlDp6cD+/fJx06bKxkJERGRsmOgQERmokyeBxESgeHGgVi2loyEiIjIuTHSIiAzUvn3yvkkTwNxc2ViIiIiMDRMdIiIDtXevvG/WTNk4iIiIjBETHSIiAyRERqLD/jlERET6Y6JDRGSAzp0D7t0DbGyAevWUjoaIiMj4MNEhIjJAmv45gYGAlZWysRARERkjJjpERAaI/XOIiIheDRMdIiIDIwQQGSkfs38OERFR3jDRISIyMFevAjduABYWQECA0tEQEREZJyY6REQGRtM/p149wN5e2ViIiIiMFRMdIiIDw/45REREr46JDhGRgeH8OURERK+OiQ4RkQG5fRu4cAFQqYDGjZWOhoiIyHgx0SEiMiAREfK+Vi2gZElFQyEiIjJqTHSIiAzInj3yvmVLZeMgIiIydkx0iIgMSHi4vG/RQtk4iIiIjB0THSIiA3HjBnDxImBmxhHXiIiIXhUTHSIiA6GpzfH3BxwdlY2FiIjI2DHRISIyEOyfQ0RElH+Y6BARGQAhMhId9s8hIiJ6dUx0iIgMQEwMcO0aYGEBNGmidDRERETGj4kOEZEB0PTPadgQsLdXNhYiIiJTwESHiMgAsNkaERFR/tIr0Zk1axbq168PBwcHuLi44PXXX8f58+dfuk1ERARUKlWm27lz514pcCIiUyFERo0OByLQH8smIiLKil6JTmRkJIYPH45Dhw5h165dSEtLQ5s2bZCUlJTjtufPn0dsbKz2VqlSpTwHTURkSs6fB2JjAWtrIDBQ6WiMD8smIiLKioU+K+/YsUPn+fLly+Hi4oJjx46hWQ6z27m4uKBEiRJ6B0hEZOo0tTmNGgE2NsrGYoxYNhERUVZeqY9OfHw8AMDJySnHdf38/ODu7o7g4GCEa0p1IiJi/5x8xrKJiIgAPWt0nieEwNixY9GkSRPUqFEj2/Xc3d2xZMkS+Pv7Izk5Gb/88guCg4MRERGR7ZW25ORkJCcna58nJCTkNUwiIoOWng7s3i0fBwcrG4spYNlEREQaKiGEyMuGw4cPx/bt27F//354enrqtW3nzp2hUqmwZcuWLF+fNm0apk+fnml5fHw8ihcvnpdwiYgM0l9/ySGlixcH7t+X8+gYkoSEBDg6OhrN+ZdlExGR6ctt2ZSnpmvvv/8+tmzZgvDwcL0LEgAICAjAxYsXs3190qRJiI+P196uX7+elzCJiAzerl3yPjjY8JIcY8OyiYiInqdXsSqEwPvvv4+NGzciIiICPj4+efrQ6OhouLu7Z/u6tbU1rK2t8/TeRETGZOdOed+6tbJxGDOWTURERcvjx7lbT69EZ/jw4Vi9ejU2b94MBwcHxMXFAQAcHR1ha2sLQF7xunnzJlasWAEAmDdvHry9vVG9enWkpKRg5cqVCA0NRWhoqD4fTURkchITgYMH5eM2bZSNxZixbCIiKjri44Fu3XK3rl6JzsKFCwEAzZs311m+fPlyDBgwAAAQGxuLa9euaV9LSUnBuHHjcPPmTdja2qJ69erYvn07OnTooM9HExGZnMhIIC0NKF8eqFBB6WiMF8smIqKi4f59oF074OjR3K2f58EICpOxdYYlIsqNkSOB774DhgwBFi1SOpqs8fybPR4bIqLCc/OmbP1w5gzg5JSABw8KaDACIiJ6dZr+OWy2RkRElL2LF4HGjWWS4+EB/P577rbjGD9ERAq4dg04fx4wMwNatlQ6GiIiIsN04gTQti1w5w5QsaIcrTQX80EDYI0OEZEiNMNKN2gAlCihaChEREQGaf9+oHlzmeTUqSOfe3vnfnsmOkRECmCzNSIiouyFhckyMj4eaNIECA8HXF31ew8mOkREhSw9HfjzT/mYiQ4REZGuNWuALl2Ap0+BDh2AP/7IW+sHJjpERIXs+HHgwQPAwUE2XSMiIiJpwQKgb185/UKfPsCmTYCdXd7ei4kOEVEh04wWExwMWFoqGwsREZEhEAKYPh0YPlw+Hj4c+OWXVysnOeoaEVEhCwuT9x07KhsHERGRIUhLA957D/jhB/l8yhSZ9KhUr/a+THSIiArR3bvAX3/Jxx06KBsLERGR0p48AXr3BrZulYnN/PnAsGH5895MdIiICtGOHbJKvk4dOekZERFRUXXvHtC5M3DoEGBjA6xeDXTtmn/vz0SHiKgQaZqtsTaHiIiKspgYoF074MIFoGRJWaPTuHH+fgYTHSKiQpKWJofIBJjoEBFR0RUdLcvBuDigXDnZ2qFq1fz/HI66RkRUSA4fBh4+lFeuGjZUOhoiIqLCt2sX0KyZTHJq1QIOHiyYJAdgokNEVGg0zdbatQMsWJ9ORERFzMqVsibn8WOgRQtg716gTJmC+zwmOkREhWT7dnnPZmtERFSUaObI6ddPNuPu3VvOKefoWLCfy2uKRESF4OZN4ORJOXRm27ZKR0NERFQ4kpOBQYNkbQ4AjB8PzJ4NmBVCdQsTHSKiQvD77/K+QQOgdGllYyEiIioM9+/L4aL37QPMzYEFC4DBgwvv85noEBEVAk2ztY4dlY2DiIioMFy8KMu8ixeB4sWBdeuANm0KNwYmOkREBezpU2DnTvmYiQ4REZm6ffuA118HHjyQw0dv3w7UqFH4cXAwAiKiArZ7N/DkCVC2LODnp3Q0REREBWf1aqBVK5nk1Ksnp1ZQIskBmOgQERW4TZvk/WuvycEIiIiITI0QwKefAn37Aikpsm9OZCTg5qZcTGy6RkRUgNLTga1b5ePXX1c0FCIiogLx5AkwcCDw22/y+bhxwJw5hTOy2ssw0SEiKkCHDwN37si5AoKClI6GiIgof924AXTpAhw/LifDnj+/cEdWexkmOkREBWjzZnnfoQNgaalsLERERPnp0CHZWuH2bcDZGQgNNayLeuyjQ0RUgDT9c9hsjYiITMmKFTKpuX0bqFkTOHLEsJIcgIkOEVGBOXcOuHBB1uS0a6d0NERERK8uPR348EMgJEQOOtClC3DgAODjo3RkmTHRISIqIJpmay1bysnSiIiIjFlCghxB9Msv5fOPPgI2bAAcHJSNKzvso0NElN/S04F9+7DpxxoASqFLZzV4XYmIiIzGf+UYYmMBd3egaVNcuGyO118Hzp4FbGyA5cuB3r2VDvTl9Cp5Z82ahfr168PBwQEuLi54/fXXcf78+Ry3i4yMhL+/P2xsbFC+fHksWrQozwETERm0DRsAb2/EteiNw5ecAACvfd5QLqcCwbKJiCgf/VeOoUULoE8foEULbHF9F/X9UnH2LFCmjMyBDD3JAfRMdCIjIzF8+HAcOnQIu3btQlpaGtq0aYOkpKRst4mJiUGHDh3QtGlTREdHY/LkyRg5ciRCQ0NfOXgiIoOyYQPwxhvAjRvYhNchYIb6+Atl4o7J5Ux2CgTLJiKifPJcOQYAaqgwFdPQ5f4yJDyxRJOq93D0KFCvnsJx5pJKCCHyuvHdu3fh4uKCyMhINGvWLMt1JkyYgC1btuDs2bPaZUOHDsXJkycRFRWVq89JSEiAo6Mj4uPjUZwN3YnIEKWnyytg/xUOwfgTexCMOfgQH+JLQKUCPD2BmBjA3FzZWPVgjOdflk1ERHnwQjn2ECXwFlYiDB0BAO/jO8z1nAerKxcUL8dye/59pT468fHxAAAnJ6ds14mKikKbNm10lrVt2xZLly5FamoqLLOYWCI5ORnJycna5wkJCa8SJuVRcjJw8qQcOeriReDmTeDRIyA+Xv4vWFvLW6lS8vebpydQvTpQq5bhdkojKjD79mkLhzsojQg0BwD0wDr5uhDA9etyvebNlYmxiGDZRESUB8+VY6dRA12xEZdRETZ4iiUYjH5YCdyAUZVjeU50hBAYO3YsmjRpgho1amS7XlxcHFxdXXWWubq6Ii0tDffu3YO7u3umbWbNmoXp06fnNTTKIyHkrLabNwN//gkcOyaHDdSXSgVUqiTHUm/VSo44VapU/sdLZFBiY7UPN6Ir1DCHP47CB1eyXY/yH8smIqI8unULAPAreuEdLMUT2MMbMdiAbvDDiYz1jKgcy3OiM2LECJw6dQr79+/PcV2VSqXzXNNa7sXlGpMmTcLYsWO1zxMSElC2bNm8hko5uHYNWLZMjp5x7Zrua6VLy0mgKlYEvLyAEiUAR0c5L0hyMvDsGXDnjrwAcOUKcPq0rPm5cEHefvgBMDOT/dl69QK6dZMz5xKZnOd+GK9DDwDP1eZksx7lP5ZNRER6EgLYtQtpn8/BBMzF1/gAANAGf2A1+sAZD3TXN6JyLE+Jzvvvv48tW7Zg79698PT0fOm6bm5uiIuL01l2584dWFhYwDmbX7zW1tawtrbOS2ikh+ho4PPPZb8zTU8te3s5sWGnTkDTpkD58rKGRh937wJ//QXs3g3s2gX8/bd8vHs38N57MtkZPly+v77vTWSwmjYFPD1x98YzhKMFgBcSHU0fnaZNFQrQ9LFsIiLSgxDyx9m0abh14F/0xq/YB9mvcRJmYgamwBzqjPWNsBzTa9Q1IQRGjBiBDRs2YM+ePfDJxRSogYGB2LVrl86ynTt3ol69elm2gaaCd/asnOypbl0gNFT+nbdoAaxeLZOU9euBAQOAChXyloiULg107Ah8/bWs4bl8GZg1C6hTB0hLA377TTZrq10b+OknIDU1n3eQSAnm5sCMGTrN1sojRr6m+UeaN0/xDpymiGUTEZEehAD27JE/xlq3xq4DtqiDE9iHZnCwSUEoumOm6uPMSQ5gfOWY0MOwYcOEo6OjiIiIELGxsdrbkydPtOtMnDhR9OvXT/v833//FXZ2dmLMmDHizJkzYunSpcLS0lKsX78+158bHx8vAIj4+Hh9wqUXPHwoxOjRQlhYCAEIYWYmRJ8+Qpw+XXgxnDghxLvvCmFnJ2MAhPDxEWLxYiGSkwsvDqIC8eGHohV2CkCI2fgw44+8bFkhQkOVji5PjOH8y7KJiCiXwsOFaNZMCECkwUxMNZ8hVEgXgBC1awtx8aKQ5ZWnZ0YZZoDlWG7Pv3olOgCyvC1fvly7TkhIiAgKCtLZLiIiQvj5+QkrKyvh7e0tFi5cqM/HsjDJB3/8IUSZMhl/r6+9JsS5c8rF8/ChEHPmCOHiopvwrFsnhFqtXFxEeRYXJ+7alhXmSBWAEJdXRQmxerUsVNLSlI4uz4zh/MuyiYgoB5GRQjRvrv3RFWfpKYI9z2l/gw0eLMRz14ZkuRUebrDlWG7Pv680j05h4VwFeZeUBIwfDyxcKJ9XrAjMnw+8MKqqYp48AZYsAebMATTN5Rs3Br76CmjYUNnYiPQyZgx+mPcYg/ED6tYVOHbMNDqg8fybPR4bIjJ4+/cDU6fKpmoAYGWFvR1mo/fBkYi9Yw47O2DxYuCtt5QNU1+5Pf/q1UeHjMvly0BgYEaSM2IEcOKE4SQ5AGBnB4weLefpmTpVPj9wAAgIAN59F3j4UOkIiXLh1i1g4UL8hp4AgB49TCPJISIiI3XwINC6tRw4YM8ewNIS6iHDMGfsbbTcOgaxd8xRrRpw5IjxJTn6YKJjonbsAOrVk4MBuLrK0c+++06OqmaIihUDpk2TQ1IPGCCX/fgj4OsLrFmTMSockUGaORNxySWwBy0BAD16KBwPEREVTYcOAW3byuYxf/4JWFgAgwfjTtRldLq2ABNnl0B6OtCvnxwht1o1pQMuWEx0TND8+UCHDsCjR7Jm5NgxOXGnMShTRs7ns28fULWqnKOnTx853LURzU9FRcm1a8APP+BX9IYa5ggIkCMWEhERFZrDh4H27WVTnp07ZYLz7rvAxYvY2X0xanUsi99/B2xs5ByHP/9suBe/8xMTHRMiBPDJJ7KJmhDAoEFARIRMHoxNkyZynp9PPwWsrICwMDlx6YYNSkdG9ILPPwdSUrDKYSgAoG9fheMhIqKi48gReXU7IEA25zE3B955B7hwASnfL8GHC7zRti1w+zZQo4ZcfdCgojOPIRMdE6FWA8OGATNmyOfTpslO/sY8t521NTBlCnD8uJyD5/59oHt32bQtIUHp6IgA/PsvsGwZLqASjib6wtwc6NlT6aCIiMjkHTsmm7s0aAD8/rtMcAYOlH0AfvwRF9N80Lgx8OWXcvX33pNN1WrUUDbswsZExwSo1cDQoXLUDJVKDj4wdarpZOvVq8sa2UmTADMzWd1at66s8SFS1IwZQFoaVpX/BIAc6MPFReGYiIjIdB0/Lmd9r1cP2L5d/jAKCQHOnQOWLQPKl8cvv8jfSUePAiVLAhs3ym4NtrZKB1/4mOgYOSGAUaNke0szM2DlSpn0mBorK2DmTCAyEvDyyhhRbskSDlRACrl4EVixAgLAqpQ3ALDZGhERFZATJ4DXXwf8/YGtW+WPvn79ZILz009AxYpISJAjqPXvDzx+DDRrBpw8KTcrqpjoGLkJE4Dvv5e1N8uWyY77pqxJE3kxo1MnIDkZGDJEXshISlI6Mipypk8H1Gr81WgMLt+wgZ0d0KWL0kEREZFJOXUK6NYN8PMDNm+WCU7fvsCZM8CKFUClSgDk1Bx16gCrVslWbDNmyFGly5ZVNnylMdExYt9+m9H2cvFi+YO/KHBykv/rs2fLf+ZffpGTi/77r9KRUZFx5gywejUAYFXZiQDkFbNixRSMiYiITMfp08AbbwC1a8u2ZyqVvJr9zz+y+U6VKgCAlBTZtL9ZMyAmRrZ6iYwEPv5Y/kYq6pjoGKnNm+VEm4D8wf/uu4qGU+jMzGRt1p49gJub/L9v0ECOMkdU4KZNA4RAapc38Ose2SmHzdaIiOiV/f23nIytVi0gNFQmOL17yx86q1bJCQafW7VhQ/k7UK2WF7xPnpRT6JDERMcIHT0qk3ohgMGDgQ8/VDoi5TRrJo9HvXpyVLbWrYFFi5SOikzaqVPAunUAgLCWX+LuXTkAQevWCsdFRETG68wZoFcvmeCsXy+X9ewpa3bWrJGTC/5HrQa+/lr+9jlxAnB2ljnRTz8Bjo6KRG+wmOgYmdu3ZT+AJ0+Adu3kKBqmMrpaXpUpA+zdC7z5JpCWJofZHj4cSE1VOjIySVOnyvuePbF8jzcA2fHT0lK5kIiIyEidPSt/wNSoAfz2m7yK/cYbMsFZu1YOPfuca9fkJPAffCD7KnfsKGt2unVTKH4Dx0THiKSlyWT/1i2Z2K9dKye+JTlk4qpVcmQ2AFiwQA5YkJiobFxkYo4dAzZtAszMcHvEDGzbJhcPHKhoVEREZGzOn5dtnqtXB379VSY43brJtmfr1mWa8EYI2Se5Zk0gPBywt5f9s7dulU34KWtMdIzIpEmyg1mxYsCGDUDx4kpHZFhUKnmMNm0C7OyAnTtl07bYWKUjI5Ohqc3p0wcr/6qM9HTZPrpaNWXDIiIiI3HhghwWulo1OaiNEHI0m+ho2f6sVq1Mm9y6JVvz9O8vJ0wPDJRN1gYPZquenDDRMRKhocDcufLxTz/p9EWjF3TpIgclKF1anggCAmTNMNErOXRITs5mbg4x5RMsWyYXszaHiIhydPGiHC2galU5appaLSf+PH5cjqpWp06mTYSQk6RXry5rbiwtgc8/l831K1Ys/F0wRkx0jMD168CgQfLx+PFA9+7KxmMM6tcHoqLk8PLXrgGNGgH79ikdFRm1Tz6R9/3748ijSjhzBrCxkYPhEBERZenyZWDAAJngrFghE5zOneVISps3y/lxsnDzpmyCP2AA8OiRHHjg+HFg8mR2W9AHEx0Dl54uqyofPZLDJ3/+udIRGY8KFYCDB2UV76NHsvNeaKjSUZFR2rcP2LVLli5Tpmhrc7p35wg3RESUhX//Bd5+W8538/PP8gddx47AkSPAli2Av3+WmwkhJ4CvXh0ICwOsrIBZs+TF2xe67VAuMNExcF99JZth2dvLmk6O7KSfUqWA3buBrl3lpFo9ewLLlysdFRkdTW3O22/jqZsPfv1V+5SIiCjDlSuyGU6VKvIHR3o60L49cPgwsG2brJrJxvXrctV33gHi42Uf0OhoYOJE1uLkFRMdAxYdLWe2BYBvvpHNsEh/trZyAJN33pE1xm+/Dcybp3RUZDT27JFXG6ysgI8/xm+/yQLI2xto3lzh2IiIyDBcvSpHB6hUCVi6VA6V27atrIoJC5PNcrKhVsvRYqtXB/74A7C2Br78EjhwgIPdvCrmhwYqNVX+IE9NlbURvHL8aszNgR9+AEqUkLVkY8YADx/KCe45YgllSwhgyhT5ePBgoGxZLFyY8dSMl4qIiIq2a9fk3BbLlmVM4Ne6NTB9umw7n4O//5blSVSUfB4YKN+Kg07lDxbTBmruXDlimJMTsHAhf4znB5VKXiH57DP5/NNPgdGj5ZUUoizt3Ck7etnYAJMmITpatj6wtJQ1hEREVERdvw68954c/mzxYpnkBAcD+/fLsiOHJOfpU+Cjj+RYBFFRgIMD8N13sksok5z8wxodA3T+vLwQAMgmVq6uioZjUlQqeWIpUQIYMQL49lvZDOnHH9n+lV7wfG3OsGGAhwcW/fd/2b074OKiXGhERKSQGzfk6AA//ig7/wJAy5ayiUjTprl6iz17gCFDgEuX5PMuXYDvvwc8PQsm5KKMNToGRq2WfdiSk2XTzrfeUjoi0zR8uBzl0dxcDobSr59sTkuktW2bHB3Hzg6YOBEJCcCqVfKloUOVDY2IiArZrVvA++/LIV0XLJBJTvPmsg/n7t25SnLu35dzrwUHyyTHw0NOAL9pE5OcgsJEx8AsXy5rPe3tZU0om6wVnH79gPXrZTOkX3+V86FomtdSEadWZ4y09v77gIsLVq4EkpLkVAjNmikbHhERFZLYWGDUKKB8eVntkpIik5o9e4DwcCAoKMe3EEJeXPX1lZO+q1TyguuZM7IfNhUcJjoG5OFDOYQgIJuueXkpG09R8Prr8mqKlZWcY6dnz4yaaCrCNm2SneSKFQPGjYMQ0A5CMHQoL0AQEZm8uDg5clH58rKde3Iy0KSJrL2JjARatMjV25w8KS+OhYQA9+7JuXAOHJA5E+dhK3hMdAzI1Knyn6BqVWDkSKWjKTo6dZKTE1tby9+33bvL8xkVUWq1/GcE5GgVpUrhwAE5Mo6dnZzAl4iITNTt28AHH8gEZ9484NkzoFEjOWn03r2yP04urnbFx8uKoLp1M1rqzJkDHDuWq8HYKJ8w0TEQp0/LJp+AvHDAiUELV7t2wNatcnCtbdtkTc/Tp0pHRYr47TeZ1Tg6AmPHAsiYd6lPHzmQBRERmZg7d4Dx4wEfH+Drr+WPgIAAObHN/v1Aq1a5SnCEAH75Rc4X+u238tpZz57AuXPAhx/KFiRUePROdPbu3YvOnTvDw8MDKpUKmzZteun6ERERUKlUmW7nzp3La8wmRwjZDSA9XdYmtGqldERFU+vWck4vOztgxw45CsqTJ0pHRYUqPV2OnAPIK3olSyImBti4US4aNUqxyOglWC4RUZ7dvSszEB8fObfH06dycs/ff5fTC7Rpk+v2yqdOyWZq/fvLiqEqVWRF0Nq1HGxAKXonOklJSahduza+//57vbY7f/48YmNjtbdKlSrp+9Ema/162dzT1lZOZknKadFCntvs7eXJqVMn2QGdiojVq+X47k5O2qzm++/lFbnWrWXbajI8LJeISG/37smO0T4+cpK9J0+AevWA7duBQ4dkU49cJjgPHug2U7OzA2bPlokPL14rS++ZQ9q3b4/27dvr/UEuLi4owTYfmaSkAJMmyccffsgBCAxBs2Zyrq927eSAKq+9Jpu12dkpHRkVqNTUjAmsxo8HihdHYqKcKgGQfVLJMLFcIqJcu39fXlX+7jvg8WO5zN9f1uZ37KjXaDOpqcCiRbJb58OHctkbb8iWb2XL5n/opL9C66Pj5+cHd3d3BAcHIzw8vLA+1uD98ANw+bKcFHTcOKWjIY1GjWSy4+AgR5Ds0oV9dkzeL7/If8bSpeVssgCWLQMSEmTzg7ZtFY6P8h3LJaIi5MED4OOPZQ3OrFkyyfHzA7ZskXOmdeqkV5KzYwdQu7YcPOrhQ1njv2sXsG4dkxxDUuCJjru7O5YsWYLQ0FBs2LABVapUQXBwMPbu3ZvtNsnJyUhISNC5maLExIwLyFOnypFsyXAEBGQ0Y/vzT6BbNzn4CpmglBTg00/l4wkTgGLFkJ4uO5ICcvA1Mw7dYjLyUi4BRadsIjIpDx/KedF8fIDPP5c/vurUkcOsHjsGdO6sV4Jz7hzQoQPQvj1w9ixQqpScfiA6ms3UDJJ4BQDExo0b9d6uU6dOonPnztm+PnXqVAEg0y0+Pv4VojU8U6cKAQhRsaIQKSlKR0PZiYwUws5OflcdOgjx7JnSEVG+W7hQfsFubkIkJQkhhNi4US4qWVK7qEiKj483qvNvQZVLQhSdsonIJDx8KH9oOTrKkzkgRK1aQmzYIER6ut5vd/++ECNHCmFhId/KwkKIsWPlx1Dhy23ZpMg1yoCAAFy8eDHb1ydNmoT4+Hjt7fr164UYXeGIi5ODewDAzJkcTtqQNWsm+yba2spR2Xr04KSiJuXZM+Czz+TjyZMBOzsIITuSAnKCUPbPMn05lUtA0SibiIxefLysoff2ls1m4uNlu7L162W1S9euelXRP3smu/RUrChr+dPSZN/df/6Ry9nNz7DpPRhBfoiOjoa7u3u2r1tbW8Pa2roQIyp8M2fK0bwaNJAd18iwNW8uByTo1Ene9+olp1thgmoCfvgBuHlTjv357rsAgIgI4PBhOa8Sh5QuGnIql4CiUTYRGa2EBJmJfPUV8OiRXFa9uuwb0L273u2P09OBVauAKVOAa9fksho15EADrVvnb+hUcPROdB4/foxLly5pn8fExODEiRNwcnJCuXLlMGnSJNy8eRMrVqwAAMybNw/e3t6oXr06UlJSsHLlSoSGhiI0NDT/9sLI3LoFLFkiH3/+uV5NQ0lBwcHA5s3ySs6mTcCbbwJr1jDZMWpPnsirDgDw0Ucys0HGonfekQOFkGFjuURUhCUmyhHUvvpKDjgAAFWrygSnRw+9Exwh5ByhEybI4aEBoEwZWUkUEgKYm+dz/FSg9E50jh49ihYtWmifj/1v5vCQkBD89NNPiI2NxTVN6gsgJSUF48aNw82bN2Fra4vq1atj+/bt6NChQz6Eb5zmzAGSk4HGjeWPZzIebdrIJKdLFyA0FOjbV069YqFI3Si9soULZTtSb2/g7bcBAH/9JQefsLCQo0yT4WO5RFQEJSbKic7mzs1IcHx9MxKcPGQkx47JqT727JHPHR3lFCAjR8rm62R8VEIIoXQQOUlISICjoyPi4+NRvHhxpcN5JbGxQPnyss3nzp2s/jRW27fLZr6pqUDv3nJkYiY7RubxY/nPePcusHSpNtHp2lUmsyEhwE8/KRqhQTCl829+47EhUsDjx8D8+XKSz/v35bIqVeTIar165SnBuXBB5ke//iqfW1nJWQYmTwacnfMxdso3uT3/csDUQvbFFzLJCQzkMITGrGNH2a/RwkKeGAcOlO15yYh8/71McipUAPr1AyA7l27aJJuTTpigbHhERPScpCSZ3Pj4ABMnyiSnUiV5pfGff4A+ffROcmJi5DWuatVkWa5SAW+9BZw/L1vCMckxfkx0ClFcnJxBF5BXDtg3x7i99pockMDcHFi5UvbnYLJjJBISZIEJyH/G/zpaff65XNStm2ziTURECnvyRGYd5cvLdmX37skh0H7+GThzRmYmeiY4N24Aw4YBlSsDy5fLsrtTJ+D4cZk3eXsXzK5Q4WOiU4i+/FLW5jRsKPt6kPHr2lVeBTI3l+fcd98F1Gqlo6IczZsn23RXqSKvAgL4+++MZguTJysXGhERAXj6FPjf/2SCM24ccOeOfLx8uZyps39/vduMx8XJCaArVpQXntPSZBeCQ4fkiKp16hTInpCC2KugkDx8CCxeLB9/8glrc0zJG2/IISj79JHnX5VKjlis50AvVFgePpTjgwLAtGnaK4GffCJH2+neHahbV7nwiIiKtKdP5dC0s2fLzASQzdU+/lg2M87DUKd378oxC777Tr49IOfImzFD3pPpYqJTSBYulM1La9UC2rdXOhrKb716yR/JffsCy5bJJGfxYiY7BunrrzMmkOvZE4AcaWfjRpmkfvqpwvERERVFz55lJDixsXKZt7dMcPr3z1OCExsrW9MsXixbwAFAQIBMcIKDedG5KGCiUwiePZNzWAGy9pX/WKapd2/ZbK1fP+DHH+X3vGgRkx2Dcu+ebLYGyBmz//tyPv5YLurbV3ZKJSKiQvLsmSw0Z82SEw0CQLly8sQcEiKHQNPT1aty8KelS+V0HgBQr56sxO/Qgb/DihImOoXgl1+A27eBsmXlj2EyXX36yJqd/v1l8zWVStbmMdkxEHPnyqFJ/fxkBysA+/cDO3bIFmxTpyocHxFRUZGcLDORmTOBmzflsrJl5eTNAwfmKcG5dEnmSytWyP43gJyzcMoU2TeaCU7Rw0SngKnVcrAQQHaAy0PNKxmZvn0zkp0lS2SSM38+kx3F3b4tG2gDsjZHpYIQskwF5BCjFSsqFx4RUZGQkiLbeM+cCVy/Lpd5espRYN5+G7C21vst//lHJjhr1mQMCNSypUxwgoKY4BRlTHQK2Natcjx2R0c5IhcVDW+9JZOdkBDZfE2lkskOT7YKmjNHNtJu0ECOIwpg82Zg717Axiaj+RoRERWAlBQ5C/PnnwPXrsllHh4ywRk0SO8ERwggMlL2wQkLy1jeoYO8gNWoUf6FTsaLiU4B00zVMWwY4OCgbCxUuPr1kyfiAQNk8zWVSs5RyWRHAbduyS8BkKMNqFRISZFTMgDA2LGySTgREeWz1FQ5/8Jnn8nOMwDg7g5MmiSvANvY6PV2aWnAhg3y99XRo3KZSiVbI0+eDPj753P8ZNSY6BSgY8eAAwdkc7WRI5WOhpTQv79MdgYOBBYskM3Xvv2WyU6hmzVLdnht1Eg7idXChcDFi4CLi5xkm4iI8lFqquyk/NlnQEyMXObmJk+4gwcDtrZ6vd2TJ3IKh6+/Bv79Vy6zsZEXE8eOBSpVyt/wyTQw0SlA338v73v0kBcvqGgKCZFtht95J6NG55tvmOwUmmvXZGcpQI4pqlLhwQPZTUeziLWtRET5JC0tI8HRZCSurjLBGTJE7wQnNlZemFqwALh/Xy5zdgaGD5c3F5d8jp9MChOdAnL3ruwUBwDvv69sLKS8gQNlzc6gQbI/vEolRzlmslMIPv9ctg1v3lz2ToVsvfbwoZxK5+23lQ2PiMgkpKXJ2bNnzAAuX5bLXFyACROAoUMBOzu93u7wYdkC4rffMkZQK19e1t4MHKj321ERxUSngGjGbvf3Bxo2VDoaMgRvv52R7Hz7LZCeLu85GlsBiomRo/sA2plAT57MqG39+mvAgmdBIqK8S0uTV3ZnzJDtgQGgdGnZCXLYMMDePtdvlZICrF8vy8bDhzOWN2kiuwB07cpzNumHfy4FIC1NVrECsjaHV+1J4513ZLIzeLAche3pU9mqytxc6chM1IwZ8h+ydWugaVOo1bLcTU+XTUpbt1Y6QCIiI5WeDvz6q7yIdOGCXFaqFDB+vGxTpkeCc/u2LAsXLpRN1QA5jc6bb8rfURxggPKKiU4B2LpVDg1fqhTQq5fS0ZChGTRIdqAMCZGVDU+fygFpOMdSPrt4Uc4aB8iEB7Ija1QUUKwY8L//KRgbEZGxSk+X7ck+/RQ4d04uc3KSCc6IEfIEmwua4aEXLZKjqKWmyuXu7vKC1JAh7H9Dr46JTgHQzEmYh1ETqYh46y35t/Hmm7LG/9kzeZ+HedIoO9OnywK5Y0egYUPcu5cxnPT06UCZMsqGR0RkVNRqYN06eQI9e1Yuc3ICPvhAVrvkclSXBw/kxb3Fi+U8gxoNGwKjRgHdu8vaHKL8wEQnn509C4SHy34XQ4cqHQ0ZsjfekMnOG28AGzcCr78ur2rpOSANZeXsWWD1avn4v+HVPvhAFrA1a3KAECKiXFOrZceZ6dOBM2fkspIlMxKc4sVzfAshgIMHZe3NunWyDzMgK3/69pW1N35+BbgPVGQx0clnP/wg7zt14gSElLNOnYBt24DXXgN27JCVD1u25Lrmn7IzbZosWV9/HfD3x7ZtshWbmZm8ishmgkREOVCr5dW36dOBv/+Wy0qUkMOejRwJODrm+BZ37siB2JYty3gLAKhTR14M7tOHw/tTwWKik4+SkzO6BLz7rrKxkPFo1Qr44w+Z5ISHy/kst2+XF8woD06dku3HAWD6dDx8KK8WArJ8DgxULjQiIoOnVgObNskE59QpuczRERgzRrYtK1HipZunpgJhYbJP5PbtGUND29rK5tpDhgD163OgJiocTHTy0caNcjKrMmWAdu2UjoaMSdOmwJ9/yr+bqCj5fMcOwNNT6ciM0LRp8r5HD6BWLYwZANy6BVSpoh1hmoiIXiQEsHmzPIeePCmXFS8OjB4tk5wcEpzTp2Vys3KlnEtQo0EDYMAAmeTk8BZE+Y6JTj7SNFt7+22O8076a9BAjkDTrh3wzz9Ao0aypqdqVaUjMyLHj8srDioVMG0aNm+WnV5VKlkAs/8TEdELhJDDxU6bBkRHy2UODrL2ZswYOeBANuLiZAX6ihXAsWMZy11dgX79ZIJTvXqBRk/0Uvw5nk8uXwb27JE/qN55R+loyFjVrCk7bLZtK0ejadJENgHgpLO59Mkn8r5PH9woXg1vvy2fjhvHJmtERDqEkG3Lpk3LyFKKFZP9b8aOBZyds9wsPl5eT1q9Gti9W7Z0A2Tfx86dgYEDZRnGvpBkCJjo5JMff5T3bdoAXl7KxkLGzcsL2L9f9tn56y+gZUs54E379kpHZuAOH5aFtrk50j+eir595Shr9eoBn32mdHBERAZCCHkFbdo04OhRuczeXo6g9sEHchLAFyQny01Wr5aVP5pR0wAgIEAOKvDmm1luSqQoJjr5IDVVNosBOAgB5Y9SpWQN4RtvyL46nTvLZHrAAKUjM2Ca2px+/fDZ2krYu1denFyzhnMyEBFBCFmgTJsmr6IBgJ2dnORz3DigdGmd1ZOTZd/R0FA5+Fp8fMZrvr5yWOg33wQqVCi8XSDSFxOdfBAWBty+LWfw7dxZ6WjIVNjby6Gm335bdu4cOBC4eBGYMUMOk0zP2b8f2LkTsLDArhYz8elAuXjRIqBiRWVDIyJSlBDy/DhtGnDokFxmZwcMHw6MH6+T4Dx9Klddv16WPwkJGW9TpoxMbPr0kcNDc9Q0MgZMdPKBZkjpt97ilWPKX5aWsjN9uXLAzJnyduGCXGZnp3R0BmTKFADA5Tc+RK/R7lCrZYLYt6/CcRERKUUIWSUzdaoczhOQI7K89x7w4Yfy6iyApCTg999lcrN9O/D4ccZbuLsD3bvL1gVNm/IiGxkfJjqv6MED2V4VAEJClI2FTJOZGfD553J45EGDZGF09aocBdTdXenoDEB4OBARgceWJfF69FQ8fCgHb5g/X+nAiIgUIIRs+zx1KnDggFxmYwMMGyYTHDc33LoFbFsia2127waePcvYvGzZjOQmMJDJDRk3vf989+7di86dO8PDwwMqlQqbNm3KcZvIyEj4+/vDxsYG5cuXx6JFi/ISq0H67TfZR6d2baBWLaWjIVPWv78skJydgSNH5I/5EyeUjkphQgBTpkANFULK7sHf563g5ibblNvYKB0cFRaWS0T/CQ8HgoLkTNQHDgDW1sCoURCX/8WJ/l/j0yVuqF9fNkMbMkTW4Dx7Bvj4yFZshw/LC2n/+x/QuDGTHDJ+etfoJCUloXbt2hg4cCC6d++e4/oxMTHo0KED3n33XaxcuRIHDhzAe++9h9KlS+dqe0OnabbWv7+ycVDR0LSpbGLdqZMcfrpRI2DxYjlfQZGSng7s2wfs2AFx4ADGmH+PDf/WgaWlTHLKlFE6QCpMLJfI5GnOebGxsiq/aVPA3Dzj9chIWYMTGSmfW1sjccD7CA+YiB1/OWNbAHD9esbqKpWcu+2112Tf4ho12OeGTJR4BQDExo0bX7rOhx9+KHx9fXWWDRkyRAQEBOT6c+Lj4wUAER8fn5cwC8yFC0IAQpiZCREbq3Q0VJQ8eCBEu3by7w8Q4r33hEhOVjqqQhIaKoSnp3bn52C89jisXq10cKbHUM+/2SmsckkI4zs2ZKReOOcJQD4PDRVi714hWrQQAhDpUIljFg3EzIDNIijgmbC01N3E1laILl2E+PFH/mYh45fb82+B99GJiopCmzZtdJa1bdsWS5cuRWpqKiyzmFEqOTkZyc8N0p7w/LAfBuSXX+R927aAm5uysVDRUrIksG0b8Omn8rZgAXD8uOy/Y9K1GRs2yIbjQgAAfkZ/TMAXAICvMRZvWjcB0E3BAMkY5KVcAoynbCIT8sI5T+vGDaB7d8TBFbvQGn+YrcJO6864+9QBOJSxWoUKcn6/jh3lnGy2toUbPpHSCrz1ZVxcHFxdXXWWubq6Ii0tDffu3ctym1mzZsHR0VF7K1u2bEGHqTe1OiPRKXLNhsggmJsD06fLwTBKlJBN2urWldMkmKT0dGDUKG2B/wvewttYBgAYhy8xRjUPGD1arkf0EnkplwDjKJvIhLxwzgOAe3DGBnTFSHyDWjgJd8ShP37BKnUf3H3qgGLFZHO0+fOBS5fkbcECmegwyaGiqFC6maleaPgp/vunfXG5xqRJkxAfH6+9XX++YamB2L8fuHIFcHAAunRROhoqyjp1kpNb164N3LkDtG8vf+8/P4qOSdi3T17FhKzJCcHPUMMcg7EYczBB/hi4fl2uR5QDfcslwDjKJjIhW7fi7o1nCEU3bWJTGvfQHRvwHUbiNOQISP6VEzF5suyec/++HJHzvfc4kScRUAjDS7u5uSEuLk5n2Z07d2BhYQFnZ+cst7G2toa1tXVBh/ZKNIMQ9OjB+UxIeRUqyGkSJkwAvvsO+OYbObro6tWyk6lJiI0FAPyAQRiCxRAww1AsxHwMhxlEpvWIspOXcgkwjrKJjFB6OnDhAsSJk7i69yqiooADF10Q+aQe/sbdTKvXwGkEIRLNEYEgRKL0tG/lTJ5ElEmBJzqBgYHYqplo5j87d+5EvXr1sm0HbeiSk2VfCIDN1shw2NoC334LtGsHDBwInD4N1K8PzJgha3gsjHzWLHH9Bj7GZ5iJjwAA72E+vscIZLr+zsmFKAemWC6RkUhIAE6dAk6cwNNjZ3DsUCqiLpVGVFo9RCEIcch8/qqB02iOCDRHBJphL0rjheaVPOcRZUvvnz6PHz/GpUuXtM9jYmJw4sQJODk5oVy5cpg0aRJu3ryJFf9VeQwdOhTff/89xo4di3fffRdRUVFYunQp1qxZk397Uch27gTi4wEPD6BZM6WjIdLVoYMsR99+GwgLk3Mj/Por8OOPQJ06SkeXB8nJSP5wCt75thZW4S0AwCeYjmmYppvkqFSAp6ccdpWKFJZLZHCEkO3bT54ETp6EOvokLh2Lx7EbLjiEAEQhECcwGKmw0tnMwiwddcreR2D9dAR1Lo5mExuhdNzpzIMRADznEeWGvsO5hYeHCwCZbiEhIUIIIUJCQkRQUJDONhEREcLPz09YWVkJb29vsXDhQr0+09CG8OzbVw7VOGqU0pEQZU+tFmLZMiFKlJB/r+bmQkycKMTjx0pHpoezZ0VM1faiHv4SgBAWqlSxFAOFUKl0x01VqeQtNFTpiE2OoZ1/s6JEuSSEcRwbKgRPngjx119C/PCDECNGiLTGzcQ/9vXFL+grRuNr0QwRwgHxOqcszc3V8al4vVWimDM7XezdK99KR2hoxvmN5zwirdyef1VCZHWZwLAkJCTA0dER8fHxKF68uKKxPH0KuLgAjx8DBw8CgYGKhkOUo7g4YORIYN06+bxMGWDOHNmk2yBmvc5qIjwzM+DHH7F1xB8ISVmCh3CCk0MKfg21QuvEDXIkov8GJgAAlC0LzJsHdOPQ0vnNkM6/hobHxgjkNNGmPoSQJ9QTJ7Q1NYnHL+Kfi1b4W1TDKdTCMfjjBOrgCewzbW5jlY7aNQXqB1ogMFD+fvD2zsVEnRt4ziN6UW7Pv0x09BQaKoe0L1dO1kpzJmEyFps3y746V67I5w0bAl99BTRurGBQWRXg7u545FENo4+9hZ8xAADQwC8F6zZZoVy5/9bJzx8v9FKGdP41NDw2Bi6r84unpxytJacEITUVOHdOm9Q8O34G5048w98PPfA3amhvV+Gd5eb2tunwq6tCXX8z1K0L+PsDvr6v0FeS5zwiHbk9/xp59+TCt3atvO/Vi0kOGZcuXeTktv/7HzBzJnD4MNCkCdCqFTB1qnxcqLKYCE8A+DU2CONi5+IWykClEhg7Bpg5ywpWzzdlNzcHmjcv5ICJyGhkN9HmzZty+fr1GcnOgwfafjS3oq7iQnQSLl6xxIX08riAyriAd3EJFaFG1omFh1s6atQyQ40aKm1SU6mSef7mITznEeUJa3T08PixbLb29Kmct8TfX7FQiF5JbKxMbpYvB9LS5LKWLYEPPpCjtuW5SVt6OhARIW9qNeDsDLi6Am5u8vU7d+TVyEaN5JjYz11pPYhAjMNcRKERAKCy+SUsD/dBo6a8aqkkQzn/GqIicWxerElo1Ei2286qZiEvtQ4FUVORni7bhD1fk/OfFFjiOsriql01XKkQjMtXzHEh0Q0XUQkXUSnLJmcaJYunoWYtFWrUMkeNGnLo/urVASenVwuXiPTHGp0CsG2bTHIqVJAz0BMZK3d3YMkSYPJkWbuzfLmcd2fPHvn3PWwY0L8/ULq0Hm/6ww/A4MG5W7d0aeDuXQgAuxGMmZiMcLQEANjjMSZhFsamfw3b9N8BNNdz74hIbw8eyCQjJkYmGnXrygsVe/fKWSg1zM1lIqGhaQoG6N9M7FWalr1IrYa4/wCPLt7FrY2HcfNGVVxFO1yBN67CC1fgjSvwxi14QMAMeALgdOa3MVepUd7tCSr7mqFSLVtUrqJCpUoyoXFzs2BLDiIjwxodPXTtCmzaJH8cfv65YmEQ5burV+UcPMuWAY8eyWXm5kBwsBy0oHNn+ZsnW5aWGVVDuXAdnliDN/EjBuEiKgMALJCKEPyM6ZiKMrglV1y9mhPhKcxQzr+GSHNsnJzi4e5eHK6ustbf1RU6j11c5FV/JyfA0dEAu1a4uQG3b+dtW5Uq66GPNa8Bus3ENLJrWvb8Nl27AvHxeHrlNu5dfIh7l+Nx9+oT3L6RiluxKty6b4XYR3a49cQRt1JK4RY8kAybHEO2wVN44Sq8PVLhU78UKjcogcq1bVGpEuDjI09nRGTYOBhBvscgC6vkZDnYSq1aioRBVKCSkmRusXgxcOxYxnKVSjbVbNUKCAiQF3s9Pf/7TWJtDaSkvPR9H6AkDqMhDqAxwtAB0cioEi2GRAzATxiPL1EO13U3DA9nu3SFGcL511Bpjg0QDyB3x0alksmOJvEpWVL3ccmSgIMDUKyYvM/qcbFi+ZgsvUqSk0vCowweRxxF/O1nSLj9FPFxT5EwaRbiE1VIQHHEwxEJKI5HKIF7KPXcrTTuwfmlzcmy4mT2EO7qmzKZwRXtveaxC+7IObh4fiEyWkx08tnKlUC/fnLUlDNnOBABmb6LF+XgG7/9BpzOoolHyZKAd5kUlPs7DG6Igy2ewgbPoIaZ9sfLdZTFZVRALDx0tjVDOhrhIAZiOXriNxRDku6baybC0zSjIcUYwvnXUGmOzcGD8UhKKo7bt2U3tNu3dW937wIPH8p+nvnFzg6wtwdsbOS1hpzuraxk3zvNzdwcMEt5BrOF38MMap0bAKTBAukwz/L++ccpsMJT2OZws3vl/bVECkpZJcDZ9glcHZ7Co1QKPNzU8ChnAffyNvCo7ACPaiXg7mUFG8v/+ujcvPnyiTZ5fiEyWuyjk882bpT3PXowyaGioVIl4OOP5e3WLWD3bnkB9Ngxmew/fAg8fGiFaLyeq/erjPMIwCG0QDg6YjtK417WK2r+webN448QMgrVqwO5yQFTUmTT0AcP5O3hw8yPHz0CEhPl7fHjjMeam6Z7zJMn8vZqbACMe9U3yTVzpMFRlQBHs0QUT38IR8SjOBK09yXwSKc+p9TYEJQa9DpKlbGGg4MVVKpSuf4kfPONbBr3YtM6nl+IihTW6OTCkydAqVJyIILjxwE/v0IPgcigPHsGXLgAXPPrgutqD9yBC57BBs9gAxUEHBEPByTCA7dQAZdREZdQEo+yfrNSpYB7zyU9nAjPoCh9/jVkhX1shJDNpzVJT1KSfJ6cLP8nnz3LePzifUqK3F6tlsmSWg2o//cN1M+SX6jPMYOA6r96G3kzR7rO/fOPLZH6Qt3NkyzrdBy2r4Vt++Yyz4iIAFq0yHmHX7VpGSfaJDJZbLqWjzZtkn0ivbxkTTdrdIj+4+QkL0Xnhab5yKVL2Q9XS4pT+vxryIz+2NSsCfz9d8F+RlbNxNILsWkZJ9okMklsupaPNM3WunZlkkOk4/Rp+YNEX883H7GyYodgIiVERuYwnGIuPN80LLfNxMwLsWkZJ9okKtLyOi1gkZGaCmzZIh+zppvoBWXKyERFX56eWQ85S0SFx8lJjn+tjxeTD09PIDRU3sqUyfxadv/n3brJ1/TZhohIT6zRyUFkpOwcWrq0nBCaiF6QnPzyIaZHjZIT8QBySCo2HyEyHHFxOQ8xXbw48PbbQJcusiDMrqlply76NRPr1k3/bYiI9MBEJweaZmtduvDcS5St5GTZ3r5mTTnplLU18OOPQM+e/MchMnRxcXLIt6ZNM/rF1K0rk5pWrWTTr+f/j7NrCpaXZmJsWkZEBYiJzkuo1br9c4joJcqUkT+WiMj4ODkB//yjdBRERPmKfXRe4q+/ZG26gwMQHKx0NERERERElFtMdF5CU5vTsaNsiUNERERERMaBic5LsNkaEREREZFxYqKTjQsXgIsXAUtLoF07paMhIiIiIiJ9MNHJxrZt8j4oSI6sSURERERExoOJTjY0iU6nTsrGQURERERE+mOik4VHj+T8ZQATHSIiIiIiY8REJws7dwJpaYCvL1ChgtLREBERERGRvpjoZIHN1oiIiIiIjBsTnRekpwNhYfIxEx0iIiIiIuPEROcFhw8D9+8DJUoAjRopHQ0REREREeUFE50XaJqttWsn59AhIiIiIiLjk6dEZ8GCBfDx8YGNjQ38/f2xTzNEWRYiIiKgUqky3c6dO5fnoAsS++cQERknUy6biIhIf3onOmvXrsXo0aPx0UcfITo6Gk2bNkX79u1x7dq1l253/vx5xMbGam+VKlXKc9AF5do14PRpwMxM1ugQEZFxMOWyiYiI8kbvROfrr7/GO++8g0GDBqFq1aqYN28eypYti4ULF750OxcXF7i5uWlv5ubmeQ66oGgGIQgMBJydlY2FiIhyz5TLJiIiyhu9Ep2UlBQcO3YMbdq00Vnepk0bHDx48KXb+vn5wd3dHcHBwQgPD9c/0kKwY4e879BB2TiIiCj3TL1sIiKivLHQZ+V79+4hPT0drq6uOstdXV0RFxeX5Tbu7u5YsmQJ/P39kZycjF9++QXBwcGIiIhAs2bNstwmOTkZycnJ2ucJCQn6hJknKSnAnj3ycdu2Bf5xRESUT0y5bCIiorzTK9HRUKlUOs+FEJmWaVSpUgVVqlTRPg8MDMT169cxd+7cbAuTWbNmYfr06XkJLc+iooDERKB0acDPr1A/moiI8oEplk1ERJR3ejVdK1WqFMzNzTNdIbtz506mK2kvExAQgIsXL2b7+qRJkxAfH6+9Xb9+XZ8w8+SPP+R969ZyMAIiIjIOplw2ERFR3un1k97Kygr+/v7YtWuXzvJdu3ahkR6za0ZHR8Pd3T3b162trVG8eHGdW0HTJDocbY2IyLiYctlERER5p3fTtbFjx6Jfv36oV68eAgMDsWTJEly7dg1Dhw4FIK943bx5EytWrAAAzJs3D97e3qhevTpSUlKwcuVKhIaGIjQ0NH/35BXcuQMcPy4fv9CXlYiIjIAplk1ERPRq9E50evXqhfv37+PTTz9FbGwsatSogbCwMHh5eQEAYmNjdeYtSElJwbhx43Dz5k3Y2tqievXq2L59OzoY0NBmO3fK+zp1AD1aORARkYEwxbKJiIhejUoIIZQOIicJCQlwdHREfHx8gTQV6NcPWLkSmDABmD0739+eiMhoFfT515jx2BARKSO3598i3+1erc6o0WH/HCIiIiIi01DkE52TJ2UfnWLFAD36rBIRERERkQEr8omOZrS1Fi0AKytlYyEiIiIiovxR5BOdHTvkPZutERERERGZjiKd6CQlAQcPysccVpqIiIiIyHQU6URn3z4gNRXw8gIqVFA6GiIiIiIiyi9FOtHZvVveBwcDKpWysRARERERUf4p0onOnj3yPjhY2TiIiIiIiCh/FdlE58EDIDpaPm7RQtlYiIiIiIgofxXZRCc8HBACqFYNcHdXOhoiIiIiIspPRTbReb5/DhERERERmRYmOkx0iIiIiIhMTpFMdG7cAC5cAMzMgKAgpaMhIiIiIqL8ViQTHU1tTr16QIkSioZCREREREQFoEgnOmy2RkRERERkmopcoiNExvw5LVsqGwsRERERERWMIpfoXLgA3LwJWFsDjRsrHQ0RERERERWEIpfoaGpzGjUCbG2VjYWIiIiIiApGkUt0IiLkfYsWioZBREREREQFqEglOkIAkZHyMYeVJiIiIiIyXUUq0blwAbh9W/bPadBA6WiIiIiIiKigFKlER1ObExAA2NgoGwsRERERERWcIpnosNkaEREREZFpKzKJzvP9c5o1UzYWIiIiIiIqWEUm0YmJkfPnWFoCgYFKR0NERERERAWpyCQ6mtqc+vUBOztlYyEiIiIiooJV5BId9s8hIiIiIjJ9THSIiIiIiMjk5CnRWbBgAXx8fGBjYwN/f3/s27fvpetHRkbC398fNjY2KF++PBYtWpSnYPPq2jXgyhXA3Bxo1KhQP5qIiAqJsZVNRERUsPROdNauXYvRo0fjo48+QnR0NJo2bYr27dvj2rVrWa4fExODDh06oGnTpoiOjsbkyZMxcuRIhIaGvnLwuaWpzfH3BxwcCu1jiYiokBhj2URERAVLJYQQ+mzQsGFD1K1bFwsXLtQuq1q1Kl5//XXMmjUr0/oTJkzAli1bcPbsWe2yoUOH4uTJk4iKisrVZyYkJMDR0RHx8fEoXry4PuECAAYNApYuBcaPB774Qu/NiYiKrFc9/xYWYyybiIgob3J7/rXQ501TUlJw7NgxTJw4UWd5mzZtcPDgwSy3iYqKQps2bXSWtW3bFkuXLkVqaiosLS0zbZOcnIzk5GTt8/j4eAByp/IiPFze168P5PEtiIiKJM15V89rYoXKWMsmIiLKm9yWTXolOvfu3UN6ejpcXV11lru6uiIuLi7LbeLi4rJcPy0tDffu3YO7u3umbWbNmoXp06dnWl62bFl9ws2kZ89X2pyIqMi6f/8+HB0dlQ4jS8ZeNhERUd4kJia+tGzSK9HRUKlUOs+FEJmW5bR+Vss1Jk2ahLFjx2qfP3r0CF5eXrh27ZrBFrS5kZCQgLJly+L69etG3czBVPYDMJ194X4YHlPZl/j4eJQrVw5OTk5Kh5Ijlk3GwVT+Nwobj1ve8LjlnSEfOyEEEhMT4eHh8dL19Ep0SpUqBXNz80xXyO7cuZPpypiGm5tblutbWFjA2dk5y22sra1hbW2dabmjo6PBHei8KF68OPfDwJjKvnA/DI+p7IuZmeHORsCyyTiZyv9GYeNxyxset7wz1GOXmwtMepVcVlZW8Pf3x65du3SW79q1C42yGbc5MDAw0/o7d+5EvXr1smwDTUREpA+WTURElBW9L9GNHTsWP/74I5YtW4azZ89izJgxuHbtGoYOHQpAVu33799fu/7QoUNx9epVjB07FmfPnsWyZcuwdOlSjBs3Lv/2goiIijSWTURE9CK9++j06tUL9+/fx6efforY2FjUqFEDYWFh8PLyAgDExsbqzFvg4+ODsLAwjBkzBvPnz4eHhwe+/fZbdO/ePdefaW1tjalTp2bZZMCYcD8Mj6nsC/fD8JjKvhjLfrBsMh48bnnD45Y3PG55ZwrHTu95dIiIiIiIiAyd4fYuJSIiIiIiyiMmOkREREREZHKY6BARERERkclhokNERERERCbHIBOdzz//HI0aNYKdnR1KlCiRq22EEJg2bRo8PDxga2uL5s2b459//inYQHPh4cOH6NevHxwdHeHo6Ih+/frh0aNHL91mwIABUKlUOreAgIDCCfg/CxYsgI+PD2xsbODv7499+/a9dP3IyEj4+/vDxsYG5cuXx6JFiwop0pzpsy8RERGZjr1KpcK5c+cKMeLM9u7di86dO8PDwwMqlQqbNm3KcRtD/E703Q9D/T5mzZqF+vXrw8HBAS4uLnj99ddx/vz5HLcztO8kL/thqN9JftP3b3XDhg1o3bo1SpcujeLFiyMwMBB//PFH4QRrQPJyrtI4cOAALCwsUKdOnQKLz5Dl5dglJyfjo48+gpeXF6ytrVGhQgUsW7as4IM1IHk5bqtWrULt2rVhZ2cHd3d3DBw4EPfv3y/4YA2IqZRjOTHIRCclJQU9evTAsGHDcr3NF198ga+//hrff/89jhw5Ajc3N7Ru3RqJiYkFGGnO+vTpgxMnTmDHjh3YsWMHTpw4gX79+uW4Xbt27RAbG6u9hYWFFUK00tq1azF69Gh89NFHiI6ORtOmTdG+fXudoVmfFxMTgw4dOqBp06aIjo7G5MmTMXLkSISGhhZazNnRd180zp8/r3P8K1WqVEgRZy0pKQm1a9fG999/n6v1DfU70Xc/NAzt+4iMjMTw4cNx6NAh7Nq1C2lpaWjTpg2SkpKy3cYQv5O87IeGoX0n+U3fv9W9e/eidevWCAsLw7Fjx9CiRQt07twZ0dHRBRypYcnr/3h8fDz69++P4ODgAorM8OXl2PXs2RO7d+/G0qVLcf78eaxZswa+vr4FGKXh0fe47d+/H/3798c777yDf/75B+vWrcORI0cwaNCgAo7UsJhKOZYjYcCWL18uHB0dc1xPrVYLNzc3MXv2bO2yZ8+eCUdHR7Fo0aICjPDlzpw5IwCIQ4cOaZdFRUUJAOLcuXPZbhcSEiK6dOlSCBFmrUGDBmLo0KE6y3x9fcXEiROzXP/DDz8Uvr6+OsuGDBkiAgICCizG3NJ3X8LDwwUA8fDhw0KILm8AiI0bN750HUP+TjRysx/G8H0IIcSdO3cEABEZGZntOsbwneRmP4zlO8lPuflbzUq1atXE9OnT8z8gI6HPcevVq5f4+OOPxdSpU0Xt2rULNC5jkJtj9/vvvwtHR0dx//79wgnKCOTmuH355ZeifPnyOsu+/fZb4enpWYCRGT5TKcdeZJA1OvqKiYlBXFwc2rRpo11mbW2NoKAgHDx4ULG4oqKi4OjoiIYNG2qXBQQEwNHRMce4IiIi4OLigsqVK+Pdd9/FnTt3CjpcALI27dixYzrHEgDatGmTbcxRUVGZ1m/bti2OHj2K1NTUAos1J3nZFw0/Pz+4u7sjODgY4eHhBRlmgTDU7ySvDP37iI+PBwA4OTllu44xfCe52Q8NQ/9OlKZWq5GYmJirY1nULV++HJcvX8bUqVOVDsWobNmyBfXq1cMXX3yBMmXKoHLlyhg3bhyePn2qdGgGrVGjRrhx4wbCwsIghMDt27exfv16dOzYUenQFGUq5diLTCLRiYuLAwC4urrqLHd1ddW+poS4uDi4uLhkWu7i4vLSuNq3b49Vq1Zhz549+Oqrr3DkyBG0bNkSycnJBRkuAODevXtIT0/X61jGxcVluX5aWhru3btXYLHmJC/74u7ujiVLliA0NBQbNmxAlSpVEBwcjL179xZGyPnGUL8TfRnD9yGEwNixY9GkSRPUqFEj2/UM/TvJ7X4Yw3diCL766iskJSWhZ8+eSodi0C5evIiJEydi1apVsLCwUDoco/Lvv/9i//79+Pvvv7Fx40bMmzcP69evx/Dhw5UOzaA1atQIq1atQq9evWBlZQU3NzeUKFEC3333ndKhKcZUyrGsFNpZZdq0aZg+ffpL1zly5Ajq1auX589QqVQ6z4UQmZblh9zuS1Yx5SauXr16aR/XqFED9erVg5eXF7Zv345u3brlMWr96Hsss1o/q+VK0GdfqlSpgipVqmifBwYG4vr165g7dy6aNWtWoHHmN0P+TnLLGL6PESNG4NSpU9i/f3+O6xryd5Lb/TCG70Rpa9aswbRp07B58+YsL3aRlJ6ejj59+mD69OmoXLmy0uEYHbVaDZVKhVWrVsHR0REA8PXXX+ONN97A/PnzYWtrq3CEhunMmTMYOXIkPvnkE7Rt2xaxsbEYP348hg4diqVLlyodniJMpRzLSqElOiNGjEDv3r1fuo63t3ee3tvNzQ2AzDTd3d21y+/cuZMp88wPud2XU6dO4fbt25leu3v3rl5xubu7w8vLCxcvXtQ7Vn2VKlUK5ubmmWo8XnYs3dzcslzfwsICzs7OBRZrTvKyL1kJCAjAypUr8zu8AmWo30l+MKTv4/3338eWLVuwd+9eeHp6vnRdQ/5O9NmPrBjSd6K0tWvX4p133sG6devQqlUrpcMxaImJiTh69Ciio6MxYsQIAPLHuxACFhYW2LlzJ1q2bKlwlIbL3d0dZcqU0SY5AFC1alUIIXDjxg2TGyAkv8yaNQuNGzfG+PHjAQC1atWCvb09mjZtis8++0znd2RRYCrlWHYKLdEpVaoUSpUqVSDv7ePjAzc3N+zatQt+fn4AZP+MyMhIzJkzJ98/L7f7EhgYiPj4ePz1119o0KABAODw4cOIj49Ho0aNcv159+/fx/Xr1wvln8/Kygr+/v7YtWsXunbtql2+a9cudOnSJcttAgMDsXXrVp1lO3fuRL169WBpaVmg8b5MXvYlK9HR0UZ34jPU7yQ/GML3IYTA+++/j40bNyIiIgI+Pj45bmOI30le9iMrhvCdGII1a9bg7bffxpo1a4p8e//cKF68OE6fPq2zbMGCBdizZw/Wr1+f57/HoqJx48ZYt24dHj9+jGLFigEALly4ADMzszxdsCgqnjx5kqmZpLm5OYCM2omiwFTKsRwV9ugHuXH16lURHR0tpk+fLooVKyaio6NFdHS0SExM1K5TpUoVsWHDBu3z2bNnC0dHR7FhwwZx+vRp8eabbwp3d3eRkJCgxC5otWvXTtSqVUtERUWJqKgoUbNmTdGpUyeddZ7fl8TERPHBBx+IgwcPipiYGBEeHi4CAwNFmTJlCm1ffv31V2FpaSmWLl0qzpw5I0aPHi3s7e3FlStXhBBCTJw4UfTr10+7/r///ivs7OzEmDFjxJkzZ8TSpUuFpaWlWL9+faHE+zL67sv//vc/sXHjRnHhwgXx999/i4kTJwoAIjQ0VKldEELIvwvN/wEA8fXXX4vo6Ghx9epVIYTxfCf67oehfh/Dhg0Tjo6OIiIiQsTGxmpvT5480a5jDN9JXvbDUL+T/Kbv3+rq1auFhYWFmD9/vs6xfPTokVK7oAh9j9uLivKoa/oeu8TEROHp6SneeOMN8c8//4jIyEhRqVIlMWjQIKV2QRH6Hrfly5cLCwsLsWDBAnH58mWxf/9+Ua9ePdGgQQOldkERplKO5cQgE52QkBABINMtPDxcuw4AsXz5cu1ztVotpk6dKtzc3IS1tbVo1qyZOH36dOEH/4L79++Lvn37CgcHB+Hg4CD69u2baVjW5/flyZMnok2bNqJ06dLC0tJSlCtXToSEhIhr164Vatzz588XXl5ewsrKStStW1dnuMGQkBARFBSks35ERITw8/MTVlZWwtvbWyxcuLBQ430ZffZlzpw5okKFCsLGxkaULFlSNGnSRGzfvl2BqHVphvR98RYSEiKEMJ7vRN/9MNTvI6t9ePGcZAzfSV72w1C/k/ym799qUFDQS9cvKvJyrnpeUU508nLszp49K1q1aiVsbW2Fp6enGDt2rM4P1aIgL8ft22+/FdWqVRO2trbC3d1d9O3bV9y4caPwg1eQqZRjOVEJUYTq6YiIiIiIqEgwieGliYiIiIiInsdEh4iIiIiITA4THSIiIiIiMjlMdIiIiIiIyOQw0SEiIiIiIpPDRIeIiIiIiEwOEx0iIiIiIjI5THSIiIiIiMjkMNEhIiIiIiKTw0SHiIiIiIhMDhMdIiIiIiIyOUx0iIiIiIjI5PwfHEWjNw6wqdQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[10,3])\n",
    "plt.subplot(1,2,1)\n",
    "plt.scatter(x_list,y_list,c=\"r\")\n",
    "plt.plot(x_list,y_list,c=\"r\")\n",
    "plt.plot(x,f(x), c=\"b\")\n",
    "plt.xlim([-1,2.5])\n",
    "plt.ylim([0,3])\n",
    "plt.title(\"Gradient descent\")\n",
    "plt.subplot(1,2,2)\n",
    "plt.scatter(x_list,y_list,c=\"r\")\n",
    "plt.plot(x_list,y_list,c=\"r\")\n",
    "plt.plot(x,f(x), c=\"b\")\n",
    "plt.xlim([1.2,2.1])\n",
    "plt.ylim([0,3])\n",
    "plt.title(\"Gradient descent (zoomed in)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SLIM MSE with Cython\n",
    "\n",
    "#### Load the usual data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Movielens10M: Verifying data consistency...\n",
      "Movielens10M: Verifying data consistency... Passed!\n",
      "DataReader: current dataset is: Movielens10M\n",
      "\tNumber of items: 10681\n",
      "\tNumber of users: 69878\n",
      "\tNumber of interactions in URM_all: 10000054\n",
      "\tValue range in URM_all: 0.50-5.00\n",
      "\tInteraction density: 1.34E-02\n",
      "\tInteractions per user:\n",
      "\t\t Min: 2.00E+01\n",
      "\t\t Avg: 1.43E+02\n",
      "\t\t Max: 7.36E+03\n",
      "\tInteractions per item:\n",
      "\t\t Min: 0.00E+00\n",
      "\t\t Avg: 9.36E+02\n",
      "\t\t Max: 3.49E+04\n",
      "\tGini Index: 0.57\n",
      "\n",
      "\tICM name: ICM_year, Value range: 1.92E+03 / 2.01E+03, Num features: 1, feature occurrences: 10681, density 1.00E+00\n",
      "\tICM name: ICM_all, Value range: 1.00 / 69.00, Num features: 10126, feature occurrences: 128384, density 1.19E-03\n",
      "\tICM name: ICM_tags, Value range: 1.00 / 69.00, Num features: 10106, feature occurrences: 106820, density 9.90E-04\n",
      "\tICM name: ICM_genres, Value range: 1.00 / 1.00, Num features: 20, feature occurrences: 21564, density 1.01E-01\n",
      "\n",
      "\n",
      "Warning: 87 (0.12 %) of 69878 users have no sampled items\n"
     ]
    }
   ],
   "source": [
    "from Data_manager.split_functions.split_train_validation_random_holdout import split_train_in_two_percentage_global_sample\n",
    "from Data_manager.Movielens.Movielens10MReader import Movielens10MReader\n",
    "\n",
    "data_reader = Movielens10MReader()\n",
    "data_loaded = data_reader.load_data()\n",
    "\n",
    "URM_all = data_loaded.get_URM_all()\n",
    "\n",
    "URM_train, URM_test = split_train_in_two_percentage_global_sample(URM_all, train_percentage = 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<69878x10681 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 8000043 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What do we need for a SLIM MSE?\n",
    "\n",
    "* Loss function\n",
    "* Item-Item similarity matrix\n",
    "* Computing prediction\n",
    "* Update rule\n",
    "* Training loop and some patience\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_users, n_items = URM_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SLIM MSE is based on a dense item-item similarity matrix $ S \\in R^{I \\times I}$ \n",
    "\n",
    "### How to compute the predictions\n",
    "$ \\hat{r}_{ui} = \\sum_{j=0}^{I} r_{uj}S_{ji}$\n",
    "\n",
    "\n",
    "### The loss function we are interested in minimizing is\n",
    "\n",
    "$L = ||R - RS||_2 + \\alpha||S||_2$\n",
    "\n",
    "#### The gradient of this with respect to S is\n",
    "$\\frac{\\partial}{\\partial S} L = -2(R - RS)R + 2 \\alpha S$\n",
    "\n",
    "#### The update is going to be (we can remove the coefficients)\n",
    "$ S = S - \\frac{\\partial}{\\partial S}$, or \n",
    "\n",
    "$ S = S + l((R - RS)R - \\alpha S)$, with $l$ the learning rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WARNING:\n",
    "During the learning process we must account for the fact that a possible solution for SLIM EN is $S=I$, which is a trivial and useless solution. We want to avoid that the diagonal is 1. This can be done by either introducing a penalty in the loss function or by simply never using and updating the diagonal. In this notebook we will never learn the self-similarity (the similarity of an item with itself, i.e., the diagonal) by ensuring it is always zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: We create a dense similarity matrix, initialized as zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "item_item_S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: We sample an interaction and compute the prediction of the current SLIM model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "714222"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "\n",
    "sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "sample_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6512, 1369, 3.0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_id = URM_train_coo.row[sample_index]\n",
    "item_id = URM_train_coo.col[sample_index]\n",
    "true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "(user_id, item_id, true_rating)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "predicted_rating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The first predicted rating is zero, of course, the model is \"empty\"\n",
    "\n",
    "### Step 3: We compute the prediction error and update the item-item similarity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prediction_error = true_rating - predicted_rating\n",
    "prediction_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The error is positive, so we need to increase the prediction our model computes. Meaning, we have to increase the values in the item-item similarity matrix\n",
    "\n",
    "### Which item similarities do we modify? Only those we used to compute the prediction, i.e., only the items in the profile of the sampled user. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   4,    5,    7,    9,   14,   18,   19,   20,   21,   22,   24,\n",
       "         30,   31,   34,   36,   37,   38,   40,   43,   47,   51,   57,\n",
       "         58,   60,   63,   64,   68,   72,   74,   75,   81,   83,   84,\n",
       "         85,   86,   91,   92,   93,   95,  133,  148,  171,  175,  177,\n",
       "        178,  179,  187,  192,  196,  198,  201,  205,  206,  210,  217,\n",
       "        228,  248,  252,  259,  260,  271,  279,  280,  296,  300,  307,\n",
       "        327,  346,  376,  381,  399,  411,  413,  415,  416,  424,  425,\n",
       "        426,  443,  458,  461,  464,  488,  489,  504,  505,  509,  519,\n",
       "        521,  536,  537,  543,  546,  547,  554,  584,  585,  598,  603,\n",
       "        620,  621,  628,  638,  642,  653,  654,  658,  687,  694,  718,\n",
       "        755,  760,  788,  793,  798,  819,  822,  823,  826,  827,  849,\n",
       "        856,  864,  866,  867,  874,  902,  927,  942,  949,  958,  962,\n",
       "        968,  972,  992,  993, 1003, 1005, 1031, 1051, 1074, 1077, 1078,\n",
       "       1084, 1093, 1094, 1096, 1119, 1122, 1133, 1176, 1182, 1183, 1187,\n",
       "       1190, 1191, 1200, 1201, 1248, 1296, 1302, 1311, 1314, 1316, 1324,\n",
       "       1330, 1334, 1358, 1366, 1367, 1369, 1418, 1426, 1427, 1434, 1448,\n",
       "       1449, 1451, 1471, 1475, 1505, 1508, 1544, 1599, 1613, 1620, 1622,\n",
       "       1653, 1657, 1658, 1716, 1726, 1762, 1782, 1808, 1813, 1814, 1815,\n",
       "       1820, 1821, 1832, 1838, 1843, 1861, 1864, 1882, 1913, 1919, 1920,\n",
       "       1963, 2008, 2019, 2053, 2071, 2084, 2088, 2109, 2151, 2157, 2177,\n",
       "       2188, 2299, 2310, 2312, 2314, 2477, 2489, 2502, 2514, 2516, 2544,\n",
       "       2574, 2578, 2587, 2594, 2595, 2606, 2694, 2714, 2741, 2747, 2942,\n",
       "       2943, 2948, 2955, 2960, 3006, 3184, 3268, 3294, 3327, 3511, 3536,\n",
       "       3542, 3561, 3602, 3787, 3807, 3810, 3858, 4066, 4394, 4926, 5131,\n",
       "       5401, 6716, 7845, 8990], dtype=int32)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "items_in_user_profile = URM_train[user_id].indices\n",
    "items_in_user_profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5. , 3.5, 3. , 5. , 4. , 4.5, 4. , 3. , 3. , 4.5, 4.5, 1.5, 3.5,\n",
       "       5. , 3. , 4. , 4.5, 3.5, 5. , 4. , 4. , 4. , 0.5, 5. , 5. , 5. ,\n",
       "       5. , 3. , 4.5, 3.5, 3. , 4. , 3. , 4. , 3. , 3.5, 2.5, 3.5, 2.5,\n",
       "       5. , 4.5, 4. , 3.5, 5. , 5. , 4.5, 4. , 5. , 4.5, 5. , 4.5, 4.5,\n",
       "       5. , 4. , 5. , 4. , 3. , 1.5, 4.5, 3.5, 4. , 4. , 4. , 4.5, 3. ,\n",
       "       2.5, 1.5, 2.5, 3.5, 4.5, 4. , 4. , 4.5, 4.5, 4. , 3. , 4. , 2. ,\n",
       "       0.5, 2. , 4.5, 2.5, 4. , 4. , 4.5, 4. , 3. , 3. , 3. , 4. , 3. ,\n",
       "       2. , 3. , 3. , 2.5, 3.5, 3. , 3.5, 4.5, 3.5, 4. , 3.5, 4. , 2.5,\n",
       "       3.5, 4. , 2. , 4.5, 3.5, 3.5, 0.5, 2. , 4. , 4.5, 5. , 3. , 5. ,\n",
       "       2. , 4.5, 1. , 4. , 4. , 4.5, 3.5, 3. , 4. , 2.5, 3. , 3. , 4. ,\n",
       "       4.5, 4. , 4.5, 4. , 3.5, 3. , 4. , 4. , 1. , 2.5, 4. , 4.5, 4. ,\n",
       "       2.5, 3.5, 3.5, 3.5, 3.5, 4.5, 3. , 2.5, 3. , 4.5, 0.5, 4.5, 3. ,\n",
       "       1.5, 5. , 3. , 1. , 2. , 2. , 4.5, 4. , 3.5, 0.5, 2.5, 2. , 2.5,\n",
       "       4. , 3. , 3. , 4. , 3.5, 2.5, 3. , 3. , 4. , 4. , 3.5, 3.5, 2.5,\n",
       "       4. , 4. , 2.5, 2.5, 4.5, 4.5, 4.5, 2. , 4. , 3.5, 4. , 3. , 3. ,\n",
       "       2.5, 3. , 3.5, 4. , 4. , 4. , 3.5, 3.5, 4. , 3. , 3.5, 1.5, 3.5,\n",
       "       4. , 3. , 3.5, 4.5, 3. , 3.5, 4.5, 4.5, 3.5, 3. , 3. , 3. , 4.5,\n",
       "       3. , 3.5, 4.5, 4. , 3. , 3.5, 2.5, 4. , 3. , 2.5, 2. , 4.5, 2.5,\n",
       "       4. , 3.5, 4. , 5. , 3. , 4.5, 4. , 4.5, 3. , 1. , 1.5, 4. , 2.5,\n",
       "       0.5, 4.5, 4. , 3.5, 3. , 3. , 4.5, 3. , 4. , 2.5, 4.5, 3.5, 3. ,\n",
       "       1.5, 3.5, 5. , 3.5, 4. , 2. , 1. , 2.5])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratings_in_user_profile = URM_train[user_id].data\n",
    "ratings_in_user_profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Apply the update rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S[items_in_user_profile,item_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-4\n",
    "regularization_1 = 1e-3\n",
    "regularization_2 = 1e-3\n",
    "\n",
    "item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\n",
    "                                                               regularization_1 * (item_item_S[items_in_user_profile, item_id] > 0) -\n",
    "                                                               regularization_2 * item_item_S[items_in_user_profile,item_id])\n",
    "\n",
    "# Ensure diagonal is always zero\n",
    "item_item_S[item_id,item_id] = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0015 , 0.00105, 0.0009 , 0.0015 , 0.0012 , 0.00135, 0.0012 ,\n",
       "       0.0009 , 0.0009 , 0.00135, 0.00135, 0.00045, 0.00105, 0.0015 ,\n",
       "       0.0009 , 0.0012 , 0.00135, 0.00105, 0.0015 , 0.0012 , 0.0012 ,\n",
       "       0.0012 , 0.00015, 0.0015 , 0.0015 , 0.0015 , 0.0015 , 0.0009 ,\n",
       "       0.00135, 0.00105, 0.0009 , 0.0012 , 0.0009 , 0.0012 , 0.0009 ,\n",
       "       0.00105, 0.00075, 0.00105, 0.00075, 0.0015 , 0.00135, 0.0012 ,\n",
       "       0.00105, 0.0015 , 0.0015 , 0.00135, 0.0012 , 0.0015 , 0.00135,\n",
       "       0.0015 , 0.00135, 0.00135, 0.0015 , 0.0012 , 0.0015 , 0.0012 ,\n",
       "       0.0009 , 0.00045, 0.00135, 0.00105, 0.0012 , 0.0012 , 0.0012 ,\n",
       "       0.00135, 0.0009 , 0.00075, 0.00045, 0.00075, 0.00105, 0.00135,\n",
       "       0.0012 , 0.0012 , 0.00135, 0.00135, 0.0012 , 0.0009 , 0.0012 ,\n",
       "       0.0006 , 0.00015, 0.0006 , 0.00135, 0.00075, 0.0012 , 0.0012 ,\n",
       "       0.00135, 0.0012 , 0.0009 , 0.0009 , 0.0009 , 0.0012 , 0.0009 ,\n",
       "       0.0006 , 0.0009 , 0.0009 , 0.00075, 0.00105, 0.0009 , 0.00105,\n",
       "       0.00135, 0.00105, 0.0012 , 0.00105, 0.0012 , 0.00075, 0.00105,\n",
       "       0.0012 , 0.0006 , 0.00135, 0.00105, 0.00105, 0.00015, 0.0006 ,\n",
       "       0.0012 , 0.00135, 0.0015 , 0.0009 , 0.0015 , 0.0006 , 0.00135,\n",
       "       0.0003 , 0.0012 , 0.0012 , 0.00135, 0.00105, 0.0009 , 0.0012 ,\n",
       "       0.00075, 0.0009 , 0.0009 , 0.0012 , 0.00135, 0.0012 , 0.00135,\n",
       "       0.0012 , 0.00105, 0.0009 , 0.0012 , 0.0012 , 0.0003 , 0.00075,\n",
       "       0.0012 , 0.00135, 0.0012 , 0.00075, 0.00105, 0.00105, 0.00105,\n",
       "       0.00105, 0.00135, 0.0009 , 0.00075, 0.0009 , 0.00135, 0.00015,\n",
       "       0.00135, 0.0009 , 0.00045, 0.0015 , 0.0009 , 0.0003 , 0.0006 ,\n",
       "       0.0006 , 0.00135, 0.0012 , 0.00105, 0.00015, 0.00075, 0.0006 ,\n",
       "       0.00075, 0.0012 , 0.     , 0.0009 , 0.0012 , 0.00105, 0.00075,\n",
       "       0.0009 , 0.0009 , 0.0012 , 0.0012 , 0.00105, 0.00105, 0.00075,\n",
       "       0.0012 , 0.0012 , 0.00075, 0.00075, 0.00135, 0.00135, 0.00135,\n",
       "       0.0006 , 0.0012 , 0.00105, 0.0012 , 0.0009 , 0.0009 , 0.00075,\n",
       "       0.0009 , 0.00105, 0.0012 , 0.0012 , 0.0012 , 0.00105, 0.00105,\n",
       "       0.0012 , 0.0009 , 0.00105, 0.00045, 0.00105, 0.0012 , 0.0009 ,\n",
       "       0.00105, 0.00135, 0.0009 , 0.00105, 0.00135, 0.00135, 0.00105,\n",
       "       0.0009 , 0.0009 , 0.0009 , 0.00135, 0.0009 , 0.00105, 0.00135,\n",
       "       0.0012 , 0.0009 , 0.00105, 0.00075, 0.0012 , 0.0009 , 0.00075,\n",
       "       0.0006 , 0.00135, 0.00075, 0.0012 , 0.00105, 0.0012 , 0.0015 ,\n",
       "       0.0009 , 0.00135, 0.0012 , 0.00135, 0.0009 , 0.0003 , 0.00045,\n",
       "       0.0012 , 0.00075, 0.00015, 0.00135, 0.0012 , 0.00105, 0.0009 ,\n",
       "       0.0009 , 0.00135, 0.0009 , 0.0012 , 0.00075, 0.00135, 0.00105,\n",
       "       0.0009 , 0.00045, 0.00105, 0.0015 , 0.00105, 0.0012 , 0.0006 ,\n",
       "       0.0003 , 0.00075])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S[items_in_user_profile,item_id]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's check what the new prediction for the same user-item interaction would be"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.049250000000002"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "predicted_rating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The value is not zero anymore, but higher, we are moving in the right direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And now? Sample another interaction and repeat... a lot of times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's put all together in a training loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 5000 in 7.05 seconds, loss is 13.25. Samples per second 708.77\n",
      "Iteration 10000 in 12.87 seconds, loss is 13.35. Samples per second 776.99\n",
      "Iteration 15000 in 18.57 seconds, loss is 13.38. Samples per second 807.90\n",
      "Iteration 20000 in 24.19 seconds, loss is 13.36. Samples per second 826.83\n",
      "Iteration 25000 in 29.57 seconds, loss is 13.35. Samples per second 845.55\n",
      "Iteration 30000 in 35.08 seconds, loss is 13.33. Samples per second 855.31\n",
      "Iteration 35000 in 40.87 seconds, loss is 13.31. Samples per second 856.33\n",
      "Iteration 40000 in 46.45 seconds, loss is 13.28. Samples per second 861.06\n",
      "Iteration 45000 in 51.94 seconds, loss is 13.25. Samples per second 866.39\n",
      "Iteration 50000 in 57.29 seconds, loss is 13.23. Samples per second 872.79\n",
      "Iteration 55000 in 62.64 seconds, loss is 13.20. Samples per second 878.05\n",
      "Iteration 60000 in 68.33 seconds, loss is 13.16. Samples per second 878.10\n",
      "Iteration 65000 in 74.38 seconds, loss is 13.13. Samples per second 873.92\n",
      "Iteration 70000 in 80.36 seconds, loss is 13.11. Samples per second 871.12\n",
      "Iteration 75000 in 86.43 seconds, loss is 13.09. Samples per second 867.78\n",
      "Iteration 80000 in 92.43 seconds, loss is 13.06. Samples per second 865.53\n",
      "Iteration 85000 in 98.49 seconds, loss is 13.03. Samples per second 863.07\n",
      "Iteration 90000 in 104.66 seconds, loss is 13.01. Samples per second 859.96\n",
      "Iteration 95000 in 110.85 seconds, loss is 12.99. Samples per second 857.03\n",
      "Iteration 100000 in 117.00 seconds, loss is 12.96. Samples per second 854.73\n"
     ]
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "\n",
    "learning_rate = 1e-6\n",
    "regularization_1 = 1e-3\n",
    "regularization_2 = 1e-3\n",
    "\n",
    "loss = 0.0\n",
    "start_time = time.time()\n",
    "for sample_num in range(100000):\n",
    "    \n",
    "    # Randomly pick sample\n",
    "    sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "    user_id = URM_train_coo.row[sample_index]\n",
    "    item_id = URM_train_coo.col[sample_index]\n",
    "    true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "    # Compute prediction\n",
    "    predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "        \n",
    "    # Compute prediction error, or gradient\n",
    "    prediction_error = true_rating - predicted_rating\n",
    "    loss += prediction_error**2\n",
    "    \n",
    "    # Update model, in this case the similarity\n",
    "    items_in_user_profile = URM_train[user_id].indices\n",
    "    ratings_in_user_profile = URM_train[user_id].data\n",
    "    item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\n",
    "                                                                   regularization_1 * (item_item_S[items_in_user_profile, item_id] > 0) -\n",
    "                                                                   regularization_2 * item_item_S[items_in_user_profile,item_id])\n",
    "    \n",
    "    # Ensure diagonal is always zero\n",
    "    item_item_S[item_id,item_id] = 0.0\n",
    "\n",
    "    # Print some stats\n",
    "    if (sample_num +1)% 5000 == 0:\n",
    "        elapsed_time = time.time() - start_time\n",
    "        samples_per_second = (sample_num +1)/elapsed_time\n",
    "        print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What do we see? The loss generally goes down but may oscillate a bit.\n",
    "### How long do we train such a model?\n",
    "\n",
    "* An epoch: a complete loop over all the train data\n",
    "* Usually you train for multiple epochs. Depending on the algorithm and data 10s or 100s of epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_one_epoch(URM_train, item_item_S, learning_rate, regularization_1, regularization_2):\n",
    "    \n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "\n",
    "    loss = 0.0\n",
    "    start_time = time.time()\n",
    "    for sample_num in range(URM_train.nnz):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "        user_id = URM_train_coo.row[sample_index]\n",
    "        item_id = URM_train_coo.col[sample_index]\n",
    "        true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "        # Compute prediction\n",
    "        predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "\n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = true_rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Update model, in this case the similarity\n",
    "        items_in_user_profile = URM_train[user_id].indices\n",
    "        ratings_in_user_profile = URM_train[user_id].data\n",
    "        item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\n",
    "                                                                       regularization_1 * (item_item_S[items_in_user_profile, item_id] > 0) -\n",
    "                                                                       regularization_2 * item_item_S[items_in_user_profile,item_id])\n",
    "\n",
    "        # Ensure diagonal is always zero\n",
    "        item_item_S[item_id,item_id] = 0.0\n",
    "    \n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 5000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n",
    "         \n",
    "            # Stop training because this implementation is too slow\n",
    "            print(\"\\tStopping the epoch early because this implementation is too slow\")\n",
    "            return item_item_S\n",
    "            \n",
    "    return item_item_S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 5000 in 6.67 seconds, loss is 13.41. Samples per second 749.76\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 5.72 seconds, loss is 13.36. Samples per second 873.67\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 5.67 seconds, loss is 13.29. Samples per second 882.19\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 5.68 seconds, loss is 13.27. Samples per second 880.88\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 6.01 seconds, loss is 13.08. Samples per second 832.17\n",
      "\tStopping the epoch early because this implementation is too slow\n"
     ]
    }
   ],
   "source": [
    "n_items = URM_train.shape[1]\n",
    "learning_rate = 1e-6\n",
    "regularization_1 = 1e-3\n",
    "regularization_2 = 1e-3\n",
    "    \n",
    "item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "\n",
    "for n_epoch in range(5):\n",
    "    item_item_S = train_one_epoch(URM_train, item_item_S, learning_rate, regularization_1, regularization_2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 8.93654229e-06, 3.24833494e-05, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [1.59628922e-05, 0.00000000e+00, 1.37651520e-04, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [1.64268941e-05, 1.30798574e-04, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       ...,\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How do we use this similarity? As in a simple item-based KNN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's estimate the training time. Say we train for 10 epochs and we have 8M interactions in the train data..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 93596.90 seconds (1559.95 minutes, 26.00 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unacceptable!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's rewrite the loop with some smarter use of the data structures. In particular:\n",
    "* Use the indptr/indices data structures to get the seen items\n",
    "* Not much else we can do with this tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 5000 in 4.63 seconds, loss is 13.30. Samples per second 1079.06\n",
      "Iteration 10000 in 8.75 seconds, loss is 13.37. Samples per second 1142.22\n",
      "Iteration 15000 in 12.89 seconds, loss is 13.34. Samples per second 1163.44\n",
      "Iteration 20000 in 17.07 seconds, loss is 13.31. Samples per second 1171.58\n",
      "Iteration 25000 in 21.14 seconds, loss is 13.27. Samples per second 1182.75\n",
      "Iteration 30000 in 25.35 seconds, loss is 13.24. Samples per second 1183.37\n",
      "Iteration 35000 in 29.40 seconds, loss is 13.25. Samples per second 1190.34\n",
      "Iteration 40000 in 33.54 seconds, loss is 13.24. Samples per second 1192.77\n",
      "Iteration 45000 in 37.60 seconds, loss is 13.22. Samples per second 1196.81\n",
      "Iteration 50000 in 41.76 seconds, loss is 13.21. Samples per second 1197.34\n",
      "Iteration 55000 in 45.82 seconds, loss is 13.17. Samples per second 1200.23\n",
      "Iteration 60000 in 49.88 seconds, loss is 13.14. Samples per second 1202.98\n",
      "Iteration 65000 in 53.94 seconds, loss is 13.12. Samples per second 1205.14\n",
      "Iteration 70000 in 58.03 seconds, loss is 13.10. Samples per second 1206.37\n",
      "Iteration 75000 in 62.17 seconds, loss is 13.07. Samples per second 1206.44\n",
      "Iteration 80000 in 66.34 seconds, loss is 13.05. Samples per second 1205.93\n",
      "Iteration 85000 in 70.45 seconds, loss is 13.03. Samples per second 1206.57\n",
      "Iteration 90000 in 74.49 seconds, loss is 13.01. Samples per second 1208.15\n",
      "Iteration 95000 in 78.67 seconds, loss is 12.99. Samples per second 1207.53\n",
      "Iteration 100000 in 82.70 seconds, loss is 12.96. Samples per second 1209.13\n"
     ]
    }
   ],
   "source": [
    "URM_train_coo = URM_train.tocoo()\n",
    "item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "\n",
    "learning_rate = 1e-6\n",
    "regularization_1 = 1e-3\n",
    "regularization_2 = 1e-3\n",
    "\n",
    "loss = 0.0\n",
    "start_time = time.time()\n",
    "for sample_num in range(100000):\n",
    "    \n",
    "    # Randomly pick sample\n",
    "    sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "    user_id = URM_train_coo.row[sample_index]\n",
    "    item_id = URM_train_coo.col[sample_index]\n",
    "    true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "    # Compute prediction\n",
    "    predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "        \n",
    "    # Compute prediction error, or gradient\n",
    "    prediction_error = true_rating - predicted_rating\n",
    "    loss += prediction_error**2\n",
    "    \n",
    "    # Update model, in this case the similarity\n",
    "    items_in_user_profile = URM_train.indices[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "    ratings_in_user_profile = URM_train.data[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "    item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\n",
    "                                                                   regularization_1 * (item_item_S[items_in_user_profile, item_id] > 0) -\n",
    "                                                                   regularization_2 * item_item_S[items_in_user_profile,item_id])\n",
    "    \n",
    "    # Ensure diagonal is always zero\n",
    "    item_item_S[item_id,item_id] = 0.0\n",
    "    \n",
    "    # Print some stats\n",
    "    if (sample_num +1)% 5000 == 0:\n",
    "        elapsed_time = time.time() - start_time\n",
    "        samples_per_second = (sample_num +1)/elapsed_time\n",
    "        print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 66163.27 seconds (1102.72 minutes, 18.38 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Just as bad as before\n",
    "\n",
    "### Let's see what we can do with Cython\n",
    "### First step, just compile it. We do not have the data at compile time, so we put the loop in a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "def do_some_training(URM_train):\n",
    "    \n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    n_items = URM_train.shape[1]\n",
    "    \n",
    "    item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "\n",
    "    learning_rate = 1e-6\n",
    "    regularization_1 = 1e-3\n",
    "    regularization_2 = 1e-3\n",
    "\n",
    "    loss = 0.0\n",
    "    start_time = time.time()\n",
    "    for sample_num in range(100000):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "        user_id = URM_train_coo.row[sample_index]\n",
    "        item_id = URM_train_coo.col[sample_index]\n",
    "        true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "        # Compute prediction\n",
    "        predicted_rating = URM_train[user_id].dot(item_item_S[:,item_id])[0]\n",
    "\n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = true_rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Update model, in this case the similarity\n",
    "        items_in_user_profile = URM_train.indices[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "        ratings_in_user_profile = URM_train.data[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "        item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\n",
    "                                                                       regularization_1 * (item_item_S[items_in_user_profile,item_id] > 0) -\n",
    "                                                                       regularization_2 * item_item_S[items_in_user_profile,item_id])\n",
    "\n",
    "        # Ensure diagonal is always zero\n",
    "        item_item_S[item_id,item_id] = 0.0\n",
    "    \n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 5000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n",
    "    \n",
    "    return loss, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 5000 in 4.86 seconds, loss is 13.32. Samples per second 1028.64\n",
      "Iteration 10000 in 8.69 seconds, loss is 13.31. Samples per second 1150.86\n",
      "Iteration 15000 in 12.46 seconds, loss is 13.31. Samples per second 1203.60\n",
      "Iteration 20000 in 16.21 seconds, loss is 13.30. Samples per second 1233.82\n",
      "Iteration 25000 in 19.88 seconds, loss is 13.29. Samples per second 1257.41\n",
      "Iteration 30000 in 23.59 seconds, loss is 13.27. Samples per second 1271.66\n",
      "Iteration 35000 in 27.33 seconds, loss is 13.24. Samples per second 1280.75\n",
      "Iteration 40000 in 31.00 seconds, loss is 13.23. Samples per second 1290.48\n",
      "Iteration 45000 in 34.52 seconds, loss is 13.20. Samples per second 1303.46\n",
      "Iteration 50000 in 38.43 seconds, loss is 13.17. Samples per second 1301.08\n",
      "Iteration 55000 in 42.28 seconds, loss is 13.14. Samples per second 1300.82\n",
      "Iteration 60000 in 46.02 seconds, loss is 13.13. Samples per second 1303.87\n",
      "Iteration 65000 in 49.75 seconds, loss is 13.08. Samples per second 1306.60\n",
      "Iteration 70000 in 53.36 seconds, loss is 13.05. Samples per second 1311.89\n",
      "Iteration 75000 in 56.89 seconds, loss is 13.03. Samples per second 1318.36\n",
      "Iteration 80000 in 60.55 seconds, loss is 13.01. Samples per second 1321.20\n",
      "Iteration 85000 in 64.22 seconds, loss is 12.98. Samples per second 1323.65\n",
      "Iteration 90000 in 67.76 seconds, loss is 12.96. Samples per second 1328.18\n",
      "Iteration 95000 in 71.69 seconds, loss is 12.93. Samples per second 1325.08\n",
      "Iteration 100000 in 75.49 seconds, loss is 12.91. Samples per second 1324.72\n"
     ]
    }
   ],
   "source": [
    "loss, samples_per_second = do_some_training(URM_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 60390.28 seconds (1006.50 minutes, 16.78 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Still far too time consuming\n",
    "### The compiler is just porting in C all operations that the python interpreter would have to perform, dynamic typing included. Have a look at the html reports in the Cython_examples folder\n",
    "\n",
    "### Now try to add some types: If you use a variable only as a C object, use primitive tipes\n",
    "\n",
    "* cdef int namevar\n",
    "* cdef double namevar\n",
    "* cdef float namevar\n",
    "* cdef double[:] singledimensionarray\n",
    "* cdef double[:,:] bidimensionalmatrix\n",
    "\n",
    "\n",
    "### We now use types for all main variables\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "def do_some_training(URM_train):\n",
    "\n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    n_items = URM_train.shape[1]\n",
    "\n",
    "    cdef double[:,:] item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "    cdef double learning_rate = 1e-6, regularization_1 = 1e-3, regularization_2 = 1e-3\n",
    "    cdef double loss = 0.0\n",
    "    cdef long start_time = time.time()\n",
    "    cdef double true_rating, predicted_rating, prediction_error, profile_rating\n",
    "    cdef int[:] items_in_user_profile\n",
    "    cdef double[:] ratings_in_user_profile\n",
    "    cdef int index, sample_num, user_id, item_id, profile_item_id\n",
    "\n",
    "    for sample_num in range(100000):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        sample_index = np.random.randint(URM_train_coo.nnz)\n",
    "\n",
    "        user_id = URM_train_coo.row[sample_index]\n",
    "        item_id = URM_train_coo.col[sample_index]\n",
    "        true_rating = URM_train_coo.data[sample_index]\n",
    "\n",
    "        # Compute prediction\n",
    "        items_in_user_profile = URM_train.indices[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "        ratings_in_user_profile = URM_train.data[URM_train.indptr[user_id]:URM_train.indptr[user_id+1]]\n",
    "        predicted_rating = 0.0\n",
    "\n",
    "        for index in range(len(items_in_user_profile)):\n",
    "            profile_item_id = items_in_user_profile[index]\n",
    "            profile_rating = ratings_in_user_profile[index]\n",
    "            predicted_rating += profile_rating * item_item_S[profile_item_id,item_id]\n",
    "\n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = true_rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Update model, in this case the similarity\n",
    "        for index in range(len(items_in_user_profile)):\n",
    "            profile_item_id = items_in_user_profile[index]\n",
    "            profile_rating = ratings_in_user_profile[index]\n",
    "            item_item_S[profile_item_id,item_id] += learning_rate * (prediction_error * profile_rating -\n",
    "                                                                     regularization_1 * (item_item_S[profile_item_id,item_id] > 0) -\n",
    "                                                                     regularization_2 * item_item_S[profile_item_id,item_id])\n",
    "\n",
    "        # Ensure diagonal is always zero\n",
    "        item_item_S[item_id,item_id] = 0.0\n",
    "    \n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 5000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n",
    "\n",
    "    return loss, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 5000 in 1.20 seconds, loss is 13.31. Samples per second 4167.44\n",
      "Iteration 10000 in 1.62 seconds, loss is 13.35. Samples per second 6179.08\n",
      "Iteration 15000 in 1.88 seconds, loss is 13.37. Samples per second 7992.60\n",
      "Iteration 20000 in 2.11 seconds, loss is 13.30. Samples per second 9465.20\n",
      "Iteration 25000 in 2.32 seconds, loss is 13.23. Samples per second 10765.29\n",
      "Iteration 30000 in 2.57 seconds, loss is 13.23. Samples per second 11685.30\n",
      "Iteration 35000 in 2.77 seconds, loss is 13.21. Samples per second 12641.41\n",
      "Iteration 40000 in 2.95 seconds, loss is 13.20. Samples per second 13562.21\n",
      "Iteration 45000 in 3.12 seconds, loss is 13.18. Samples per second 14407.45\n",
      "Iteration 50000 in 3.30 seconds, loss is 13.15. Samples per second 15160.34\n",
      "Iteration 55000 in 3.47 seconds, loss is 13.13. Samples per second 15827.43\n",
      "Iteration 60000 in 3.66 seconds, loss is 13.11. Samples per second 16391.67\n",
      "Iteration 65000 in 3.82 seconds, loss is 13.08. Samples per second 17029.99\n",
      "Iteration 70000 in 3.98 seconds, loss is 13.07. Samples per second 17604.08\n",
      "Iteration 75000 in 4.14 seconds, loss is 13.03. Samples per second 18135.68\n",
      "Iteration 80000 in 4.29 seconds, loss is 13.00. Samples per second 18634.18\n",
      "Iteration 85000 in 4.45 seconds, loss is 12.98. Samples per second 19095.09\n",
      "Iteration 90000 in 4.61 seconds, loss is 12.94. Samples per second 19510.13\n",
      "Iteration 95000 in 4.77 seconds, loss is 12.92. Samples per second 19911.86\n",
      "Iteration 100000 in 4.94 seconds, loss is 12.89. Samples per second 20256.18\n"
     ]
    }
   ],
   "source": [
    "loss, samples_per_second = do_some_training(URM_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 3949.41 seconds (65.82 minutes, 1.10 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This is why we use it for machine learning algorithms...\n",
    "\n",
    "### Some operations are still done with sparse matrices, those cannot be correctly optimized because the compiler does not know how what is the type of the data.\n",
    "\n",
    "### To address this, we create typed arrays in which we put the URM_train data\n",
    "For example, this operation: user_id = URM_train_coo.row[sample_index]\n",
    "\n",
    "Becomes:\n",
    "\n",
    "cdef int user_id\n",
    "\n",
    "cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "\n",
    "user_id = URM_train_coo_row[sample_index]\n",
    "\n",
    "\n",
    "### We can also skip the creation of the items_in_user_profile array and replace the np.random call with the faster native C function rand()\n",
    "\n",
    "## To obtain a more reliable speed estimate we increase the number of samples and the print step by a factor of 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "from libc.stdlib cimport rand, srand, RAND_MAX\n",
    "\n",
    "def do_some_training(URM_train):\n",
    "\n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    cdef int n_items = URM_train.shape[1]\n",
    "    cdef int n_interactions = URM_train.nnz\n",
    "    cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "    cdef int[:] URM_train_coo_col = URM_train_coo.col\n",
    "    cdef double[:] URM_train_coo_data = URM_train_coo.data\n",
    "    cdef int[:] URM_train_indices = URM_train.indices\n",
    "    cdef int[:] URM_train_indptr = URM_train.indptr\n",
    "    cdef double[:] URM_train_data = URM_train.data\n",
    "\n",
    "    cdef double[:,:] item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "    cdef double learning_rate = 1e-6, regularization_1 = 1e-3, regularization_2 = 1e-3\n",
    "    cdef double loss = 0.0\n",
    "    cdef long start_time = time.time()\n",
    "    cdef double true_rating, predicted_rating, prediction_error, profile_rating\n",
    "    cdef int start_profile, end_profile\n",
    "    cdef int index, sample_num, user_id, item_id, profile_item_id\n",
    "\n",
    "    for sample_num in range(1000000):\n",
    "\n",
    "        # Randomly pick sample\n",
    "        index = rand() % n_interactions\n",
    "\n",
    "        user_id = URM_train_coo_row[index]\n",
    "        item_id = URM_train_coo_col[index]\n",
    "        true_rating = URM_train_coo_data[index]\n",
    "\n",
    "        # Compute prediction\n",
    "        start_profile = URM_train_indptr[user_id]\n",
    "        end_profile = URM_train_indptr[user_id+1]\n",
    "        predicted_rating = 0.0\n",
    "\n",
    "        for index in range(start_profile, end_profile):\n",
    "            profile_item_id = URM_train_indices[index]\n",
    "            profile_rating = URM_train_data[index]\n",
    "            predicted_rating += item_item_S[profile_item_id,item_id] * profile_rating\n",
    "\n",
    "        # Compute prediction error, or gradient\n",
    "        prediction_error = true_rating - predicted_rating\n",
    "        loss += prediction_error**2\n",
    "\n",
    "        # Update model, in this case the similarity\n",
    "        for index in range(start_profile, end_profile):\n",
    "            profile_item_id = URM_train_indices[index]\n",
    "            profile_rating = URM_train_data[index]\n",
    "            item_item_S[profile_item_id,item_id] += learning_rate * (prediction_error * profile_rating - \n",
    "                                                                     regularization_1 * (item_item_S[profile_item_id,item_id] > 0) -\n",
    "                                                                     regularization_2 * item_item_S[profile_item_id,item_id])\n",
    "\n",
    "        # Ensure diagonal is always zero\n",
    "        item_item_S[item_id,item_id] = 0.0\n",
    "        \n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 50000 == 0:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds, loss is {:.2f}. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, loss/(sample_num +1), samples_per_second))\n",
    "\n",
    "    return loss, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 50000 in 3.12 seconds, loss is 13.19. Samples per second 16024.97\n",
      "Iteration 100000 in 4.21 seconds, loss is 12.93. Samples per second 23758.03\n",
      "Iteration 150000 in 5.26 seconds, loss is 12.70. Samples per second 28526.27\n",
      "Iteration 200000 in 6.29 seconds, loss is 12.47. Samples per second 31805.45\n",
      "Iteration 250000 in 7.29 seconds, loss is 12.25. Samples per second 34310.53\n",
      "Iteration 300000 in 8.28 seconds, loss is 12.04. Samples per second 36221.45\n",
      "Iteration 350000 in 9.25 seconds, loss is 11.86. Samples per second 37856.85\n",
      "Iteration 400000 in 10.23 seconds, loss is 11.67. Samples per second 39112.58\n",
      "Iteration 450000 in 11.19 seconds, loss is 11.50. Samples per second 40225.88\n",
      "Iteration 500000 in 12.17 seconds, loss is 11.33. Samples per second 41085.30\n",
      "Iteration 550000 in 13.20 seconds, loss is 11.18. Samples per second 41661.56\n",
      "Iteration 600000 in 14.15 seconds, loss is 11.02. Samples per second 42389.78\n",
      "Iteration 650000 in 15.11 seconds, loss is 10.88. Samples per second 43027.64\n",
      "Iteration 700000 in 16.05 seconds, loss is 10.74. Samples per second 43618.52\n",
      "Iteration 750000 in 17.15 seconds, loss is 10.61. Samples per second 43722.56\n",
      "Iteration 800000 in 18.13 seconds, loss is 10.48. Samples per second 44128.69\n",
      "Iteration 850000 in 19.09 seconds, loss is 10.36. Samples per second 44526.81\n",
      "Iteration 900000 in 20.06 seconds, loss is 10.25. Samples per second 44867.47\n",
      "Iteration 950000 in 21.03 seconds, loss is 10.14. Samples per second 45175.27\n",
      "Iteration 1000000 in 22.12 seconds, loss is 10.03. Samples per second 45208.94\n"
     ]
    }
   ],
   "source": [
    "loss, samples_per_second = do_some_training(URM_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 1769.56 seconds (29.49 minutes, 0.49 hours)\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds ({:.2f} minutes, {:.2f} hours)\".format(estimated_seconds, estimated_seconds/60, estimated_seconds/3600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As the source code gets less readable due to the addition of types and native C functions, it also gets remarkably faster\n",
    "\n",
    "### We started from a naive python implementation which took 30 hours (1k samples per second) and we now have a Cython one with takes only a few minutes (50k samples per second)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "from libc.stdlib cimport rand, srand, RAND_MAX\n",
    "\n",
    "def train_multiple_epochs(URM_train, learning_rate_input, regularization_1_input, regularization_2_input, n_epochs):\n",
    "\n",
    "    URM_train_coo = URM_train.tocoo()\n",
    "    cdef int n_items = URM_train.shape[1]\n",
    "    cdef int n_interactions = URM_train.nnz\n",
    "    cdef int[:] URM_train_coo_row = URM_train_coo.row\n",
    "    cdef int[:] URM_train_coo_col = URM_train_coo.col\n",
    "    cdef double[:] URM_train_coo_data = URM_train_coo.data\n",
    "    cdef int[:] URM_train_indices = URM_train.indices\n",
    "    cdef int[:] URM_train_indptr = URM_train.indptr\n",
    "    cdef double[:] URM_train_data = URM_train.data\n",
    "\n",
    "    cdef double[:,:] item_item_S = np.zeros((n_items, n_items), dtype = float)\n",
    "    cdef double learning_rate = learning_rate_input\n",
    "    cdef double regularization_1 = regularization_1_input\n",
    "    cdef double regularization_2 = regularization_2_input\n",
    "    cdef double loss = 0.0\n",
    "    cdef long start_time\n",
    "    cdef double true_rating, predicted_rating, prediction_error, profile_rating\n",
    "    cdef int start_profile, end_profile\n",
    "    cdef int index, sample_num, user_id, item_id, profile_item_id\n",
    "    \n",
    "    for n_epoch in range(n_epochs):\n",
    "        \n",
    "        loss = 0.0\n",
    "        start_time = time.time()\n",
    "        \n",
    "        for sample_num in range(n_interactions):\n",
    "\n",
    "            # Randomly pick sample\n",
    "            index = rand() % n_interactions\n",
    "\n",
    "            user_id = URM_train_coo_row[index]\n",
    "            item_id = URM_train_coo_col[index]\n",
    "            true_rating = URM_train_coo_data[index]\n",
    "\n",
    "            # Compute prediction\n",
    "            start_profile = URM_train_indptr[user_id]\n",
    "            end_profile = URM_train_indptr[user_id+1]\n",
    "            predicted_rating = 0.0\n",
    "\n",
    "            for index in range(start_profile, end_profile):\n",
    "                profile_item_id = URM_train_indices[index]\n",
    "                profile_rating = URM_train_data[index]\n",
    "                predicted_rating += item_item_S[profile_item_id,item_id] * profile_rating\n",
    "\n",
    "            # Compute prediction error, or gradient\n",
    "            prediction_error = true_rating - predicted_rating\n",
    "            loss += prediction_error**2\n",
    "\n",
    "            # Update model, in this case the similarity\n",
    "            for index in range(start_profile, end_profile):\n",
    "                profile_item_id = URM_train_indices[index]\n",
    "                profile_rating = URM_train_data[index]\n",
    "                item_item_S[profile_item_id,item_id] += learning_rate * (prediction_error * profile_rating - \n",
    "                                                                         regularization_1 * (item_item_S[profile_item_id,item_id] > 0) -\n",
    "                                                                         regularization_2 * item_item_S[profile_item_id,item_id])\n",
    "\n",
    "            # Ensure diagonal is always zero\n",
    "            item_item_S[item_id,item_id] = 0.0\n",
    "        \n",
    "#             if sample_num % 1000000 == 0:\n",
    "#                 print(\"Epoch {}: {:.2f}%\".format(n_epoch+1, sample_num/n_interactions*100))\n",
    "            \n",
    "            \n",
    "        elapsed_time = time.time() - start_time\n",
    "        samples_per_second = (sample_num+1)/elapsed_time\n",
    "     \n",
    "        print(\"Epoch {} complete in in {:.2f} seconds, loss is {:.3E}. Samples per second {:.2f}\".format(n_epoch+1, time.time() - start_time, loss/(sample_num+1), samples_per_second))\n",
    "\n",
    "    return np.array(item_item_S), loss/(sample_num+1), samples_per_second"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now let's compare the loss function if we use no regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1 complete in in 160.86 seconds, loss is 5.754E+00. Samples per second 49733.88\n",
      "Epoch 2 complete in in 177.63 seconds, loss is 3.707E+00. Samples per second 45037.10\n",
      "Epoch 3 complete in in 175.62 seconds, loss is 3.240E+00. Samples per second 45552.25\n",
      "Epoch 4 complete in in 178.31 seconds, loss is 2.986E+00. Samples per second 44866.18\n",
      "Epoch 5 complete in in 177.54 seconds, loss is 2.824E+00. Samples per second 45060.94\n",
      "Epoch 6 complete in in 181.85 seconds, loss is 2.701E+00. Samples per second 43993.20\n",
      "Epoch 7 complete in in 177.91 seconds, loss is 2.608E+00. Samples per second 44967.43\n",
      "Epoch 8 complete in in 182.84 seconds, loss is 2.536E+00. Samples per second 43755.38\n",
      "Epoch 9 complete in in 186.96 seconds, loss is 2.475E+00. Samples per second 42790.43\n",
      "Epoch 10 complete in in 172.62 seconds, loss is 2.420E+00. Samples per second 46345.60\n"
     ]
    }
   ],
   "source": [
    "n_items = URM_train.shape[1]\n",
    "learning_rate = 1e-6\n",
    "regularization_1 = 1e-3\n",
    "regularization_2 = 1e-3\n",
    "    \n",
    "item_item_S, loss, samples_per_second = train_multiple_epochs(URM_train, learning_rate, regularization_1, regularization_2, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1 complete in in 182.39 seconds, loss is 5.753E+00. Samples per second 43862.25\n",
      "Epoch 2 complete in in 183.98 seconds, loss is 3.707E+00. Samples per second 43482.14\n",
      "Epoch 3 complete in in 181.38 seconds, loss is 3.244E+00. Samples per second 44107.56\n",
      "Epoch 4 complete in in 178.94 seconds, loss is 2.987E+00. Samples per second 44709.07\n",
      "Epoch 5 complete in in 179.09 seconds, loss is 2.821E+00. Samples per second 44670.98\n",
      "Epoch 6 complete in in 181.81 seconds, loss is 2.701E+00. Samples per second 44003.44\n",
      "Epoch 7 complete in in 180.49 seconds, loss is 2.607E+00. Samples per second 44323.22\n",
      "Epoch 8 complete in in 176.34 seconds, loss is 2.536E+00. Samples per second 45368.11\n",
      "Epoch 9 complete in in 179.19 seconds, loss is 2.472E+00. Samples per second 44644.40\n",
      "Epoch 10 complete in in 180.35 seconds, loss is 2.419E+00. Samples per second 44358.80\n"
     ]
    }
   ],
   "source": [
    "item_item_S, loss, samples_per_second = train_multiple_epochs(URM_train, learning_rate, 0.0, 0.0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## How to use Cython outside a notebook\n",
    "\n",
    "### Step1: Create a .pyx file and write your code\n",
    "\n",
    "### Step2: Create a compilation script \"compileCython.py\" with the following content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# This code will not run in a notebook cell\n",
    "\n",
    "try:\n",
    "    from setuptools import setup\n",
    "    from setuptools import Extension\n",
    "except ImportError:\n",
    "    from distutils.core import setup\n",
    "    from distutils.extension import Extension\n",
    "\n",
    "\n",
    "from Cython.Distutils import build_ext\n",
    "import numpy\n",
    "import sys\n",
    "import re\n",
    "\n",
    "\n",
    "if len(sys.argv) != 4:\n",
    "    raise ValueError(\"Wrong number of paramethers received. Expected 4, got {}\".format(sys.argv))\n",
    "\n",
    "\n",
    "# Get the name of the file to compile\n",
    "fileToCompile = sys.argv[1]\n",
    "\n",
    "# Remove the argument from sys argv in order for it to contain only what setup needs\n",
    "del sys.argv[1]\n",
    "\n",
    "extensionName = re.sub(\"\\.pyx\", \"\", fileToCompile)\n",
    "\n",
    "\n",
    "ext_modules = Extension(extensionName,\n",
    "                [fileToCompile],\n",
    "                extra_compile_args=['-O3'],\n",
    "                include_dirs=[numpy.get_include(),],\n",
    "                )\n",
    "\n",
    "setup(\n",
    "    cmdclass={'build_ext': build_ext},\n",
    "    ext_modules=[ext_modules]\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step3: Compile your code with the following command \n",
    "\n",
    "python compileCython.py Cython_examples\\SLIM_MSE_fastest.pyx build_ext --inplace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step4: Generate cython report and look for \"yellow lines\". The report is an .html file which represents how many operations are necessary to translate each python operation in cython code. If a line is white, it has a direct C translation. If it is yellow it will require many indirect steps that will slow down execution. Some of those steps may be inevitable, some may be removed via static typing.\n",
    "\n",
    "### IMPORTANT: white does not mean fast!! If a system call is involved that part might be slow anyway.\n",
    "\n",
    "cython -a Cython_examples\\SLIM_MSE_fastest.pyx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step5: Add static types and C functions to remove \"yellow\" lines.\n",
    "\n",
    "#### If you use a variable only as a C object, use primitive tipes \n",
    "cdef int namevar\n",
    "\n",
    "cdef double namevar\n",
    "\n",
    "cdef float namevar\n",
    "\n",
    "#### If you call a function only within C code, use a specific declaration \"cdef\"\n",
    "\n",
    "cdef function_name(self, int param1, double param2):\n",
    "...\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step6: Iterate step 4 and 5 until you are satisfied with how clean your code is, then compile. An example of non optimized code can be found in the source folder of this notebook with the _SLOW suffix\n",
    "\n",
    "## Step7: the compilation generates a file whose name is something like \"SLIM_MSE_fastest.cp36-win_amd64.pyd\" and tells you the source file, the architecture it is compiled for and the OS\n",
    "\n",
    "## Step8: Import and use the compiled file as if it were any python object, function or class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cython_examples.SLIM_MSE_fastest import train_multiple_epochs\n",
    "\n",
    "train_multiple_epochs(URM_train, 1e-3, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A few warnings on ML algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - Why do we bother with KNNs if we have ML?\n",
    "Because sometimes ML algorithms work better than heuristic ones, sometimes they do not\n",
    "\n",
    "### - ML algorithms are always best because they learn from the data\n",
    "Not really... Yes they learn from the data but the data is sometimes too noisy, too sparse and does not yeld good results. There is plenty of examples of cases where ML algorithms are not the best choice.\n",
    "\n",
    "### - We should use this complex model because it can in theory approximate any function!!\n",
    "Theory is important but... does it work in practice? Often complex modes are veeeery difficult to train and you need to use lots of tricks: adaptive gradients, data normalization, careful initialization and crafted batches...\n",
    "\n",
    "### - I have trained my model for 2 epochs but the result is not great\n",
    "Have you just used the default learning rate or have you optimized it? Why did you stop after 2 epochs? You may need 100s of epochs\n",
    "\n",
    "### - If I select a high learning rate (maybe 1e-3) after 5 epochs I get a result wich is not very good, if I use a lower learning rate (maybe 1e-6) the result is much worse\n",
    "Of course, the lower the learning rate the slower the training process but the better the solution you may find. Again, you may need 100s of epochs\n",
    "\n",
    "### - Training and optimizing the hyperparameters of this ML model takes several hours, what am I doing wrong?\n",
    "Probably nothing, ML is computationally expensive and takes time... A few hours is a normal timespan. Sometimes the end result is still not satisfactory."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RecSysFramework",
   "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.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}