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 prime 2? True\n",
      "Is prime 3? True\n",
      "Is prime 5? True\n",
      "Is prime 15? False\n",
      "Is prime 20? False\n"
     ]
    }
   ],
   "source": [
    "print(\"Is prime 2? {}\".format(isPrime(2)))\n",
    "print(\"Is prime 3? {}\".format(isPrime(3)))\n",
    "print(\"Is prime 5? {}\".format(isPrime(5)))\n",
    "print(\"Is prime 15? {}\".format(isPrime(15)))\n",
    "print(\"Is prime 20? {}\".format(isPrime(20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Is Prime 80000023? True, time required 7.22 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": [
    "#### 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 3.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 tiping"
   ]
  },
  {
   "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 Prime 80000023? True, time required 0.68 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": [
    "#### Cython code with two tipe declaration, 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 opertions 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": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApEUlEQVR4nO3dd1hUZ97/8fd3KCIgAgI27GI3NuwxMaaYuibZ9J5s1iRPspvsJtuyfZ/8su1JdjdNY8omm2STrDHFTTTGFHsFFRURREBBihSlCgJz//5g8jwsARlw4MzM+b6uay4H5jB+OODHM/fc5z5ijEEppZT/c1gdQCmlVPfQwldKKZvQwldKKZvQwldKKZvQwldKKZvQwldKKZtot/BFJEREdohIioikishvW9lGROQZEckUkb0iMrVr4iqllOqsQDe2qQMWGGOqRCQI2CQiq40x25ptcxmQ4LrNBJa4/lRKKeUl2j3CN02qXB8GuW4tz9ZaBPzDte02IFJE+ns2qlJKqbPhzhE+IhIAJAMjgeeNMdtbbDIQyG32cZ7rcwUtnmcxsBggLCxs2pgxYzoZWymlrGUMpBaUEx0WzIDePbvt701OTi4xxsR25mvdKnxjTCMwWUQigQ9EZIIxZn+zTaS1L2vleZYBywASExNNUlJSxxMrpZQX2H30BNe8sIUlt07lsondN6AhIkc6+7UdmqVjjDkJrAMubfFQHjCo2cfxQH5nQymllLdLyjkBwLShURYncZ87s3RiXUf2iEhP4CLgYIvNVgJ3uGbrzALKjTEFKKWUn9qZU8bQPqHE9QqxOorb3BnS6Q+87hrHdwD/MsZ8LCL3AxhjlgKrgMuBTKAGuLuL8iqllOWMMSQdOcGCMXFWR+mQdgvfGLMXmNLK55c2u2+ABz0bTSmlvFNWSTVl1aeZ7kPDOaBn2iqlVIcl5ZQBkDg02uIkHaOFr5RSHbQz5wTRYcEMjwmzOkqHaOErpVQHJeWUkTgkCpHWZqR7Ly18pZTqgOOVteSU1jDdx4ZzQAtfKaU6JNk1/z7Rx96wBS18pZTqkJ05JwgJcjB+QG+ro3SYFr5SSnVA0pEyJg+KJDjQ9+rT9xIrpZRFqusaSM2v8Mnxe9DCV0opt+0+epJGp2HaEN8bvwctfKWUctv27FICHOJzJ1x9TQtfKaXctC2rlIkDexPew62V5b2OFr5SSrnh1OlG9uSeZNbwPlZH6TQtfKWUcsOuoyeobzTMGu6bwzmgha+UUm7ZluXb4/egha+UUm7x9fF70MJXSql2+cP4PWjhK6VUu/xh/B608JVSql3+MH4PWvhKKdUufxi/By18pZQ6I38ZvwctfKWUOiN/Gb8HLXyllDojfxm/By18pZQ6I38ZvwctfKWUalN1XYPfjN+DFr5SSrVpR04Z9Y2Gc0fGWB3FI7TwlVKqDZsPlRAc6PDJC5a3RgtfKaXasCmzhOlDowgJCrA6ike0W/giMkhEvhKRNBFJFZGHW9lmvoiUi8ge1+1XXRNXKaW6R3FlHQcLK5nrJ8M5AO687dwAPGqM2SUivYBkEVlrjDnQYruNxpgrPR9RKaW635bDJQDMGxlrcRLPafcI3xhTYIzZ5bpfCaQBA7s6mFJKWWnToRIiQ4MYNyDC6ige06ExfBEZCkwBtrfy8GwRSRGR1SIy3hPhlFLKCsYYNmeWMGdEHwIcYnUcj3G78EUkHFgBPGKMqWjx8C5giDFmEvAs8GEbz7FYRJJEJKm4uLiTkZVSqmtllVSTX17rV+P34Gbhi0gQTWX/ljHm/ZaPG2MqjDFVrvurgCAR+caeMsYsM8YkGmMSY2P9Z1xMKeVfNmc2jd/7y/z7r7kzS0eAV4A0Y8zTbWzTz7UdIjLD9bylngyqlFLdZdOhEuKjejI4OtTqKB7lziyducDtwD4R2eP63OPAYABjzFLgOuABEWkATgE3GWOM5+MqpVTXamh0sjWrlCsm9sd1HOs32i18Y8wm4IzftTHmOeA5T4VSSimr7D1WTmVtg9+N34OeaauUUv9hQ0YxImjhK6WUv1uXXsyk+Eiiw4KtjuJxWvhKKeVSVn2alLyTnD/KP2cRauErpZTLxkPFGAPzR2vhK6WUX1ufXkxUaBDnxEdaHaVLaOErpRTgdBo2HCpmXkKsXy2n0JwWvlJKAan5FZRUnfbb4RzQwldKKQDWZxwHYF6CFr5SSvm1denFTBzYm9hePayO0mW08JVStldeU8+uoyf8djrm17TwlVK2tymzBKcfT8f8mha+Usr21qUfJyIkkMmDIq2O0qW08JVStuZ0Gr5KP875o+MIDPDvSvTv704ppdqxJ+8kJVWnuWhsnNVRupwWvlLK1j4/UESAQ5g/SgtfKaX82hdpx5k+NIreoUFWR+lyWvhKKdvKLashvaiSi8b2tTpKt9DCV0rZ1udpRQBa+Eop5e8+TytiZFw4Q2PCrI7SLbTwlVK2VFFbz/asMi60weycr2nhK6VsaX16MQ1OY5vhHNDCV0rZ1BdpRUSFBjF1cJTVUbqNFr5SynbqG518lV7MBWPi/PZiJ63RwldK2c72rDLKT9Vzybh+VkfpVlr4SinbWb2/gJ5BAX6/HHJLWvhKKVtpdBrWpBZxwZhYegYHWB2nW2nhK6VsJfnICUqq6rh0Qn+ro3Q7LXyllK18ur+Q4AAHC8bYZ/7919otfBEZJCJfiUiaiKSKyMOtbCMi8oyIZIrIXhGZ2jVxlVKq84wxrEktZF5CDOE9Aq2O0+3cOcJvAB41xowFZgEPisi4FttcBiS4bouBJR5NqZRSHrA3r5xjJ09x6QR7zc75Wrv/xRljCoAC1/1KEUkDBgIHmm22CPiHMcYA20QkUkT6u75WWaCqroH9x8o5WFBBelElheW1FFfVcaK6HmMMBugR6CAqLJjo0GAGRYcyIi6cUXHhnBMfabs3s5Q9rN5fSKBDuHicfc6uba5Dr2lEZCgwBdje4qGBQG6zj/Ncn/uPwheRxTS9AmDw4MEdjKrak11Szap9BazPKGbXkRM0OA0AUaFBxEeFEhveg1FxvXA4BAFqG5ycrDlNfnktW7NKqTndCECgQ5gY35uZw/pw8bg4pgyKwmGjk1OUfzLG8On+AmaP6ENkaLDVcSzhduGLSDiwAnjEGFPR8uFWvsR84xPGLAOWASQmJn7jcdVxtfWNfLj7GMuT80g+cgKACQMj+O55w5kxLJrx/SOI7dUDkTMXtjGGgvJaDhZWsDPnBDuzy3hlUxZL1x8mtlcPLh3fj+umxXNOfO92n0spb5ReVElOaQ3fPW+41VEs41bhi0gQTWX/ljHm/VY2yQMGNfs4Hsg/+3iqLZW19by1/Sgvb8ympKqOhLhwfnbZGK6eMpC+ESEdfj4RYUBkTwZE9mTBmKaXuxW19Xx18DhrUgtZnpzLG9uOMK5/BDfPHMy3pw4kNNh+b3op3/VxSgEOwXZn1zbX7r9YaTqcewVIM8Y83cZmK4GHROQdYCZQruP3XaPRafhXUi7/syad0urTzEuI4YH5k5k9vI/Hj7wjQoJYNHkgiyYPpKK2no92H+OfO3L55Yf7+cvaDO6eM5Q7Zg+1xaXhlG8zxrAyJZ85I2KI7dXD6jiWcecQbS5wO7BPRPa4Pvc4MBjAGLMUWAVcDmQCNcDdHk+qSMk9yeMf7CM1v4IZQ6N55a6xTB4U2S1/d0RIELfPHspts4aQdOQES9Yd5qm1GSxdf5h7zh3G4vOG0ytEi195p5S8co6W1fDQBSOtjmIpd2bpbKL1Mfrm2xjgQU+FUv+pvtHJs19m8vxXmcSG9+DZm6dw5Tn9LRlLFxGmD41m+l3RpBVU8NxXmTz7ZSb/3H6Uhy9K4OYZgwkK0PP5lHdZuSef4AAHC206HfNr+i/Ty+WW1XDtC1t45otDLJo8gDU/OI+rJg3wijdOx/aP4PlbpvLRg3MZGRfOrz5KZeFfN7DlcInV0ZT6X41Ow8d785k/OpbePe39KlQL34ttOlTCVc9t4khpNUtuncrTN0z2yl/YSYMieWfxLF65M5GGRsMtL23nh+/uoaSqzupoSrE9u5TjlXV8a/IAq6NYTqdZeKmXN2bx5Ko0RsaFs+z2RK+/yLKIcOHYvswdGcPzX2WydP1hvjh4nF9eOY5vTx3oFa9IlD39OyWfsOAALhxjz5OtmtMjfC/jdBqeXJXGE5+kccm4fnzwX3O9vuybCwkK4NFLRrP64XmM7tuLx5ancN8byXq0ryxxusHJqn2FXDyur549jha+V6lvdPKj9/aybEMWd84ewgu3TiXMRxd4GhnXi3cWz+IXV4xlXUYxC/+ygc9SC62OpWxm46Fiyk/V63COixa+l2hodPLwO7tZsSuPH148it98a7zPL2fgcAj3zhvOvx86l74RISx+I5lffLiP2vpGq6Mpm/hoTz6RoUGcO9JeV7Zqixa+F2h0Gh5bnsKqfYX84oqxfP/CBL8a8x7drxcfPjiXxecN581tR7lu6RaOlFZbHUv5uYraej47UMgVE/sTHKhVB1r4lnM6DY+/v48P9+Tzo4WjuXeef67zERzo4PHLx/LyHYnklp3iymc28el+PRlbdZ1VewuorXdyfeKg9je2CS18iz21Np13k3L53oKRPGiDswAvGteXj793LsPjwrn/zV38fnUajU5dR0953nvJeYyMC2dSfG+ro3gNLXwL/WtnLs9/dZibZwzihxePsjpOtxkUHcry+2Zz68zBvLg+i3tf30lFbb3VsZQfyS6pJunICb49Nd6vhkfPlha+RTYdKuHxD/YxLyGG3y2aYLtfyuBAB//vmok8cfUENh4q4ernN5NVXGV1LOUn3t+Vh0PgmikDrY7iVbTwLXCktJoH3kpmZFw4L9w61dZrz9w2awhv3TuTkzX1LHp+M+sziq2OpHyc02lYkZzHvIRY+vXu+FLh/sy+TWORU6cbue+NZAIcwkt3JOoKk8DM4X1Y+dBc4qNCuee1nbyz46jVkZQP25pVSn55LddNi7c6itfRwu9Gxhh+/sE+0osq+euNkxkUHWp1JK8RHxXK8vtnM3dkDD99fx9PfZZO0yKsSnXMe8l59AoJtO11a89EC78bvb0jl/d3H+ORC0cxf3Sc1XG8TniPQF65M5EbEwfx7JeZPPqvFE43OK2OpXxIRW09n+4v5KpJAwgJ0qUUWvLN8/Z9UObxKn73cSrzEmL43gL/n37ZWUEBDv7w7YnER/XkqbUZFFbUsuS2aV65SqjyPh/tPsap+kZu1Ln3rdIj/G5wusHJI+/uJjQ4kKeun+TzSyZ0NRHhexcm8PQNk9iRXcaNL27leGWt1bGUlzPG8Nb2o4wfEME5Ove+VVr43eDptRnsP1bBH66dSFwnLjBuV9dOjefvd0/naFkN1y/dSm5ZjdWRlBfbnXuSg4WV3DJzsO2mObtLC7+LJeWU8eKGppOrLhlv78urdca8hFjedE3bvG7pFjKKKq2OpLzUP7cfJSw4gEWTde59W7Twu1BtfSM/WbGXAb178osrxlkdx2dNHRzFu/fNwmnghhe3sif3pNWRlJcpP1XPx3vzWTRlIOE+uqR4d9DC70LPf5XJ4eJqfn/tRJ9d195bjOkXwYr759ArJJBbX9rGlky9bq76Px/syqO23sktMwZbHcWraeF3kbSCCpasO8y1Uwdy3ihdi9sTBvcJ5b375xAfFcpdf9/J2gNFVkdSXsAYwz93HGVSfG8mDNQ3a89EC78LNDoNP12xl949g/ilDuV4VN+IEN69bxZjB0TwwJvJrN6nSyzb3basMjKKqrh15hCro3g9Lfwu8M/tR0jJK+dXV40jKizY6jh+JzI0mDe+M4NJgyJ56O3d/Dsl3+pIykJ/35xNVGiQXsbQDVr4Hnai+jT/81kGc0b04VuT9Bewq0SEBPH6PTOYNiSKh9/Zzfu78qyOpCyQW1bD2rQibpk5WM+sdYMWvoc9vTaDqroGfn3VeJ0L3MXCewTy2t3TmTW8D48uT+FfO3OtjqS62etbcggQ4fZZQ62O4hO08D3oQH4Fb20/wu2zhjC6Xy+r49hCaHAgr941nXNHxvDjFXt5a/sRqyOpblJd18C7SblcNrG/LoPspnYLX0ReFZHjIrK/jcfni0i5iOxx3X7l+ZjezxjDb/6dSu+eQfzgIvtcvcobhAQF8NIdiSwYE8fPP9jPa5uzrY6kusGKXXlU1jZw99yhVkfxGe4c4b8GXNrONhuNMZNdt9+dfSzfsya1kB3ZZTy2cDS9Q3Whr+4WEhTA0tumccm4vvzm3wd4eWOW1ZFUF3I6Da9tyWHSoEimDo6yOo7PaLfwjTEbgLJuyOKzGhqd/GlNOglx4dw0XU/8sEpwoIPnb53K5RP78cQnaSzbcNjqSKqLfHHwOFnF1dyjR/cd4qkx/NkikiIiq0VkfFsbichiEUkSkaTiYv+5lN3y5Dyyiqv50cLRBOhKmJYKCnDwt5umcMU5/Xly1UGWrNPS9zfGGF5Yl0l8VE+umNjf6jg+xRPn++8ChhhjqkTkcuBDIKG1DY0xy4BlAImJiX5xOaNTpxv56+cZTBsSpVfY8RJBAQ7+duNkAkT446cHcRrDgxfoNQj8xY7sMnYfPcl/LxpPoI2vB90ZZ134xpiKZvdXicgLIhJjjLHFYievbcmhqKKOZ2+eqtMwvUhggIOnb5iEQ+DPa9JxOg3fu7DV4xDlY5asP0yfsGCu14ucdNhZF76I9AOKjDFGRGbQNExUetbJfEB5TT1L1mWyYEwcM4ZFWx1HtRAY4OCpGybjEOGptRk0GsMjOoPKpx3Ir2BdejGPXTJKT7TqhHYLX0TeBuYDMSKSB/waCAIwxiwFrgMeEJEG4BRwk7HJ1adf2ZxNRW0DP1o42uooqg0BDuHPrquM/fXzQzgN/OCiBH015qOWrj9MWHCAnmjVSe0WvjHm5nYefw54zmOJfET5qXr+vjmbyyb0Y2z/CKvjqDMIcAh/+vY5OASe+eIQxhh+ePEoLX0fk1NSzcd78/nOucN06nMn6SLtnfT6lhwqaxt4SC9I7hMcDuEP156DQ4Rnv8yk0Wn40cLRWvo+5JkvDxEc6OC75w23OorP0sLvhMrael7ZlM3F4/oyfoCuv+0rHA7hyWsm4nAIL6w7jNPATy7V0vcFh4ur+HD3Me6dN5y4XrqMQmdp4XfCP7YeofxUPd9foLM+fI3DITyxaAIOaRoPdhrDzy4bo6Xv5Z754hAhQQHcp0f3Z0ULv4Oq6xp4eWMWF4yOZWK8Ht37IodD+O9FE3CIsGxDFo1Owy+uGKul76UOFVWyMiWf+88fQZ/wHlbH8Wla+B309o6jnKip1zndPk5E+O23xuMQ4ZVN2TiN4VdXjtPS90J//eIQoUEBLJ6nR/dnSwu/A+obnfx9cw4zhkXrgk1+QET49VXjcIjw6uZsjIFfX6Wl7032Hyvnk70FPHTBSL16nAdo4XfAqn0FHDt5it9+q83lgpSPERF+eeVYAhzw0sZsGp2m6chf10SynDGGJ1elERUaxOLz9ejeE7Tw3WSMYdmGLEbEhrFgTJzVcZQHiQiPXz4Wh0N4cX0WTmOaxvi19C21Lr2YLYdL+c1V44gI0Xn3nqCF76ath0tJza/gD9dO1CLwQyLCTy8dg0OEJeuaZu/8v6v1Z22VhkYnT65KY1hMGLfMHGJ1HL+hhe+mFzdkERPeg6unDLQ6iuoiIsKPF44mQITnvsrE6YTf63/wllienMeh41UsvW0qwYG6IqanaOG7Ib2wkvUZumCTHYgIj14yCodDeOaLQzQawx+/fY5e56AbVdU18PTaDBKHRLFwfD+r4/gVLXw3vL41hx6BDm7Vl5a2ICL88OJROATXgmuGP183SUu/m/zt8wxKqup46Y5EnTHlYVr47Sg/Vc8Hu46xaPIAnRZmM49cNAqHCE+vzcAY+J/rtfS7WnphJa9uzuGm6YOYPCjS6jh+Rwu/HSuS8zhV38gds4daHUVZ4PsXJjQtsbwmHacxPHX9JL3KUhcxxvDLD/cTERLIjxeOsTqOX9LCPwOn0/DGtiNMHRzJhIG6jIJdPXjBSETgT5+m4zTwlxu09LvC+7uOsSOnjD9cO1FfTXcRLfwz2JRZQnZJNQ/fONnqKMpi/zV/JAEi/H71QZxOw19vmkyQlr7HlNfU8/vVaUwZHMkNeunCLqOFfwb/2JpDTHgwl03UmQIK7jt/BAEO4YlP0nAawzM3T9HS95DffpzKiZp6Xr9HT3jrSvrb2obcshq+OHicm6YPpkegTsVUTe6dN5xfXjmO1fsLeeifuzjd4LQ6ks/7Iq2I93cd48H5I/T6El1MC78Nb20/ikOEW2YOtjqK8jLfOXcYv75qHGtSi1j8RhI1pxusjuSzymvqefyDfYzp14uH9PoSXU4LvxX1jU7eS85lwZg4BkT2tDqO8kJ3zx3G76+dyIaMYm55aTtl1aetjuSTfvfxAUqqTvPn6ybpGbXdQPdwK75IO05J1Wlumq5vHqm23TxjMEtum0ZaQQXXLd1CblmN1ZF8ysqUfFbsyuPB+SP0YkLdRAu/Ff9KyiWuVw/OHxVrdRTl5RaO78eb986kpLKOby/ZQlpBhdWRfMLR0hp+/v4+pg6O5Pt6MaFuo4XfQmF5LevSj3N9YrzOtVZumT40mvcemEOAQ7hh6Va2Hi61OpJXq2908v13doPA326aov/OupHu6RZW7MrDadC5wKpDRvXtxYoH5tCvdwh3vrqDFcl5VkfyWn/69CB7ck/yh2vPYVB0qNVxbEULvxmn0/DuzlxmDY9mSJ8wq+MoHzMgsifv3T+H6cOieHR5Cn/6tOkkLfV/PtpzjJc2ZnPH7CFccU5/q+PYjhZ+M9uySzlaVsON+mat6qTeoUG8dvcMbp4xmBfWHea/3tql0zZd9h8r5ycr9jJjWDS/vHKc1XFsSQu/mX/tzKVXSCCXTdAjD9V5QQEOnrxmAr+4YixrDhRy44vbKCyvtTqWpYor67jvjWSiQoN54dapeoayRdrd6yLyqogcF5H9bTwuIvKMiGSKyF4Rmer5mF2vorae1fsLWTR5gF7kRJ01EeHeecN5+Y5EsoqruPLZTWzPsuebudV1Ddzz2k7Kqk/z4u3TiAnvYXUk23Lnv9nXgEvP8PhlQILrthhYcvaxut+n+wqpa3By3TQdzlGec+HYvnzw4FwiQgK55eXtvLwxC2PsM65f3+jkgbd2caCggudvncI58ZFWR7K1dgvfGLMBKDvDJouAf5gm24BIEfG5MZEPdh9jWEwYk/QEEOVho/r24qOH5nLR2Die+CSNh/65m6o6/x/XdzoNP3lvLxsyinnymgksGNPX6ki254mBtIFAbrOP81yf+wYRWSwiSSKSVFxc7IG/2jPyT55iW3YpV08eqJdUU12iV0gQS2+bxs8uG8Pq/QUsem4TqfnlVsfqMk6n4Scr9vL+7mM8evEobpyua1J5A08UfmsN2eprVmPMMmNMojEmMTbWe85iXZmSjzGwaPIAq6MoPyYi3Hf+CN68dyaVtQ1c8/wWXt6Y5XdTN78u++XJeTx8YQLf0zNpvYYnCj8PaD7wHQ/ke+B5u82Hu48xZXAkQ2N07r3qenNGxPDpI+dx/uhYnvgkjTte3UFRhX/M4qlraOSRd/f8b9n/4OJRVkdSzXii8FcCd7hm68wCyo0xBR543m6RVlDBwcJKrpnS6iiUUl0iOiyYZbdP48lrJpJ0pIxL/rKB95LzfPoN3fJT9dz56g5WpuTz40tHa9l7oXaveCUibwPzgRgRyQN+DQQBGGOWAquAy4FMoAa4u6vCdoUP9xwj0CFcMdHn3mdWPk5c11uYOTyaH7+3l8eWp/DRnmM8ec1En1ty4EhpNd/9RxLZJdX89cbJXK0HUF5JrDqiSExMNElJSZb83V9zOg1z/vAl4wdE8Mpd0y3NouzN6TS8uf0If1x9EKeBRy5K4K65Q33iamtrUgt5bHkKAiy9bRpzRsZYHcmviUiyMSaxM19r69PdtmWXUlhRq0cjynIOh3DH7KGs/eH5zB3Zh9+vPsglf9nAmtRCrx3mqa1v5ImPD3DfG8kMiwnjk+/P07L3crYu/I/3FhAaHMBFY3V+sPIOAyJ78vKd03n9nhkEBzi4741kbn5pG7uOnrA62n9Iyinj8mc28vKmbG6bNZjl98/2uWEoO2p3DN9fNTQ6+XR/IReO7UvPYO9/2azs5fxRscx9eB5v78zlL2szuPaFLcxLiOF7CxKYMSzaslzHK2v5y9pDvLPzKAN69+SN78xgXoL3TLFWZ2bbwt+eXUZZ9WmumNjP6ihKtSowwMHts4Zw7ZSBvLX9CMs2ZHHDi1tJHBLFHXOGcun4ft12Hdjymnpe3ZzNSxuzON3g5K45Q3nsktGE9bBthfgk2/60vh7OmT86zuooSp1RWI9AFp83gttnDeXtHUd5bUsO3397NzHhPbhuWjxXTerPuP4RXXKWeFZxFa9vyWF5ch41pxu5YmJ/frRwtJ6z4qNsWfgNjU7WpDYN5+jKmMpX9AwO4J5zh3HXnKGsP1TMm1uP8NLGLJauP8zw2DAWju/H3BExJA6N6vTvtTGGI6U1fJ5WxL9T8knJKycoQPjWpIF859xhjBsQ4eHvSnUnWxb+tqyvh3N07r3yPQ6HcMHoOC4YHUdZ9WlW7y/g45QCXtqQxZJ1hwkOdHDOwN6M7R/BmP69GBQVSr/eIcSE9yA40EGgQ2h0GspP1VN+qp4jpTVklVRxsKCSHdllFLrO+p0wMIKfXTaGa6YMJC4ixOLvWnmCLQv/k335hAUHMH+0vtmkfFt0WDC3zhzCrTOHUFXXwI7sUjZnlrI37yQf7j5G5Tb3V+XsFxFC4tAoZg3vw7kjY3TYxg/ZrvCbz87R4RzlT8J7BLJgTN//XYbYGMOxk6fIP1lLUUUtJVV11Dc6qW80BDiE3j2D6N0ziPiongyLCaNXSJDF34HqarYr/K1ZpZyoqdcLKCu/JyLER4USH6Xz41UT2514tWpfAWHBAZw/SodzlFL2YqvCb3Qa1h4o4oIxcTqco5SyHVsV/u6jJyipOs3C8XqylVLKfmxV+J8dKCIoQHR2jlLKlmxT+MYY1qQWMmdEjM5GUErZkm0K/9DxKo6U1nDJeF0ZUyllT7Yp/M9SCwG4WJdCVkrZlG0Kf+2BIqYMjtRTxJVStmWLwi8oP0VKXjmXjNPZOUop+7JF4X9+oAiAi8fpcI5Syr5sUfifHShieGwYI+PCrY6ilFKW8fvCr6itZ+vhUj26V0rZnt8X/saMEhqcRmfnKKVsz+8L/8uDx+ndM4jJgyKtjqKUUpby68J3Og3rM45z/qhYAgP8+ltVSql2+XUL7jtWTknVaRaM0QuVK6WUXxf+lwePI4Kufa+UUrhZ+CJyqYiki0imiPy0lcfni0i5iOxx3X7l+agd91X6caYMiiQqLNjqKEopZbl2L3EoIgHA88DFQB6wU0RWGmMOtNh0ozHmyi7I2CnFlXXszSvnsUtGWR1FKaW8gjtH+DOATGNMljHmNPAOsKhrY529denHAbhAx++VUgpwr/AHArnNPs5zfa6l2SKSIiKrRWS8R9Kdha/Sj9M3ogfj+kdYHUUppbyCO4UvrXzOtPh4FzDEGDMJeBb4sNUnElksIkkiklRcXNyhoB1R3+hkY0YJF4yOQ6S1+EopZT/uFH4eMKjZx/FAfvMNjDEVxpgq1/1VQJCIxLR8ImPMMmNMojEmMTa262bOJOWcoLKuQYdzlFKqGXcKfyeQICLDRCQYuAlY2XwDEeknrkNpEZnhet5ST4d11/qMYoIChLkjv/F/jlJK2Va7s3SMMQ0i8hCwBggAXjXGpIrI/a7HlwLXAQ+ISANwCrjJGNNy2KfbbDxUzNTBUYT3aPfbU0op23CrEV3DNKtafG5ps/vPAc95NlrnlFTVkZpfwY8WjrY6ilJKeRW/O9N2c2YJAPMSdDhHKaWa87vC35BRQlRoEOMH9LY6ilJKeRW/KnxjDBsPFTN3ZAwBDp2OqZRSzflV4WcUVXG8so7zEnSxNKWUasmvCn/joaaTuc7V8XullPoGvyr8DYdKGBkXzoDInlZHUUopr+M3hV9b38j2rFKdnaOUUm3wm8JPyjlBXYNTx++VUqoNflP4GzObllOYOTza6ihKKeWV/KfwM0pIHBJNaLAup6CUUq3xi8IvqarjQEGFzs5RSqkz8IvC33q4aWFOXR1TKaXa5h+Fn1VKrx6BTBigV7dSSqm2+EfhHy5lxrBoAgP84ttRSqku4fMNWVheS3ZJNbNH9LE6ilJKeTWfL/ytWU3LIWvhK6XUmfl84W/JLCUyNIix/XT8XimlzsTnC39rVikzh0Xj0OWQlVLqjHy68HPLasg7cYo5I3Q6plJKtcenC//r+fc6fq+UUu3z7cLPKiUmPJiEuHCroyillNfz2cI3xrD1cCmzhvdBRMfvlVKqPT5b+Nkl1RRW1OpwjlJKuclnC39rVtP4vb5hq5RS7vHdwj9cSr+IEIb2CbU6ilJK+QSfLHxjDNuzy5g1PFrH75VSyk0+Wfg5pTUUV9YxY5iO3yullLt8svB3ZpcBMGNYlMVJlFLKd7hV+CJyqYiki0imiPy0lcdFRJ5xPb5XRKZ6Pur/2ZFTRnRYMCNidf69Ukq5q93CF5EA4HngMmAccLOIjGux2WVAguu2GFji4Zz/YUd2GYlDonT8XimlOsCdI/wZQKYxJssYcxp4B1jUYptFwD9Mk21ApIj093BWAIoqajlaVsOMYdFd8fRKKeW3At3YZiCQ2+zjPGCmG9sMBAqabyQii2l6BQBQJyL7O5S2me/+Eb7b2S/2jBigxNoIZ0XzW8uX8/tydvD9/KM7+4XuFH5r4yamE9tgjFkGLAMQkSRjTKIbf79X0vzW0vzW8eXs4B/5O/u17gzp5AGDmn0cD+R3YhullFIWcqfwdwIJIjJMRIKBm4CVLbZZCdzhmq0zCyg3xhS0fCKllFLWaXdIxxjTICIPAWuAAOBVY0yqiNzvenwpsAq4HMgEaoC73fi7l3U6tXfQ/NbS/Nbx5exg4/xizDeG2pVSSvkhnzzTVimlVMdp4SullE10W+GLyPUikioiThFpc0pUe8s4WEVEokVkrYgccv3Z6kI+IpIjIvtEZM/ZTJ/yBG9bEqOj3Mg/X0TKXft6j4j8yoqcbRGRV0XkeFvnm3jz/ncju7fv+0Ei8pWIpLl65+FWtvHm/e9O/o7/DIwx3XIDxtJ0wsA6ILGNbQKAw8BwIBhIAcZ1V8Z28v8J+Knr/k+BP7axXQ4Q4wV5292XNL3Rvpqm8yhmAdutzt3B/POBj63Oeobv4TxgKrC/jce9ef+3l93b931/YKrrfi8gw8d+/93J3+GfQbcd4Rtj0owx6e1s5s4yDlZZBLzuuv86cLV1UdziVUtidII3/y64xRizASg7wyZeu//dyO7VjDEFxphdrvuVQBpNZ/8358373538HeZtY/htLdHgDfoa17kFrj/j2tjOAJ+JSLJrKQmruLMvvXl/u5tttoikiMhqERnfPdE8xpv3vzt8Yt+LyFBgCrC9xUM+sf/PkB86+DNwZ2mFjgT7HOjXykM/N8Z85M5TtPK5bps3eqb8HXiaucaYfBGJA9aKyEHX0VJ389iSGBZxJ9suYIgxpkpELgc+pGnFVl/hzfu/PT6x70UkHFgBPGKMqWj5cCtf4lX7v538Hf4ZeLTwjTEXneVTWLpEw5nyi0iRiPQ3xhS4XvYdb+M58l1/HheRD2gamrCi8H19SYx2szX/B2CMWSUiL4hIjDHGVxbG8ub9f0a+sO9FJIimsnzLGPN+K5t49f5vL39nfgbeNqTjzjIOVlkJ3Om6fyfwjVcsIhImIr2+vg9cAnR6RdCz5OtLYrSbX0T6iTRdFEFEZtD0+1za7Uk7z5v3/xl5+753ZXsFSDPGPN3GZl67/93J36mfQTe+63wNTf+j1gFFwBrX5wcAq5ptdzlN70gfpmkoyPJ3zF25+gBfAIdcf0a3zE/TjJIU1y3V6vyt7UvgfuB+132h6eI2h4F9tDF7yovzP+TazynANmCO1Zlb5H+bpiXC612/+9/xlf3vRnZv3/fn0jQ8sxfY47pd7kP73538Hf4Z6NIKSillE942pKOUUqqLaOErpZRNaOErpZRNaOErpZRNaOErpZRNaOErpZRNaOErpZRN/H96r+4SpQq3WgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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 gonna be near 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": "iVBORw0KGgoAAAANSUhEUgAAAlMAAADSCAYAAABuMkW8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6r0lEQVR4nO3dd3hUZfbA8e8hgFQB6SJNRRRQUJBiAXSxYUERFUVRV+WHZa1rX3XtvSEqogioEVdpFnARFxBURAGRIipFKYJU6aGEnN8f50ZCCMmETOZOOZ/nmSeTuXfmnrmZeXPuW0VVcc4555xz+6ZE2AE455xzziUyT6acc84554rAkynnnHPOuSLwZMo555xzrgg8mXLOOeecKwJPppxzzjnnisCTKbcbEflNRDoF9+8RkTdCiqOjiCwN49jOuYJ5WRE9IjJERM4NO459UdD5F5FNInJwBK+zn4j8JCI1ohthbHgylUBEpLuITBGRzSKyMrh/nYhIcRxPVR9T1auL+joi0kBEVERKRiOusInIIBF5JOw4nNsbLyviQyRlhYgcBTQHPoxNVLGlqhVUdWEE+20D3gTuLP6oos+TqQQhIrcBLwJPA7WAmkBv4Hig9F6ekxazAJ1zccHLioTzf0C6+gzaAO8Cl4vIfmEHUmiq6rc4vwGVgM3A+QXsNwh4FRgd7N8JOBP4HtgALAH+nes5lwGLgDXAvcBvQKdg27+Bd3Ls2xb4GlgH/AB0zLFtAvAw8BWwEfgMqBZsWwwosCm4tcsj9rJB/H8CPwK3A0tzbD8QGAasAn4FbsyxrTUwNXiPK4Dncmw7IUfMS4Argsf3A54JYlsB9APKBts6AkuB24CVwHLgymBbL2AHsD14Lx+H/fnwm9+yb15WJF5ZASwETsjx+w853v+m4Hx0DLadA8wJYpwAHJHjeUcEj60L9jkn19/7FeDT4DW/whLtF4Lz+BNwdITnMN/zn8f7U+DQHHG8DIwK/vZTgENy7T8P6BD2d6nQ372wA/BbBH8kOB3IBEoWsN8gYD12BVoCKBN82Y8Mfj8qKAzODfZvEnyx2gcFxnPBcfYoIIE6WCHaOXitU4LfqwfbJwALgMOCL9sE4IlgW4PgC7XX+IEngEnAAUBdYHb2FzQ43jTgfuzK+mCsADot2D4ZuCy4XwFoG9yvF3xhLwZKAVWBFsG2F4CPguNVBD4GHg+2dQzOw0PB8zoDW4AqOc7zI2F/Lvzmt9w3LysSq6wAygfvt/petvfCEp39g/O1OTifpYA7gPnB+ywV3L8n+P3k4P00zhHHaqBl8LcehyVJPYE04BFgfITncK/nfy/vIXcytRZLaksC6cB7ufb/iBzJW6LcvJkvMVQDVqtqZvYDIvK1iKwTkQwRaZ9j3w9V9StVzVLVrao6QVVnBb/PBIYAHYJ9uwGfqOpEtfbq+4CsvcRwKTBaVUcHrzUWu8LrnGOfgar6i6pmAO8DLQrxHi8EHlXVtaq6BOiTY9uxWGHzkKpuV2t/fx3oHmzfARwqItVUdZOqfhM83gP4XFWHqOoOVV2jqjOCfiPXALcEx9sIPJbj9bJf86HgeaOxfySNC/F+nAuDlxWJVVZUDn5uzL1BRE7AkpxzVHUDcBEwSlXHquoOrLasLHAcVhNYAUtKt6vqOOATLDnMNkJVp6nqVmAEsFVV31LVncB/gKMjPIf5nf9IDFfVb4PPaDp7/u035jgvCSMpOvmlgDVANREpmV1IqupxAMEoipxJ8ZKcTxSRNtiVRDPsKmM/4INg84E591fVzSKyZi8x1AcuEJGzczxWChif4/c/ctzfgn25I7VbLFhzQs5jHygi63I8loZdHQFchV0Z/iQivwIPquon2FXTgjyOVR0oB0zL0R9XgtfMtibnP6R9eD/OhcHLisQqK7LjrAhs/esAInWxJPNyVf0lePhAcrxXVc0SkSVYTWAmsERVcya4i4Jt2VbkuJ+Rx+/ZMRd0DvM7/5Eo6G9fkV3nJWF4MpUYJgPbgC5YO3Z+cndifBfoC5yhqltF5AXs6hWsff+I7B1FpBxWvZ2XJcDbqnpN4ULPM6a8LMcKtDnB7/VyHftXVW2U54urzgMuFpESQFdgqIhUDZ7XOo+nrMYKj6aq+ntkb2H3Q+7Dc5yLBS8rEqisCJLS7CbPVQAiUhYYCbygqp/m2H0Z1gxLsJ9g5+F3YCdQV0RK5Eio6gG/UHj5nkPyP//RcATwbJRfs9h5M18CUNV1wIPAKyLSTUQqiEgJEWmBtbnnpyKwNigcWwOX5Ng2FDhLRE4QkdLYFdvePhPvAGeLyGkikiYiZYL5RQ6K4C2swpoE8ptr5H3gbhGpErzmP3Js+xbYICJ3ikjZ4PjNRORYABG5VESqB4XIuuA5O7Eq5E4icqGIlBSRqiLSItjvdeD57DlNRKSOiJwWwXsBu6IrcN4U52LNy4qELCtGs6s5FWx6gJ9U9ak83veZIvI3ESmFdXrfhnWan4L1p7pDREqJSEfgbOC9COPMKd9zSP7nv0hEpA7WF+ubgvaNN55MJYjgi3Ur1ulwJfYlfQ2bk+PrfJ56HfCQiGzEOhS+n+M15wDXY1eky7HRGXlOvha0jXfBOjiuwq5ebieCz5CqbgEeBb4K+m60zWO3B7Hq4l+x0T1v53j+TqxgaBFsXw28gY1cAut0O0dENmFDwrsHfUAWY/00bsM6Pc7A5nMBO2/zgW9EZAPwOZH3cxgANAney8gIn+NcTHhZkXBlRX+gh+xqR+wOnCc22WX27URV/Rnrj/ZS8L7OBs4O+jVtx0b6nRFsewXoqao/RRjnXyI4h3s9/1FwCTA46JeXUETVWyycc865sIjIu8D7qjoy7FjCIja31A9Ae1VdGXY8heXJlHPOOedcERRY7Rq0d38rIj+IyBwReTCPfURE+ojIfBGZKSLHFE+4zjkXOS+/nHOxEMlovm3Ayaq6Kej09qWIfJpjfg6wdtpGwa0NNrNum6hH65xzhePll3Ou2EXSIVBVdVPwa6nglrttsAvwVrDvN0BlEakd3VCdc65wvPxyzsVCRKP5gqGRM7CRIWNVdUquXeqw+yReS9l9sjDnnAuFl1/OueIW0aSdwVDJFiJSGRghIs1UdXaOXSSvp+V+QER6YWsNUb58+ZaHH3544SN2zhWrP/6A33+H5s2hZJSn9Z02bdpqVa0e3VfNX7TKL/AyzLlUll/5VaiiUlXXicgEbK6OnIXRUmxG1GwHYbO15n5+f2xODVq1aqVTp04tzOGdczFw5pmw//4wY0b0X1tECrv0RNQUtfwKXsPLMOdSVH7lVySj+aoHV3TZ09x3wlaxzukjoGcwKqYtsF5Vl+97yM65MOzcCV99BSeeGHYk0eHll3MuFiKpmaoNDBaRNCz5el9VPxGR3gCq2g+bDr8zNkvsFuDKYorXOVeM5syB9evhhBPCjiRqvPxyzhW7ApMpVZ0JHJ3H4/1y3FdsqQHnXAKbFKwLnyw1U15+Oediwdfmc8795csvoU4dqF8/7Eiccy5xeDLlnANA1WqmTjwRJK/xbc455/LkyZRzDoBFi2xKhCTqL+WcczHhyZRzDki+/lLOORcrnkw55wDrL1WpEjRtGnYkzjmXWDyZcs4BMHEiHH88pKWFHYlzziUWT6accyxfDj/9BB07hh2Jc84lHk+mnHNMmGA/Tzop1DCccy4heTLlnGP8eOsvdfQe01s655wriCdTzjkmTID27b2/lHPO7QtPppxLcb//DvPmeX8p55zbV55MOZfixo+3n95fyjnn9o0nU86luPHjoUoVaN487Eiccy4xeTLlXIobPx46dIASXho459w+8eLTuRS2aBH8+qs38TnnXFF4MuVcCvP5pZxzrug8mXIuhY0fD1Wr+np8zjlXFJ5MOZeiVGHcOJsSwftLOefcviuwCBWRuiIyXkTmisgcEbkpj306ish6EZkR3O4vnnCdc9Hy88+wZAmcemrYkRQfL7+cc7FQMoJ9MoHbVHW6iFQEponIWFX9Mdd+k1T1rOiH6JwrDmPH2s9TTgk3jmLm5ZdzrtgVWDOlqstVdXpwfyMwF6hT3IE554rXZ5/BoYdCw4ZhR1J8vPxyzsVCoXpKiEgD4GhgSh6b24nIDyLyqYh4d1bn4tj27db5PJmb+HLz8ss5V1wiaeYDQEQqAMOAm1V1Q67N04H6qrpJRDoDI4FGebxGL6AXQL169fY1ZudcEU2eDJs3p04yFY3yK3gdL8Occ3uIqGZKREphBVG6qg7PvV1VN6jqpuD+aKCUiFTLY7/+qtpKVVtVr169iKE75/bV2LGQlpYaixtHq/wKtnsZ5pzbQySj+QQYAMxV1ef2sk+tYD9EpHXwumuiGahzLno++wzatoVKlcKOpHh5+eWci4aPPsp/eyTNfMcDlwGzRGRG8Ng9QD0AVe0HdAOuFZFMIAPorqq6byE754rTmjUwdSr8+99hRxITXn4554pk4EC4+ur89ykwmVLVLwEpYJ++QN/CBOecC8e4cTZhZyr0l/LyyzlXFE8/DXfcYVPIZE8nkxef99i5FPPZZ9a816pV2JE451x8UrUk6o474KKL4JNP8t8/4tF8zrnEpwpjxsDJJ0NJ//Y759weMjOhVy9r3rvuOujTxwbs5MdrppxLIbNn2xIynTuHHYlzzsWfjAzo1s0SqQcegL59C06kwGumnEspo0bZT0+mnHNud6tXQ5cuNg/fSy/BDTdE/lxPppxLIaNGwTHHwIEHhh2Jc87Fj4UL4fTTYfFieP99q50qDG/mcy5FrFkDX38NZ54ZdiTOORc/vv3W5t1bswY+/7zwiRR4MuVcyhgzBrKyPJlyzrlsH39sK0FUqGAXmyecsG+v48mUcyli1CioXh2OPTbsSJxzLnyvvgrnngtNm1o/qcaN9/21PJlyLgVkZsJ//2sdz0v4t945l8KysuCuu2zag86dYcIEqFmzaK/pHdCdSwHffANr13oTn3MutW3ZApdfDkOHQu/eNmovGnPueTLlXAoYNcoKjFRYQsY55/Ly++829cH06fDMM3DrrSD5LjYVOU+mnEsBo0ZZx8pKlcKOxDnnYm/aNDjnHNiwAT76CM46K7qv770nnEtyCxfCrFlw9tlhR+Kcc7E3dCiceKLVzn/1VfQTKfBkyrmkN2KE/TzvvHDjcM65WFKFRx+FCy6AFi1sPqmjjiqeY3kzn3NJbsQIK0gaNgw7Eueci42MDLjmGkhPhx494I03oEyZ4jue10w5l8T++MMmovNaKedcqliyxJr10tPhkUfg7beLN5ECr5lyLql9+KFVdXsy5ZxLBV98Yc1627ZZR/NY9RX1minnktiIEXDIIdCsWdiROOdc8VGFPn3gb3+DqlWtf1QsB914MuVcklq/HsaNs1qpaM2l4pxz8SYjA664Am66yUbqTZlStKVh9kWByZSI1BWR8SIyV0TmiMhNeewjItJHROaLyEwROaZ4wnXORSQ9nVGH3siOHXDe212t80AK8vLLuSSUng4NGkCJEiw5qB0nNl3DW2/Bgw/C8OGw//6xDymSPlOZwG2qOl1EKgLTRGSsqv6YY58zgEbBrQ3wavDTORdr6enQqxfDtwymFstpu2Ik9Bpj23r0CDW0EHj55VwyCco3tmxhLJ3o8Xs62yjJx7dN4Kz7O4YWVoHJlKouB5YH9zeKyFygDpCzMOoCvKWqCnwjIpVFpHbwXFeAZctgxgyYM8dGX61aZQvTpqVBhQpw4IFQt64Nb2/SBEqXDjtiF9fuvZeNW0owijO5igGUQG1BqnvvTblkyssv55LMvfeyc8tWHuF+HuQBmvAjwzifxkO3wTO/hRZWoUbziUgD4GhgSq5NdYAlOX5fGjy2W2EkIr2AXgD16tUrZKjJQ9VGHIwYAaNHw/z5u7aVLQvVq1vClJVlU9+vXr1re+nS0Lq1rXTdubNNQOb9YdxuFi/mY7qzlbJ0573dHk9lRS2/gtfwMsy5sKxZw6pFm+nBp4zlVHoymFe4jvJsgcXh/iOMuAO6iFQAhgE3q+qG3JvzeIru8YBqf1VtpaqtqlevXrhIk8DmzfD883D44XDSSdC/Pxx2mD32xRewdq1VICxaBPPmwYIFVku1dSv8/DO89x7ceKPtc889VlN11FHw7LOwYkXY787FjYMO4j26U4elHMfXux5P4X/+0Si/wMsw50KxahXcdRdf1e3O0XzPRNrzOlcziCsskYLQy7eIkikRKYUVROmqOjyPXZYCdXP8fhCwrOjhJYdt2+CFF+Dgg22V6urVYfBgWLPGFqC9+WZo3x6qVMn7+fvtZ0nXRRfB00/bgo3Ll0O/flCxIvzzn9YMeNVV8NNPsXxnLh79eVQH/svpXMR/rIkPoFw5W1chBXn55VyCWrECbr8drd+AZ5/MpMPW/1Kmanm+KXMSVzNg11VQHJRvkYzmE2AAMFdVn9vLbh8BPYNRMW2B9d7fwHz5pdUg3XILHHmkzUb95ZfQs6f9/fdVrVrwf/9nrzd3rvXHe/dd61N1wQVWs+VS0OrVjPy8AjsozUW1JlobcP36Vg2aYv2lwMsv5xLS8uVW89CwIWufHUjXql/wT56hy3lpTFtQhRZv3GDlWjyVb6qa7w04AavyngnMCG6dgd5A72AfAV4GFgCzgFYFvW7Lli01mW3dqnrjjaqg2qCB6ujRxX/MlStV//Uv1fLlVUuWVL3hBnvMpZDbbtPT+K82PGibZmWFHcyegKlaQNkQzVtxlV+aAmWYczG3dKn94yxTRjUtTSec+qgeVGu7liql+vzzGnqZll/5FbNCLfctmQuihQtVW7Wys3vjjaqbNsX2+MuXq/burZqWplq5suobb4T/IXQx8PvvunK/gzRNMvXuu8MOJm+xTqaK85bMZZhzMbV4ser116vut59qWpruuPwqve/6NVqihGqjRqpTp4YdoMmv/PIZ0KPsq6+gVStrZhsxAl58EcqXj20MtWrBq6/CrFnWQf3qq63D+88/xzYOF2OPPMKwHeewU9O46KKwg3HOuQIsWgTXXmtrXr32Glx2GYvGL6TDvDd4+OUD6NkTpk+Hli3DDrRgnkxF0bBhti5QtWrWSfzcc8ON54gjYPx4eP11+OEHaN7ckjvNc5ySS2gLF8Lrr/NW9dto0sSSaOeci0u//grXXAOHHgoDBtjoqfnz+eDU12l+dj1mz7Y+wAMH2lyLicCTqSgZPNg6fh9zjNVOHXJI2BGZEiWsZmruXDjlFBs52LmzT6WQdB58kHklGjN5xcFcfrnPPeaci0MLFsDf/w6NGsFbb9koqgUL2PjUq1z9cH0uvNAqAWbMgIsvDjvYwvFkKgreeQeuvBI6dYLPP7eaqXhTqxZ89BG8/DJMmGAjCz/9NOyoXFT8+CO88w5vHfMCJUrApZeGHZBzzuXwyy9w+eW2+vCQIXD99Vab3rcvk36rS/PmVgt1zz0wcSI0bBh2wIXnyVQRDR1qn5GTToKRI4s23UFxE4HrroOpU6F2bTjzTHj4YZtp3SWw++8nq1wF3lp6EqecYssPOedc6ObOtau7I46ADz6wWacXLoQXX2Rr1Trcfjt06GAtKBMn2lRRpUqFHfS+8WSqCL76yj4n7dpZrU88J1I5NW0Kkydb7Pffb3271q0LOyq3T6ZPh2HD+KLriyxemsbll4cdkHMu5c2ZY+10TZvaSKxbb7V+Us89B7Vr8/33NlDrmWespW/GDDj++LCDLhpPpvbRvHnQpYvNYP/hh7EfsVdU5cpZP6+XXrLmvmOPtdYil2D+9S844AAGZ17C/vuHP+jBOZfCZs6ECy+0fiSffAJ33gm//WZLd9SsSWYmPPYYtGljy6eNHm0jzxOlk3l+PJnaBxs2wNlnW7PZ6NFQtWrYEe0bEbjhBhvxt3EjHHec9flyCeLLL+HTT9l0y30M/bA0F15oC2U751xMzZgBXbvakPH//tc6P/32Gzz+uK2fhlVWHX883Huv7Tp7NpxxRqhRR5UnU4WkaqPj5s+3/lKHHhp2REV3wgkwZYqt73fGGTZS1cU5VSuVatXig2rXsnkz3sTnnIutadOsieboo2HcOOs38ttv8Mgjf9UybN8ODz1kuyxcaP3P33sPDjgg3NCjrWTYASSal16yfnRPPmkd55JF/fpW0XHhhZYsLlhg34cSnm7Hp7Fjrcdm3768Nmg/jjgi8fscOOcSxLffWoY0ahRUrgwPPmidyytX3m23qVNtCqmZM6F7d+jT56+KqqTj/yoL4dtv4bbbLBG//fawo4m+SpWsmbtXL6ud7dHDripcnMmulWrQgBmtezFliv3NfG4p51yxmjzZmi/atLH7jzxis5jff/9uiVRGhnWXatMGVq2yfsVDhiRvIgVeMxWxLVvgssts2PmgQcn7j6tUKejXzyYdvfNO6yQ4bFhydBBMGiNH2iXfwIG89mYpypSBnj3DDso5l7S+/NJqosaOtYkUn3jC5tmpWHGPXSdOtMnNf/nFaqWeeWaPCquk5DVTEbr7bvtwDByY/B8MEbjjDnjzTeuQ3qkTrFkTdlQOgJ074b77oHFjNp17Kenp1jSbbP0PnHNx4Isv4OST4cQTbU2yp5+2KQ7uvHOPRGr1apvcvEMHa9EYOxbeeCP5/19m82QqAuPGWVvvP/5hn6tUceWVVis1Ywa0bw+//x52RI4hQ2xYzEMPMeSDkmzcCL17hx2Ucy5pqNo/vQ4doGNHm3jzuecsifrnP/dopsjKsgvvxo3h7bctz5o92y7CU4knUwXYvHnXUkJPPBF2NLF37rk2D9WSJdbBed68sCNKYTt2wAMPQIsW0K0b/frZdC5t24YdmHMu4aladVL79vC3v9mQ9RdftCF4t9yS56zUs2dbznXVVdCkCXz/vf2fTLR5F6PBk6kCZPeve+ONxJnhPNpOOsnmotq82aZRmD077IhS1MCBVrA98gjffFuC6dNt9uBk7b/nnIsBVZsb6rjj4NRTbWqDvn1tSPeNN+Y5ed3mzXDXXTbdwY8/2nQ6X3wBzZrFPvx44clUPn780TrPXXGFJeuprGVLmDQJSpa0mt8ZM8KOKMVs3WodQNu1g86def55G33pc0s55/aJqg3fbtPGRugtW2bTkc+fbwsRlymT51Peew8OP9ymB+rZE37+2VpvUn0anRR/+3unumuwwlNPhR1NfDj8cLv6KFfO+o5NnRp2RCnk1Vet09pjj7F4iTBsmE2H4KMsnXOFompzFbRqZUt5rFoFr79ufTh694b99svzad9/b016F19sUxxMmmQ1UtWqxTj+OFVgMiUib4rIShHJs3FHRDqKyHoRmRHc7o9+mLE3ZIglDk8+mdxzYxTWoYfa0NdKlaxZffLksCNKARs32sRfnTpBx4689JI9/I9/hBtWokjVMsy53WRlwfDh1jZ37rmwfr31HP/lF5upuXTpPJ+2apV1J2jZ0vqi9+8P331nXT7cLpHUTA0CTi9gn0mq2iK4PVT0sMKVkWFTIbRsaR3r3O4aNLCEqkYNa2KfODHsiJLciy9aifboo2zaZBeR3brZ8j8uIoNIsTLMub9kZdmyHS1awPnn26SJgwfDTz/ZkO1SpfJ82o4dVvQcdpjlXDfdZJVX11wDaWmxfQuJoMBkSlUnAmtjEEvc6NMHFi+2/lKp3g68N3XrWs3dQQdZc/v//hd2REnqzz/tg9ilC7RuzcCBdkF5yy1hB5Y4UrEMc46dO62D05FH2mR0O3bAO+9YZ+CePa0DbB5UbV7gZs3g5puhdWtbDub551Nnzqh9Ea1UoZ2I/CAin4pI073tJCK9RGSqiExdtWpVlA4dXatWwWOPwTnnWEdrt3cHHggTJsDBB8NZZ8GYMWFHlISefho2bICHHyYzE154wfqgt2kTdmBJJ2nKMJfiMjMhPd2yoYsvtseGDLFh2D167DWJAvjmGxtsdd55Vvv00Uc20O+II2IUewKLRjI1Haivqs2Bl4CRe9tRVfuraitVbVU9TjsiPfSQDft88smwI0kMNWvatAmHH26VJ6NHhx1REvnjD6tn794djjyS//zHZka4446wA0s6SVWGuRSVmWnNd02awKWXWvPdBx/ArFlWhuTTNjdvHlxwgV2ozZsHr71mtVFnn+1Tr0SqyMmUqm5Q1U3B/dFAKRFJyP79v/1m69Jdc40lBy4y1apZM1/TpnZF8/HHYUeUJB5/HLZtgwcfJCvLakyPPNJqTV30JFMZ5lLQjh27piC/4gqbMXP4cJu/plu3fPuqrFxpU0k1aWKTM//73zYzQq9e+VZguTwUOZkSkVoilruKSOvgNRNyJbdHH7Xk/V//CjuSxHPAAbaOX/Pm1sdxxIiwI0pwixdbZn/lldCoESNHWleHe+7xfnzRlkxlmEsh27fbaJTDDrORUlWq2JQH06fbVW0+BcXatVaWHHwwvPKKPX3+fFtgwadb2TcF5p4iMgToCFQTkaXAA0ApAFXtB3QDrhWRTCAD6K6qWmwRF5Nff4VBg+Daa6FOnbCjSUxVqthqBKefbv0dhwyxCyO3Dx5+2H7edx+qNhN/o0ZWFe8KJ1XKMJcitm2z1RAef9wuuo491mYs79y5wDa5DRus3+Wzz9qMKxddZLVRjRvHJPKkVmAypaoXF7C9L9A3ahGFJLtW6s47w44ksVWqZB3RO3e2Zvr0dPvCukKYN88Ky+uvh3r1GD3KJswbONCHJO+LVCnDXJLbutVmyXziCVi61BblfO01OO20ApOozZvh5ZetL/DatVZx9eCD1m3ARYc3GGCdegcNsnZir5Uquv33txEgxx8Pl1xiCZUrhAcesFmI77mHrCy4916rju/RI+zAnHMxl5Fh8/UccgjccINN9Dd2LHz9tTUD5JNIbdpktVCHHGIVBW3a2ISbw4d7IhVt3sUMq5UqWdIWbnTRUaGCjew7+2y47DIbaOLryEVg5kybG+auu6BmTd5/D374waaH2cvces65ZLRli9U8PfWUjezt0MEKgo4dC6yJWrfOaqKefx7WrLHVKoYNswtcVzxSPplatgzefttG8B14YNjRJJfy5W0dzS5drB91ZqbPKF+g++6zqr3bb2fHDhsMcdRRu6aLcc4luU2bbC3OZ56x4XYnn2wXWB06FPjU1autT9RLL1n/qLPOsprttm2LP+xUl/LJVJ8+NlHsrbeGHUlyKlfOJn7r2tWWf8rMtHWeXB6mTLGT9cgjUKUKA/rBggWWkPoIPueS3MaNVp307LOWFZ1yCtx/f0SL4C1darVQ/fpZq+D559tovaOPjkHcDkjxZGrjRvvwde1qbcqueJQta1MldOtmi5JnZlrfapfLv/5lq2rfdBMbN1oH0eOPt878zrkktX69jcZ77jnrHX7GGVZD3a5dgU/9/nvLvf7zH1sG5pJLbF1Zn7E89lI6mRowwD7H//xn2JEkvzJlrM3+oousD2Vmpi2c6QLjx9tEXc89BxUq8Njd1k1ixAifgdi5pLRunTWNPP+83T/rLEuiWrfO92mqNsDnmWdg3Djrn/qPf1h5Wr9+TCJ3eUjZZGrHDvsMn3iir3MWK/vtB++/b/1/br7ZEqrbbgs7qjigah0bDjoIrr2W+fMtp+rZ0/s6OJd01q61ZaJefNGu5rt0sSSqZct8n5aRAe++a2XDjz/ayPOnnrL+vr4AcfhSNpkaOtTmO+vrs8vEVOnS1peyRw+rEczM9Lm9GD0aJk+2kTtlynDrrXaenngi7MCcc1GzZo1dwffpY31Muna1JKpFi3yf9uuv1h99wADLw5o3t0FTF15o5YSLDymbTPXpY7Pwn3lm2JGknlKl7AorezqK7FFrKSl7IqlDDoErr+TTT21tw6eegtq1ww7OOVdkq1ZZx6aXX7bZMy+4wAq8fCZ6ysqyqaT69oVRo2wAyrnnWl/TCGZGcCFIyWRq+nT45hsbQuqjpMJRsqRdXZUsaRdnmZk2V2XKFRJDh/41kdSmbaW49lpb2sH7kzmX4FassI5Nr7xibXTdu9uFU9Ome33K6tVWLr7yiq2VV6OGPeX//s96Abj4lZLJ1Kuv2pB9n0QyXGlpu5ZIefBBS6gefjiFEqrMTBv63LQpdO/Ovbda0/OkSV5971zCWr4cnn7ahopv22ZD7O69Fw4/PM/ds7Js7MmAATBypK1f3K6dlYnnn299TV38S7lk6s8/bXmTSy/1TnvxIC3NCpGSJW0m+vXrrcYwJdage/tt+PlnGD6cyd+m8dJLVo3vsxQ7l4B+/93a5/v3t74Ll15qkz0ddlieuy9aZBeTAwfaRdQBB9jUMVddZRP1usSScsnU4MFW43rddWFH4rKVKGF9rytVsq4FK1ZYnpHUV2TbttmlZ6tWZJx2Lle1grp14bHHwg7MOVcoS5bYCsJvvGEzQPfsaUlUHpMXZmTYvLxvvml9ogA6dbIcrEsXm0LGJaaUSqaysqwtul27AgdQuBgrUcK6F9SubaP8Vq2yKu9KlcKOrJi88YZdmvbvz+13CHPnwpgxULFi2IE55yKyaJENuR0wwH6/4gqbMbNhw912y8y0+aDS023euI0b7cLpvvtsma0GDWIeuSsGKZVMjRsH8+ZZR2cXn267DWrWtEKmQwf49NMkHNW2ZYstGdO+PR9uOYWXX7bljE49NezAnHMF+vVXq0IeNMiuAq++2uZ3yTFjpip8950lUP/5j9W2V6pkA/l69LCyLSW6MqSQlEqmBgyAKlWsU5+LX5deaquqnH8+HHecJVR76buZmPr2hT/+YOnLH/L3q4Sjj/bmPefi3oIF1rHzrbesk2fv3nDHHVbNhCVQM2faSg9DhthovNKlbWLzHj1sWShvxkteKZNM/fmnVbFec41/oBPBaafBhAk2D1jbtnZ1d9ppYUcVBevXw5NPsvXUc+j6RGu2b7eCN6n7hzmXyH75xZKo9HSbJO+GGyyJOvBAsrJgymQYPtxuCxdaZVXHjtbi17WrD3RKFSmTTA0ZYn1+//73sCNxkWrVCr79Fs45x67qnn/e1qBK6KkTnnsOXbuW3qXf5LvvLMFv3DjsoJxze5g715Ko7Kudm26C228ns1otvvgChj9q39/lyy3H6tTJEqhzzrH5oVxqKXDKShF5U0RWisjsvWwXEekjIvNFZKaIHBP9MIvuzTdtGv6jjw47ElcY9evDV1/B2WdbWda7t406TkirV8Nzz/HCUQMZ/ElVHnjAZjV2xStZyjAXI3Pm2AKiTZtatnTbbaz4dhGDj3qWi26qRfXqljgNGmTdENLTbcDM6NHWfcoTqdQUSc3UIKAv8NZetp8BNApubYBXg59xY+ZMmDbN1pV0iadCBatCv/deGzzz44+2vl+dOmFHVkhPPsl7m8/m1plX0LWrzdfpYmIQCV6GuRiYOdMGhgwdys5yFZl6WV9GV7+c0ePLM/Vp26VWLWu6O+ss63ZQrly4Ibv4UWAypaoTRaRBPrt0Ad5SVQW+EZHKIlJbVZdHK8iiGjjQqmEvuSTsSNy+KlECHn/cJrO75hqb2iI9PUFGwKWnw513Mvb3I+jJKNofvoL09Jq+lFGMJEMZ5qIoPd2uzBYvhnr1oFcv9Lup/DpyBuPKnMm4JjMY+0czVr9VghIlrM/mI49YV4MWLRK8m4ErNtHoM1UHWJLj96XBY3sURCLSC+gFUK9evSgcumDbt8M779iEaNWqxeSQrhhdfLE11XbrBqefbuuFPvBAHA8zTk+HXr3435a2nMtIjmAuHy46gzLDnrYhPi4exHUZ5qIo+D6yZQvLqM24RScw7t5ajJPnWUR92Aq11sDpZ1jydOqpULVq2EG7RBCNZCqvPF3z2lFV+wP9AVq1apXnPtE2Zox1VbniilgczcXC4YfDlCm29MrDD8P//mf9Fxo1CjuyPNxzD6O3dKArwzmMX/iMU6mcsdKujD2ZihdxXYa5Itq5E503n/mfLeTru77kq4znmUh7fsbmWzmANZxUZgp3PFOfk0+2ASFe++QKKxrJ1FKgbo7fDwKWReF1o+Ldd+3KIiGag1zEype3BKpTJxvh17y5rehw/fXET/PZqlUMXtyRa3idI5nFZ5xKVdbatsWLw43N5RTXZZgrhLVrYdYsMr6bzbQJG/nqhwp8vawBX2e1YTWNgTOozJ8cx9dcw+uczDia8wMltgLXZYUdvUtg0UimPgJuEJH3sE6b6+Olr8GmTfDhh1YrVapU2NG44nDppXDSSdaP6sYbbcK8l1+2gThh2vnJp9xz0QKeYjB/43OG0o3KrN+1gzcRxZO4LcPcXmRm2vxPM2eybfocZn29kWlzyjBt3cFMoyWzuIYdlAagUeWVnNVsHcd1zOC4N/7OEX+Mo0Tuisd69fM4iHORKzCZEpEhQEegmogsBR4ASgGoaj9gNNAZmA9sAa4srmAL68MPbWHJiy8OOxJXnOrUgVGjbPqL22+3Wqp//MP6UsV8wrwtW/j9uke5fPBJ/I8b6N16On1mnU+pjA279ilXzuavcTGRyGWYw/ppzJwJM2ey7ttfmDNtK7MWlmNaZnOm0ZLZdP0rcapSNoNjGm/h1uO3c9wppWh3nFC9eg0gmK/g8Cuh12Rb0imbfx9dNKhqKLeWLVtqcevcWbVuXdWdO4v9UC5OrFql2quXqohqjRqqL7ygumVLMR3snXdU69e3g9Wvr1kPPazv1blFD2C1liu5VV9/ZbtmZe25n77zTjEFFP+AqRpSmRPtWyzKsKQQ6ed/+3bVWbNU09N1483/0iltb9Q3K9+st/KMnsanWoclaou22O2A8lv1lLYb9K7bM/WDD1QXLlT7vkUrHudyya/8Etsee61atdKpU6cW2+uvXm0L5N56q/Wlcall2jRbNPmLL+xzcPfdNvt9+fJROkCOUUEAs2nKjfRhPCfT6rD1pH9cicMOi9KxkoiITFPVVmHHEQ3FXYYlhVzfE8Bqgp5+mi11GzN/wlLmT1vPvF+ymPfH/szXg5lHI5axaxK5MiV30KRhBk2bl6LZsWVp1gyaNbMl8byjuIul/MqvpF1OZuhQa1b3uaVSU8uWtrbfhAnW3HfjjXDffXD55TaL+hFHFPAC110Hr70GWbk6paalwc6df/2cxjE8zt0MpytV+JNXKt/DNXMeo2TSfrNc3MueR2nRol2f1+zx/WvXWn+97GatnPMtPfroniNMc8/JlNc+uami6zew5qdVLL5xEEu2dGIx9XbdttRj0fX1Wc6Buz2tRrmNNDoog1OPKMWhx+yk6VFpNGsGDRuWIi3NO726+Ja0NVPt28OaNTB7tl+9OFuS5pVX4IMPbDmao46yuao6d7Y+VrslP2lpeyZROfxGfUZyLu9wKdNoRSXWcT0vcyvPUVX+zPe5qS6ZaqZEWmmtWlOpUYM9btWrQ5UqdqtcedetUqVimhOtShVYty7y/UuXttaynGszlSsH/fvvSpZy1CopsJnyrNuvFut63siqWkeyYvE2Vi7PZMUKYcXaUqzYUJYVGRVZub0KK6jBVsrudsj92PpXSlWXJRz69/Yc2rY6jVruz6GHwv77F/ksOFes8iu/kjKZWrbMOiU/9JDVRjiXbcUKmy5j2DD4+mv7f1KxIrRpY7VVh/a7jdo7FlOBTZRmu/0DoTILOZifacxk2rGIBgC0ZCqX8g5XMpBKBB3M69eH334L7f3Fu2RKpmrXbqVnnz2VlSvZ7bZxY/7Pq1hx9+SqfHnLY/K7lSljCX/JkpaM7Xb/wq6U3LwOQcmiBFmUYCdped7fTmkyKLv3W1pFNlZryLqtZVi3QViv+7OOyqynEjv30pCRRiY1Sq6lZtkN1KiQQc0q26lZPYuD6kK9j1+m3vpZ1GMx1Vi9a0Iv/564BJRyzXwjRtjP888PNw4Xf2rWhFtusdvy5danatIk+OYbm7dq445n83yekEV9FtGSafyTZziFsTTml9138lFBKaVOHavIyS0jwxa+XbcO/vzTfua+ZT++fr3d//1361aUfcvIKMyC3sOL9D5KsmNXOrUzg4rrM6lcegsH6mKaMIfKrKMS66nMur/uVx/wJDUaVaJm48ocUK0kJUrkGDGXU/op0OsDHz3nkl5SJlPDhtks2U2ahB2Ji2e1a0P37nYDq6VaWaImq6jOJiqwndJUYBMV2UhdllCGbXu+SHaTYKT9SVzSK1vWPg5FnUpsxw5LqrKTq507rR9oZubu9zPbHv9X7VMaO4M6qKw875fKmTgFt5Ls3HXQnDVGDc6xfle51a8Pf28Y2ZvI/j4Utt+Vcwkm6ZKp1auttuHuu8OOxCUaEajJSmqyMrIn5O5j4lwUlSplt4L7En1d+Bf/q89UjkQqd43Ro4/mPRKvsLVKPXr4d8QlvXhZeCNqPvzQKgq8ic/tk7Jl89+e3Xu4fn1PpFx8iGRm2qpV7SZin90334SBA+1+9mO5P889ethj+e3jnAOSsGZq2DBo2BBatAg7EpeQtmyxq++MjN0fv/ZaGw7oXLz588+8R/OVL2/Te+SX/BSUGHmtknMRSapkat06+PxzuOkmnw7BFUHOZg3nEsGff4YdgXMpLama+T75xDptehOfc84552IlqZKpYcNsuHLr1mFH4pxzzrlUkTTJ1JYtMGYMnHcelEiad+Wcc865eJc0acf//md9hrt0CTsS55xzzqWSpEmmPvnElmlo3z7sSJxzzjmXSpIimVK1ZOq002wuOuecc865WEmKZOr7721x47POCjsS55xzzqWapEimPv7Y5pXq3DnsSJxzzjmXaiJKpkTkdBH5WUTmi8hdeWzvKCLrRWRGcLs/+qHu3SefQNu2UL16LI/qnEsE8V5+OecSX4EzoItIGvAycAqwFPhORD5S1R9z7TpJVWPe0LZsGUydCo89FusjO+fiXbyXX8655BBJzVRrYL6qLlTV7cB7QNxMQDBqlP30/lLOuTzEdfnlnEsOkSRTdYAlOX5fGjyWWzsR+UFEPhWRplGJLgKffGKLmTdrFqsjOucSSFyXX8655BBJMpXXksGa6/fpQH1VbQ68BIzM84VEeonIVBGZumrVqkIFmpetW21h47PO8oWNnXN5ilr5BdEvw5xzySGSZGopUDfH7wcBy3LuoKobVHVTcH80UEpEquV+IVXtr6qtVLVV9Sj0Fp80yZaROeOMIr+Ucy45Ra38CrZHtQxzziWHSJKp74BGItJQREoD3YGPcu4gIrVErG5IRFoHr7sm2sHmNmaMTdLZsWNxH8k5l6DitvxyziWPAkfzqWqmiNwAjAHSgDdVdY6I9A629wO6AdeKSCaQAXRX1dxV6VH32WdwwglQvnxxH8k5l4jiufxyziWPApMp+Kvqe3Sux/rluN8X6Bvd0PK3bBnMmgVPPhnLozrnEk08ll/OueSSsDOgf/aZ/TzttHDjcM4551xqS9hkaswYqFkTjjoq7Eicc845l8oSMpnauRPGjoVTT/UpEZxzzjkXroRMpqZPhzVrvInPOeecc+FLyGQqu7/UKaeEG4dzzjnnXEImU2PGwDHHQI0aYUfinHPOuVSXcMnUhg0webI38TnnnHMuPiRcMjVpEmRmQqdOYUfinHPOOZeAydT48bDfftCuXdiROOecc84lYDI1bhwcdxyULRt2JM4555xzCZZMrVkDM2bASSeFHYlzzjnnnEmoZOqLL0AVTj457Eicc84550xCJVPjx0P58nDssWFH4pxzzjlnEiqZGjcOTjwRSpcOOxLnnHPOOZMwydQff8CPP3oTn3POOefiS8IkU+PH209PppxzzjkXTxIqmapcGVq0CDsS55xzzrldEiaZGjcOOnSAtLSwI3HOOeec2yWiZEpETheRn0Vkvojclcd2EZE+wfaZInJMNINcsgQWLPD5pZxzhRd2+eWcS34FJlMikga8DJwBNAEuFpEmuXY7A2gU3HoBr0YzyIkT7WeHDtF8VedcsouH8ss5l/wiqZlqDcxX1YWquh14D+iSa58uwFtqvgEqi0jtaAU5aRLsvz8ceWS0XtE5lyJCL7+cc8kvkmSqDrAkx+9Lg8cKu88+mzQJjj/e+0s55wot9PLLOZf8Skawj+TxmO7DPohIL6waHWCbiMyO4PiAzTEleR0leqoBq4v1CIUXbzF5PPmLt3gg/mJqHOPjRa38gqKVYQko3j470Zbs7w+S/z3G+v3V39uGSJKppUDdHL8fBCzbh31Q1f5AfwARmaqqrSI4fkzEWzwQfzF5PPmLt3gg/mISkakxPmTUyi+I7zIs2vz9Jb5kf4/x9P4iaeb7DmgkIg1FpDTQHfgo1z4fAT2DUTFtgfWqujzKsTrnXGF5+eWcK3YF1kypaqaI3ACMAdKAN1V1joj0Drb3A0YDnYH5wBbgyuIL2TnnIuPll3MuFiJp5kNVR2MFTs7H+uW4r8D1hTx2/0LuX9ziLR6Iv5g8nvzFWzwQfzHFPJ5iKr8g/s5ttPn7S3zJ/h7j5v2JlSPOOeecc25fJMxyMs4555xz8ShmyZSIXCAic0QkS0T22vu+oKUfohjPASIyVkTmBT+r7GW/30RklojMKI6RSPG41EUEMXUUkfXBOZkhIvcXYyxvisjKvQ1Bj/X5iSCemJ2b4Hh1RWS8iMwNvl835bFPrM9RJDHF9Dztqwj+3j2CczpTRL4WkeaxjrEoCnp/OfY7VkR2iki3WMUWDZG8v+CzOCP4rH4Ry/iiIYLPaCUR+VhEfgjeY0L1CYzHMi5PqhqTG3AENsfMBKDVXvZJAxYABwOlgR+AJsUUz1PAXcH9u4An97Lfb0C1YoqhwPeLdYz9FJsLpy0wpZj/TpHE1BH4JEafm/bAMcDsvWyP9fkpKJ6YnZvgeLWBY4L7FYFf4uAzFElMMT1Pxfj3Pg6oEtw/o7jPbazfX7BPGjAO63fWLeyYo/z3qwz8CNQLfq8RdszF8B7vyf7/BlQH1gKlw467EO8v7sq4vG4xq5lS1bmq+nMBu0Wy9EO0dAEGB/cHA+cW03HyE49LXcTyb1AgVZ2Iffn3JqbnJ4J4YkpVl6vq9OD+RmAue87eHetzFElMCaGgv7eqfq2qfwa/foPNUZUwIvw8/wMYBqws/oiiK4L3dwkwXFUXB/sn43tUoKKICFAh2DczFrFFQzyWcXmJtz5TsVzWoaYGc8kEP2vsZT8FPhORaWKzH0dTPC51Eenx2gXVxp+KSNNijKcg8bgUSCjnRkQaAEcDU3JtCu0c5RMTxM9nKFquwq6Ok4aI1AHOA/oVtG+COgyoIiITgjK+Z9gBFYO+WMvQMmAWcJOqZoUb0r6JxzIuW0RTI0RKRD4HauWx6V5V/TCSl8jjsX0ebphfPIV4meNVdZmI1ADGishPwZVANER1qYsoieR404H6qrpJRDoDI4FGxRhTfmJ9fgoSyrkRkQpY7cHNqroh9+Y8nlLs56iAmOLpM1RkInISlkydEHYsUfYCcKeq7pRiXs8rJCWBlsDfgLLAZBH5RlV/CTesqDoNmAGcDByC/R+blMd3Mq7FYxmXU1STKVXtVMSXiHhZh6LGIyIrRKS2qi4PqgPzrN5V1WXBz5UiMgJrBotWMhXVpS5iFVPOD7KqjhaRV0SkmqqGsQZUrM9PvsI4NyJSCitk0lV1eB67xPwcFRRTnH2GikREjgLeAM5Q1TVhxxNlrYD3gkSqGtBZRDJVdWSoUUXPUmC1qm4GNovIRKA51i8nWVwJPKHWuWi+iPwKHA58G25YkYvHMi63eGvmi2Tph2j5CLg8uH85sEfNmYiUF5GK2feBU4FoLmwaj0tdFBiTiNQK2t8RkdbY5yisfyJxtRRIrM9NcKwBwFxVfW4vu8X0HEUSU5x9hvaZiNQDhgOXJVltBgCq2lBVG6hqA2AocF0SJVJg5f6JIlJSRMoBbbA+OclkMVbzhojUxAaCLQw1okKIxzIuL1GtmcqPiJwHvISNJhglIjNU9TQRORB4Q1U7616WfiimkJ4A3heRq7AP2wVBnH/FA9QERgRlfkngXVX9b7QC2Nv7lRCXuogwpm7AtSKSCWQA3YOrnqgTkSHYyK9qIrIUeAAolSOWmJ6fCOKJ2bkJHA9cBswSkRnBY/cA9XLEFOvlUiKJKdbnaZ9E8Pe+H6gKvBKUE5kaJwuvRiKC95fQCnp/qjpXRP4LzASysLI/mhfMxS6Cv+HDwCARmYU1h92ZYDXA8VjG7cFnQHfOOeecK4J4a+ZzzjnnnEsonkw555xzzhWBJ1POOeecc0XgyZRzzjnnXBF4MuWcc845VwSeTDnnnHPOFYEnU84555xzReDJlHPOOedcEfw/LPBKwF33tX8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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_genres, Value range: 1.00 / 1.00, Num features: 20, feature occurrences: 21564, density 1.01E-01\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_year, Value range: 6.00E+00 / 2.01E+03, Num features: 1, feature occurrences: 10681, density 1.00E+00\n",
      "\tICM name: ICM_tags, Value range: 1.00 / 69.00, Num features: 10106, feature occurrences: 106820, density 9.90E-04\n",
      "\n",
      "\n",
      "Warning: 77 (0.11 %) 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": [
       "7064025"
      ]
     },
     "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": [
       "(61820, 72, 4.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": [
       "4.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 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([   1,    3,    7,    9,   10,   11,   16,   18,   19,   21,   22,\n",
       "         24,   27,   29,   31,   32,   33,   34,   35,   37,   43,   46,\n",
       "         47,   48,   50,   51,   52,   56,   57,   58,   60,   62,   63,\n",
       "         64,   71,   72,   73,   74,   77,   78,   80,   82,   84,   85,\n",
       "         88,   90,   95,   96,   97,   98,  101,  102,  104,  105,  108,\n",
       "        110,  111,  112,  117,  118,  119,  130,  131,  133,  134,  137,\n",
       "        138,  139,  141,  142,  144,  145,  146,  147,  148,  151,  154,\n",
       "        156,  163,  165,  167,  170,  175,  176,  178,  179,  183,  187,\n",
       "        188,  189,  192,  193,  194,  197,  198,  199,  201,  203,  204,\n",
       "        205,  208,  211,  212,  213,  216,  218,  221,  228,  236,  242,\n",
       "        244,  246,  248,  252,  258,  259,  261,  264,  277,  278,  279,\n",
       "        284,  296,  297,  300,  305,  306,  307,  309,  310,  311,  312,\n",
       "        315,  316,  319,  321,  329,  335,  337,  338,  339,  349,  351,\n",
       "        355,  368,  369,  379,  380,  382,  384,  385,  389,  391,  399,\n",
       "        401,  412,  413,  414,  415,  416,  417,  418,  422,  423,  429,\n",
       "        436,  437,  439,  446,  448,  450,  453,  454,  461,  466,  467,\n",
       "        468,  469,  473,  474,  475,  476,  481,  483,  484,  485,  487,\n",
       "        488,  489,  499,  500,  504,  513,  519,  523,  529,  531,  533,\n",
       "        535,  536,  537,  540,  545,  555,  560,  561,  563,  565,  579,\n",
       "        588,  590,  603,  604,  605,  606,  610,  614,  616,  625,  626,\n",
       "        630,  632,  646,  651,  655,  660,  663,  666,  669,  670,  676,\n",
       "        686,  687,  697,  705,  707,  712,  713,  714,  716,  717,  731,\n",
       "        738,  744,  745,  746,  747,  750,  774,  777,  782,  784,  786,\n",
       "        787,  793,  795,  798,  799,  801,  821,  822,  826,  832,  841,\n",
       "        844,  863,  864,  866,  868,  871,  896,  897,  902,  915,  958,\n",
       "        961,  966,  968,  982,  989,  995,  997,  998, 1002, 1005, 1006,\n",
       "       1007, 1008, 1009, 1028, 1032, 1033, 1035, 1036, 1038, 1039, 1041,\n",
       "       1042, 1048, 1051, 1055, 1057, 1059, 1069, 1080, 1081, 1085, 1086,\n",
       "       1093, 1098, 1100, 1102, 1105, 1108, 1109, 1115, 1122, 1131, 1141,\n",
       "       1147, 1153, 1161, 1162, 1171, 1189, 1196, 1204, 1213, 1216, 1229,\n",
       "       1239, 1242, 1251, 1255, 1280, 1281, 1283, 1284, 1287, 1292, 1293,\n",
       "       1294, 1300, 1303, 1304, 1309, 1311, 1313, 1316, 1317, 1319, 1320,\n",
       "       1323, 1325, 1327, 1331, 1333, 1334, 1337, 1344, 1346, 1347, 1351,\n",
       "       1353, 1355, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367,\n",
       "       1368, 1371, 1376, 1385, 1404, 1405, 1412, 1418, 1423, 1425, 1426,\n",
       "       1449, 1464, 1467, 1472, 1473, 1478, 1482, 1484, 1487, 1494, 1500,\n",
       "       1501, 1512, 1522, 1524, 1525, 1526, 1529, 1531, 1535, 1538, 1540,\n",
       "       1544, 1547, 1572, 1573, 1576, 1588, 1593, 1601, 1612, 1613, 1618,\n",
       "       1619, 1623, 1626, 1627, 1634, 1652, 1653, 1654, 1655, 1665, 1666,\n",
       "       1667, 1673, 1674, 1675, 1679, 1682, 1691, 1692, 1695, 1698, 1704,\n",
       "       1705, 1726, 1727, 1728, 1738, 1739, 1741, 1743, 1760, 1762, 1766,\n",
       "       1769, 1772, 1792, 1793, 1795, 1800, 1805, 1812, 1815, 1822, 1823,\n",
       "       1826, 1836, 1848, 1852, 1857, 1861, 1872, 1878, 1879, 1882, 1896,\n",
       "       1906, 1910, 1913, 1919, 1920, 1921, 1922, 1925, 1930, 1931, 1953,\n",
       "       1954, 1977, 1987, 1989, 1991, 1995, 1997, 1999, 2004, 2005, 2006,\n",
       "       2007, 2011, 2015, 2017, 2023, 2025, 2027, 2034, 2037, 2038, 2040,\n",
       "       2045, 2047, 2048, 2049, 2051, 2053, 2055, 2056, 2057, 2058, 2107,\n",
       "       2138, 2158, 2192, 2240, 2257, 2287, 2292, 2302, 2306, 2452, 2472,\n",
       "       2494, 2496, 2498, 2499, 2529, 2536, 2547, 2553, 2554, 2568, 2573,\n",
       "       2598, 2606, 2618, 2653, 2673, 2674, 2677, 2687, 2695, 2712, 2721,\n",
       "       2735, 2782, 2789, 2810, 2823, 2848, 2886, 2894, 2945, 2947, 2966,\n",
       "       2980, 3001, 3007, 3013, 3014, 3021, 3049, 3071, 3073, 3074, 3093,\n",
       "       3094, 3140, 3197, 3204, 3213, 3215, 3222, 3228, 3231, 3233, 3237,\n",
       "       3284, 3310, 3316, 3317, 3329, 3348, 3355, 3387, 3395, 3399, 3421,\n",
       "       3442, 3465, 3489, 3498, 3516, 3570, 3572, 3606, 3626, 3633, 3634,\n",
       "       3669, 3690, 3780, 3790, 3835, 3862, 3870, 3889, 3924, 3952, 3953,\n",
       "       3966, 3970, 4002, 4013, 4030, 4031, 4041, 4052, 4055, 4070, 4074,\n",
       "       4082, 4088, 4181, 4564, 4587, 4617, 5075, 5250, 5283, 5315, 5423,\n",
       "       5442, 5544, 5899, 5901, 6082, 6090, 6170, 6207, 6300, 7025, 7559,\n",
       "       8235], 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([2. , 3. , 3. , 3. , 1. , 4. , 3. , 3. , 3. , 2. , 3. , 4. , 3. ,\n",
       "       3. , 3. , 2. , 3. , 5. , 3. , 3. , 3.5, 4. , 4. , 3. , 4. , 3. ,\n",
       "       4. , 4. , 3. , 3. , 4.5, 4. , 4. , 3. , 3. , 4. , 5. , 3. , 3. ,\n",
       "       3. , 3. , 4. , 3. , 4. , 3. , 3. , 3. , 4. , 3. , 4. , 5. , 1.5,\n",
       "       3. , 5. , 3. , 3. , 3. , 4. , 5. , 4. , 4. , 3. , 4. , 3. , 2. ,\n",
       "       2. , 3. , 4. , 3. , 3. , 3. , 5. , 5. , 5. , 5. , 4. , 3. , 3. ,\n",
       "       5. , 4. , 4. , 3. , 4. , 5. , 5. , 3. , 3. , 4. , 3. , 3. , 3. ,\n",
       "       3. , 2. , 3. , 3. , 3. , 1. , 4. , 2. , 4. , 3. , 3. , 3. , 5. ,\n",
       "       5. , 5. , 4. , 3. , 3. , 4. , 5. , 4. , 2. , 3. , 2. , 4. , 5. ,\n",
       "       4. , 4. , 4. , 4. , 4. , 4. , 3. , 3. , 3. , 3. , 4. , 1. , 1. ,\n",
       "       4. , 4. , 3. , 3. , 2. , 3. , 5. , 2. , 4. , 3. , 4. , 2. , 4. ,\n",
       "       3. , 2. , 3. , 3. , 2. , 4. , 5. , 4. , 3. , 3. , 4. , 4. , 4. ,\n",
       "       4. , 4. , 3. , 3. , 4. , 1. , 5. , 2. , 4. , 1. , 4. , 2. , 3. ,\n",
       "       1. , 1. , 3. , 4. , 5. , 4. , 3. , 3. , 3. , 4. , 2. , 2. , 3. ,\n",
       "       3. , 4. , 3. , 3. , 4. , 2.5, 3. , 3. , 3. , 2. , 3. , 2. , 3. ,\n",
       "       2. , 2. , 3. , 3. , 3. , 3. , 5. , 4. , 3. , 3. , 2. , 3. , 4. ,\n",
       "       4. , 3. , 3. , 5. , 3.5, 5. , 1. , 3. , 2. , 3. , 4. , 3. , 4. ,\n",
       "       3. , 3. , 3. , 4. , 3.5, 4. , 4. , 3. , 3. , 2. , 3. , 4. , 3. ,\n",
       "       3. , 3. , 3. , 3. , 3. , 3. , 3. , 4. , 1. , 3. , 3. , 1. , 3. ,\n",
       "       4. , 3. , 2. , 3. , 2. , 3. , 3. , 4. , 4. , 5. , 2. , 4. , 4. ,\n",
       "       2. , 3.5, 3. , 3.5, 4. , 2. , 2.5, 2. , 4. , 3.5, 4. , 4. , 2.5,\n",
       "       3. , 2.5, 3. , 4. , 4. , 3.5, 3. , 4. , 5. , 4.5, 3. , 3. , 2. ,\n",
       "       5. , 3. , 3. , 5. , 2. , 4. , 4. , 4. , 3. , 3. , 3. , 3. , 3. ,\n",
       "       3. , 4. , 3. , 3. , 3. , 2. , 3. , 3. , 3. , 4. , 1. , 3. , 3. ,\n",
       "       3. , 2. , 4. , 1. , 1. , 3. , 1. , 4. , 3. , 2. , 3. , 3. , 4. ,\n",
       "       4. , 3. , 4. , 2.5, 3.5, 3. , 3. , 3. , 3. , 3. , 3. , 2. , 2. ,\n",
       "       1. , 3. , 5. , 4. , 3. , 3. , 3. , 3. , 3. , 4. , 5. , 4. , 5. ,\n",
       "       3. , 3. , 3. , 2. , 3. , 3. , 3. , 3. , 4. , 2. , 1. , 4. , 5. ,\n",
       "       2. , 2. , 4. , 2. , 3. , 4. , 3. , 3. , 4. , 3. , 3. , 3. , 3. ,\n",
       "       3.5, 3. , 3. , 3. , 3. , 3. , 3. , 3. , 3. , 4. , 3. , 3. , 3. ,\n",
       "       3. , 2. , 4. , 3. , 3. , 4. , 3. , 3. , 3. , 3. , 1. , 1. , 2. ,\n",
       "       4. , 3. , 3. , 3. , 3. , 3. , 4. , 2. , 3. , 3. , 3. , 3. , 2. ,\n",
       "       3. , 2. , 4. , 2. , 2. , 2. , 3. , 3. , 4. , 4. , 3. , 3. , 3. ,\n",
       "       4. , 3. , 3. , 5. , 3. , 3. , 4. , 3. , 3. , 3. , 3. , 3. , 3. ,\n",
       "       3. , 4. , 3. , 3. , 3. , 4. , 2. , 3. , 3. , 4. , 3. , 4. , 4. ,\n",
       "       3. , 1. , 4. , 3. , 4. , 4. , 4. , 4. , 2. , 3.5, 2. , 4. , 3.5,\n",
       "       3. , 4. , 4. , 3. , 3.5, 2. , 3. , 3. , 2.5, 3. , 3. , 3. , 3. ,\n",
       "       3.5, 4. , 4. , 4. , 4. , 3. , 4. , 4. , 4. , 4. , 3. , 2. , 3. ,\n",
       "       4. , 2. , 3. , 4. , 4. , 3. , 3. , 4. , 3. , 3. , 3. , 3. , 3. ,\n",
       "       2. , 4. , 3. , 5. , 3. , 4. , 5. , 4. , 4. , 3. , 4. , 3. , 4. ,\n",
       "       3. , 4. , 3. , 1. , 3. , 4. , 3. , 4.5, 4. , 3. , 5. , 2. , 3. ,\n",
       "       4. , 4. , 3. , 1. , 3. , 3. , 3. , 5. , 3. , 3. , 4. , 4. , 4. ,\n",
       "       5. , 3. , 4. , 4. , 3. , 4. , 2. , 4. , 4. , 3. , 3. , 3. , 2. ,\n",
       "       3. , 3. , 3. , 2. , 3. , 3. , 3. , 4. , 4. , 3. , 3.5, 3. , 1. ,\n",
       "       3. , 5. , 4. , 3. , 3. , 4. , 4. , 4. , 4. , 4. , 3.5, 4.5, 4. ,\n",
       "       3. , 4. , 3. , 4. , 4. , 3. , 4. , 3. , 3. , 3. , 2. , 3. , 4. ,\n",
       "       4. , 3. , 3. , 3. , 3. , 4. , 3.5, 1. , 3. , 1. , 4. , 3. , 3.5,\n",
       "       4. , 3. , 3. , 3. , 3. , 4. , 3. , 4. , 3. , 2. , 2. , 3. , 3. ,\n",
       "       4. , 1. , 3. , 4. , 3. , 3. , 4. , 3. , 4. , 1. , 3.5, 3. , 4. ,\n",
       "       2. , 1. , 1. , 2. , 4. , 4. , 4. , 3. , 2. , 3. , 4. , 3. , 1. ])"
      ]
     },
     "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., 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., 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.])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S[items_in_user_profile,item_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-4\n",
    "regularization_2 = 1e-3\n",
    "\n",
    "# let's simplify and only use the l2 norm for regularization\n",
    "item_item_S[items_in_user_profile,item_id] += learning_rate * (prediction_error * ratings_in_user_profile -\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.0008, 0.0012, 0.0012, 0.0012, 0.0004, 0.0016, 0.0012, 0.0012,\n",
       "       0.0012, 0.0008, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0008,\n",
       "       0.0012, 0.002 , 0.0012, 0.0012, 0.0014, 0.0016, 0.0016, 0.0012,\n",
       "       0.0016, 0.0012, 0.0016, 0.0016, 0.0012, 0.0012, 0.0018, 0.0016,\n",
       "       0.0016, 0.0012, 0.0012, 0.    , 0.002 , 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0016, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0016,\n",
       "       0.0012, 0.0016, 0.002 , 0.0006, 0.0012, 0.002 , 0.0012, 0.0012,\n",
       "       0.0012, 0.0016, 0.002 , 0.0016, 0.0016, 0.0012, 0.0016, 0.0012,\n",
       "       0.0008, 0.0008, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.002 ,\n",
       "       0.002 , 0.002 , 0.002 , 0.0016, 0.0012, 0.0012, 0.002 , 0.0016,\n",
       "       0.0016, 0.0012, 0.0016, 0.002 , 0.002 , 0.0012, 0.0012, 0.0016,\n",
       "       0.0012, 0.0012, 0.0012, 0.0012, 0.0008, 0.0012, 0.0012, 0.0012,\n",
       "       0.0004, 0.0016, 0.0008, 0.0016, 0.0012, 0.0012, 0.0012, 0.002 ,\n",
       "       0.002 , 0.002 , 0.0016, 0.0012, 0.0012, 0.0016, 0.002 , 0.0016,\n",
       "       0.0008, 0.0012, 0.0008, 0.0016, 0.002 , 0.0016, 0.0016, 0.0016,\n",
       "       0.0016, 0.0016, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012, 0.0016,\n",
       "       0.0004, 0.0004, 0.0016, 0.0016, 0.0012, 0.0012, 0.0008, 0.0012,\n",
       "       0.002 , 0.0008, 0.0016, 0.0012, 0.0016, 0.0008, 0.0016, 0.0012,\n",
       "       0.0008, 0.0012, 0.0012, 0.0008, 0.0016, 0.002 , 0.0016, 0.0012,\n",
       "       0.0012, 0.0016, 0.0016, 0.0016, 0.0016, 0.0016, 0.0012, 0.0012,\n",
       "       0.0016, 0.0004, 0.002 , 0.0008, 0.0016, 0.0004, 0.0016, 0.0008,\n",
       "       0.0012, 0.0004, 0.0004, 0.0012, 0.0016, 0.002 , 0.0016, 0.0012,\n",
       "       0.0012, 0.0012, 0.0016, 0.0008, 0.0008, 0.0012, 0.0012, 0.0016,\n",
       "       0.0012, 0.0012, 0.0016, 0.001 , 0.0012, 0.0012, 0.0012, 0.0008,\n",
       "       0.0012, 0.0008, 0.0012, 0.0008, 0.0008, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.002 , 0.0016, 0.0012, 0.0012, 0.0008, 0.0012, 0.0016,\n",
       "       0.0016, 0.0012, 0.0012, 0.002 , 0.0014, 0.002 , 0.0004, 0.0012,\n",
       "       0.0008, 0.0012, 0.0016, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012,\n",
       "       0.0016, 0.0014, 0.0016, 0.0016, 0.0012, 0.0012, 0.0008, 0.0012,\n",
       "       0.0016, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0016, 0.0004, 0.0012, 0.0012, 0.0004, 0.0012, 0.0016,\n",
       "       0.0012, 0.0008, 0.0012, 0.0008, 0.0012, 0.0012, 0.0016, 0.0016,\n",
       "       0.002 , 0.0008, 0.0016, 0.0016, 0.0008, 0.0014, 0.0012, 0.0014,\n",
       "       0.0016, 0.0008, 0.001 , 0.0008, 0.0016, 0.0014, 0.0016, 0.0016,\n",
       "       0.001 , 0.0012, 0.001 , 0.0012, 0.0016, 0.0016, 0.0014, 0.0012,\n",
       "       0.0016, 0.002 , 0.0018, 0.0012, 0.0012, 0.0008, 0.002 , 0.0012,\n",
       "       0.0012, 0.002 , 0.0008, 0.0016, 0.0016, 0.0016, 0.0012, 0.0012,\n",
       "       0.0012, 0.0012, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012,\n",
       "       0.0008, 0.0012, 0.0012, 0.0012, 0.0016, 0.0004, 0.0012, 0.0012,\n",
       "       0.0012, 0.0008, 0.0016, 0.0004, 0.0004, 0.0012, 0.0004, 0.0016,\n",
       "       0.0012, 0.0008, 0.0012, 0.0012, 0.0016, 0.0016, 0.0012, 0.0016,\n",
       "       0.001 , 0.0014, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0008, 0.0008, 0.0004, 0.0012, 0.002 , 0.0016, 0.0012, 0.0012,\n",
       "       0.0012, 0.0012, 0.0012, 0.0016, 0.002 , 0.0016, 0.002 , 0.0012,\n",
       "       0.0012, 0.0012, 0.0008, 0.0012, 0.0012, 0.0012, 0.0012, 0.0016,\n",
       "       0.0008, 0.0004, 0.0016, 0.002 , 0.0008, 0.0008, 0.0016, 0.0008,\n",
       "       0.0012, 0.0016, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0014, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012, 0.0008,\n",
       "       0.0016, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0004, 0.0004, 0.0008, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0016, 0.0008, 0.0012, 0.0012, 0.0012, 0.0012, 0.0008,\n",
       "       0.0012, 0.0008, 0.0016, 0.0008, 0.0008, 0.0008, 0.0012, 0.0012,\n",
       "       0.0016, 0.0016, 0.0012, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012,\n",
       "       0.002 , 0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012, 0.0012, 0.0016,\n",
       "       0.0008, 0.0012, 0.0012, 0.0016, 0.0012, 0.0016, 0.0016, 0.0012,\n",
       "       0.0004, 0.0016, 0.0012, 0.0016, 0.0016, 0.0016, 0.0016, 0.0008,\n",
       "       0.0014, 0.0008, 0.0016, 0.0014, 0.0012, 0.0016, 0.0016, 0.0012,\n",
       "       0.0014, 0.0008, 0.0012, 0.0012, 0.001 , 0.0012, 0.0012, 0.0012,\n",
       "       0.0012, 0.0014, 0.0016, 0.0016, 0.0016, 0.0016, 0.0012, 0.0016,\n",
       "       0.0016, 0.0016, 0.0016, 0.0012, 0.0008, 0.0012, 0.0016, 0.0008,\n",
       "       0.0012, 0.0016, 0.0016, 0.0012, 0.0012, 0.0016, 0.0012, 0.0012,\n",
       "       0.0012, 0.0012, 0.0012, 0.0008, 0.0016, 0.0012, 0.002 , 0.0012,\n",
       "       0.0016, 0.002 , 0.0016, 0.0016, 0.0012, 0.0016, 0.0012, 0.0016,\n",
       "       0.0012, 0.0016, 0.0012, 0.0004, 0.0012, 0.0016, 0.0012, 0.0018,\n",
       "       0.0016, 0.0012, 0.002 , 0.0008, 0.0012, 0.0016, 0.0016, 0.0012,\n",
       "       0.0004, 0.0012, 0.0012, 0.0012, 0.002 , 0.0012, 0.0012, 0.0016,\n",
       "       0.0016, 0.0016, 0.002 , 0.0012, 0.0016, 0.0016, 0.0012, 0.0016,\n",
       "       0.0008, 0.0016, 0.0016, 0.0012, 0.0012, 0.0012, 0.0008, 0.0012,\n",
       "       0.0012, 0.0012, 0.0008, 0.0012, 0.0012, 0.0012, 0.0016, 0.0016,\n",
       "       0.0012, 0.0014, 0.0012, 0.0004, 0.0012, 0.002 , 0.0016, 0.0012,\n",
       "       0.0012, 0.0016, 0.0016, 0.0016, 0.0016, 0.0016, 0.0014, 0.0018,\n",
       "       0.0016, 0.0012, 0.0016, 0.0012, 0.0016, 0.0016, 0.0012, 0.0016,\n",
       "       0.0012, 0.0012, 0.0012, 0.0008, 0.0012, 0.0016, 0.0016, 0.0012,\n",
       "       0.0012, 0.0012, 0.0012, 0.0016, 0.0014, 0.0004, 0.0012, 0.0004,\n",
       "       0.0016, 0.0012, 0.0014, 0.0016, 0.0012, 0.0012, 0.0012, 0.0012,\n",
       "       0.0016, 0.0012, 0.0016, 0.0012, 0.0008, 0.0008, 0.0012, 0.0012,\n",
       "       0.0016, 0.0004, 0.0012, 0.0016, 0.0012, 0.0012, 0.0016, 0.0012,\n",
       "       0.0016, 0.0004, 0.0014, 0.0012, 0.0016, 0.0008, 0.0004, 0.0004,\n",
       "       0.0008, 0.0016, 0.0016, 0.0016, 0.0012, 0.0008, 0.0012, 0.0016,\n",
       "       0.0012, 0.0004])"
      ]
     },
     "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": [
       "2.8873000000000184"
      ]
     },
     "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 3.22 seconds, loss is 13.42. Samples per second 1553.39\n",
      "Iteration 10000 in 6.02 seconds, loss is 13.34. Samples per second 1659.86\n",
      "Iteration 15000 in 8.82 seconds, loss is 13.34. Samples per second 1701.45\n",
      "Iteration 20000 in 11.48 seconds, loss is 13.35. Samples per second 1742.38\n",
      "Iteration 25000 in 14.25 seconds, loss is 13.31. Samples per second 1754.70\n",
      "Iteration 30000 in 17.03 seconds, loss is 13.31. Samples per second 1761.57\n",
      "Iteration 35000 in 19.79 seconds, loss is 13.28. Samples per second 1768.34\n",
      "Iteration 40000 in 22.53 seconds, loss is 13.25. Samples per second 1775.78\n",
      "Iteration 45000 in 25.17 seconds, loss is 13.21. Samples per second 1788.18\n",
      "Iteration 50000 in 27.62 seconds, loss is 13.19. Samples per second 1810.25\n",
      "Iteration 55000 in 30.32 seconds, loss is 13.17. Samples per second 1814.05\n",
      "Iteration 60000 in 33.06 seconds, loss is 13.15. Samples per second 1815.14\n",
      "Iteration 65000 in 35.77 seconds, loss is 13.12. Samples per second 1817.31\n",
      "Iteration 70000 in 38.48 seconds, loss is 13.10. Samples per second 1818.89\n",
      "Iteration 75000 in 41.05 seconds, loss is 13.06. Samples per second 1827.13\n",
      "Iteration 80000 in 43.56 seconds, loss is 13.04. Samples per second 1836.74\n",
      "Iteration 85000 in 46.25 seconds, loss is 13.02. Samples per second 1837.96\n",
      "Iteration 90000 in 48.98 seconds, loss is 12.99. Samples per second 1837.42\n",
      "Iteration 95000 in 51.65 seconds, loss is 12.97. Samples per second 1839.18\n",
      "Iteration 100000 in 54.37 seconds, loss is 12.94. Samples per second 1839.29\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_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_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_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_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 3.13 seconds, loss is 13.40. Samples per second 1597.37\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 2.80 seconds, loss is 13.25. Samples per second 1787.86\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 2.52 seconds, loss is 13.30. Samples per second 1983.19\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 2.82 seconds, loss is 13.28. Samples per second 1773.25\n",
      "\tStopping the epoch early because this implementation is too slow\n",
      "Iteration 5000 in 2.75 seconds, loss is 13.18. Samples per second 1821.49\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_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_2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 4.87044671e-05, 3.26732943e-05, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [1.59282357e-05, 0.00000000e+00, 2.04141814e-04, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [3.74430719e-05, 2.11474715e-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 43495.02 seconds (724.92 minutes, 12.08 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 2.42 seconds, loss is 13.51. Samples per second 2069.01\n",
      "Iteration 10000 in 4.48 seconds, loss is 13.48. Samples per second 2229.87\n",
      "Iteration 15000 in 6.18 seconds, loss is 13.44. Samples per second 2428.33\n",
      "Iteration 20000 in 8.27 seconds, loss is 13.39. Samples per second 2417.81\n",
      "Iteration 25000 in 10.12 seconds, loss is 13.32. Samples per second 2469.33\n",
      "Iteration 30000 in 12.34 seconds, loss is 13.31. Samples per second 2431.83\n",
      "Iteration 35000 in 14.55 seconds, loss is 13.28. Samples per second 2405.12\n",
      "Iteration 40000 in 16.26 seconds, loss is 13.24. Samples per second 2459.41\n",
      "Iteration 45000 in 18.38 seconds, loss is 13.20. Samples per second 2448.57\n",
      "Iteration 50000 in 20.14 seconds, loss is 13.20. Samples per second 2482.22\n",
      "Iteration 55000 in 22.48 seconds, loss is 13.17. Samples per second 2446.70\n",
      "Iteration 60000 in 24.85 seconds, loss is 13.15. Samples per second 2414.96\n",
      "Iteration 65000 in 26.70 seconds, loss is 13.12. Samples per second 2434.16\n",
      "Iteration 70000 in 28.93 seconds, loss is 13.07. Samples per second 2419.59\n",
      "Iteration 75000 in 30.58 seconds, loss is 13.05. Samples per second 2452.57\n",
      "Iteration 80000 in 32.81 seconds, loss is 13.03. Samples per second 2438.03\n",
      "Iteration 85000 in 35.00 seconds, loss is 12.99. Samples per second 2428.35\n",
      "Iteration 90000 in 36.81 seconds, loss is 12.97. Samples per second 2445.31\n",
      "Iteration 95000 in 38.85 seconds, loss is 12.95. Samples per second 2445.23\n",
      "Iteration 100000 in 40.48 seconds, loss is 12.92. Samples per second 2470.11\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_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_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 32387.17 seconds (539.79 minutes, 9.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": [
    "### 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_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_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 2.40 seconds, loss is 13.61. Samples per second 2085.15\n",
      "Iteration 10000 in 4.13 seconds, loss is 13.50. Samples per second 2422.04\n",
      "Iteration 15000 in 6.43 seconds, loss is 13.43. Samples per second 2332.24\n",
      "Iteration 20000 in 8.99 seconds, loss is 13.43. Samples per second 2223.61\n",
      "Iteration 25000 in 11.28 seconds, loss is 13.35. Samples per second 2216.89\n",
      "Iteration 30000 in 13.63 seconds, loss is 13.32. Samples per second 2201.76\n",
      "Iteration 35000 in 15.54 seconds, loss is 13.28. Samples per second 2252.95\n",
      "Iteration 40000 in 17.85 seconds, loss is 13.24. Samples per second 2240.54\n",
      "Iteration 45000 in 19.71 seconds, loss is 13.23. Samples per second 2283.59\n",
      "Iteration 50000 in 21.67 seconds, loss is 13.21. Samples per second 2307.76\n",
      "Iteration 55000 in 23.69 seconds, loss is 13.17. Samples per second 2321.75\n",
      "Iteration 60000 in 25.38 seconds, loss is 13.15. Samples per second 2364.23\n",
      "Iteration 65000 in 27.43 seconds, loss is 13.14. Samples per second 2369.40\n",
      "Iteration 70000 in 29.09 seconds, loss is 13.11. Samples per second 2406.09\n",
      "Iteration 75000 in 31.14 seconds, loss is 13.08. Samples per second 2408.47\n",
      "Iteration 80000 in 33.15 seconds, loss is 13.05. Samples per second 2413.36\n",
      "Iteration 85000 in 34.95 seconds, loss is 13.02. Samples per second 2432.23\n",
      "Iteration 90000 in 36.90 seconds, loss is 12.99. Samples per second 2438.75\n",
      "Iteration 95000 in 38.61 seconds, loss is 12.98. Samples per second 2460.31\n",
      "Iteration 100000 in 40.79 seconds, loss is 12.95. Samples per second 2451.40\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 32634.40 seconds (543.91 minutes, 9.07 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 tiping 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_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_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.06 seconds, loss is 13.44. Samples per second 4699.29\n",
      "Iteration 10000 in 1.24 seconds, loss is 13.35. Samples per second 8091.25\n",
      "Iteration 15000 in 1.37 seconds, loss is 13.32. Samples per second 10988.24\n",
      "Iteration 20000 in 1.49 seconds, loss is 13.30. Samples per second 13466.54\n",
      "Iteration 25000 in 1.62 seconds, loss is 13.32. Samples per second 15434.82\n",
      "Iteration 30000 in 1.73 seconds, loss is 13.29. Samples per second 17295.51\n",
      "Iteration 35000 in 1.84 seconds, loss is 13.27. Samples per second 18999.01\n",
      "Iteration 40000 in 1.95 seconds, loss is 13.26. Samples per second 20532.35\n",
      "Iteration 45000 in 2.05 seconds, loss is 13.24. Samples per second 21991.17\n",
      "Iteration 50000 in 2.15 seconds, loss is 13.20. Samples per second 23225.82\n",
      "Iteration 55000 in 2.25 seconds, loss is 13.17. Samples per second 24439.81\n",
      "Iteration 60000 in 2.35 seconds, loss is 13.15. Samples per second 25569.12\n",
      "Iteration 65000 in 2.45 seconds, loss is 13.13. Samples per second 26537.76\n",
      "Iteration 70000 in 2.55 seconds, loss is 13.10. Samples per second 27468.26\n",
      "Iteration 75000 in 2.66 seconds, loss is 13.08. Samples per second 28191.22\n",
      "Iteration 80000 in 2.80 seconds, loss is 13.05. Samples per second 28569.15\n",
      "Iteration 85000 in 2.90 seconds, loss is 13.02. Samples per second 29296.53\n",
      "Iteration 90000 in 3.03 seconds, loss is 12.99. Samples per second 29734.23\n",
      "Iteration 95000 in 3.16 seconds, loss is 12.96. Samples per second 30018.67\n",
      "Iteration 100000 in 3.28 seconds, loss is 12.93. Samples per second 30445.66\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 2627.63 seconds (43.79 minutes, 0.73 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_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_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 1.35 seconds, loss is 13.19. Samples per second 37018.11\n",
      "Iteration 100000 in 1.95 seconds, loss is 12.93. Samples per second 51217.62\n",
      "Iteration 150000 in 2.51 seconds, loss is 12.70. Samples per second 59681.39\n",
      "Iteration 200000 in 3.07 seconds, loss is 12.47. Samples per second 65108.81\n",
      "Iteration 250000 in 3.60 seconds, loss is 12.25. Samples per second 69408.28\n",
      "Iteration 300000 in 4.18 seconds, loss is 12.06. Samples per second 71718.46\n",
      "Iteration 350000 in 4.83 seconds, loss is 11.87. Samples per second 72490.65\n",
      "Iteration 400000 in 5.38 seconds, loss is 11.68. Samples per second 74334.95\n",
      "Iteration 450000 in 5.95 seconds, loss is 11.50. Samples per second 75675.33\n",
      "Iteration 500000 in 6.47 seconds, loss is 11.34. Samples per second 77283.80\n",
      "Iteration 550000 in 6.98 seconds, loss is 11.18. Samples per second 78774.79\n",
      "Iteration 600000 in 7.53 seconds, loss is 11.03. Samples per second 79656.44\n",
      "Iteration 650000 in 8.20 seconds, loss is 10.88. Samples per second 79306.37\n",
      "Iteration 700000 in 8.79 seconds, loss is 10.74. Samples per second 79647.21\n",
      "Iteration 750000 in 9.36 seconds, loss is 10.61. Samples per second 80143.93\n",
      "Iteration 800000 in 9.92 seconds, loss is 10.49. Samples per second 80604.87\n",
      "Iteration 850000 in 10.46 seconds, loss is 10.36. Samples per second 81246.30\n",
      "Iteration 900000 in 11.05 seconds, loss is 10.25. Samples per second 81456.86\n",
      "Iteration 950000 in 11.75 seconds, loss is 10.13. Samples per second 80864.44\n",
      "Iteration 1000000 in 12.30 seconds, loss is 10.03. Samples per second 81326.37\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 983.69 seconds (16.39 minutes, 0.27 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 10 hours (2k samples per second) and we now have a Cython one with takes only a few minutes (100k 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_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_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_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": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1 complete in in 97.21 seconds, loss is 5.751E+00. Samples per second 82295.58\n",
      "Epoch 2 complete in in 94.24 seconds, loss is 3.707E+00. Samples per second 84889.81\n",
      "Epoch 3 complete in in 95.06 seconds, loss is 3.240E+00. Samples per second 84156.00\n",
      "Epoch 4 complete in in 106.15 seconds, loss is 2.987E+00. Samples per second 75365.70\n",
      "Epoch 5 complete in in 83.18 seconds, loss is 2.821E+00. Samples per second 96174.70\n",
      "Epoch 6 complete in in 86.58 seconds, loss is 2.703E+00. Samples per second 92405.60\n",
      "Epoch 7 complete in in 89.73 seconds, loss is 2.612E+00. Samples per second 89159.01\n",
      "Epoch 8 complete in in 86.16 seconds, loss is 2.536E+00. Samples per second 92856.15\n",
      "Epoch 9 complete in in 87.00 seconds, loss is 2.474E+00. Samples per second 91955.81\n",
      "Epoch 10 complete in in 85.55 seconds, loss is 2.422E+00. Samples per second 93515.19\n"
     ]
    }
   ],
   "source": [
    "n_items = URM_train.shape[1]\n",
    "learning_rate = 1e-6\n",
    "regularization_2 = 1e-3\n",
    "    \n",
    "item_item_S, loss, samples_per_second = train_multiple_epochs(URM_train, learning_rate, regularization_2, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1 complete in in 84.74 seconds, loss is 5.749E+00. Samples per second 94401.75\n",
      "Epoch 2 complete in in 85.36 seconds, loss is 3.708E+00. Samples per second 93725.61\n",
      "Epoch 3 complete in in 86.08 seconds, loss is 3.241E+00. Samples per second 92935.58\n",
      "Epoch 4 complete in in 98.70 seconds, loss is 2.988E+00. Samples per second 81056.75\n",
      "Epoch 5 complete in in 118.07 seconds, loss is 2.824E+00. Samples per second 67755.27\n",
      "Epoch 6 complete in in 118.41 seconds, loss is 2.703E+00. Samples per second 67561.51\n",
      "Epoch 7 complete in in 118.37 seconds, loss is 2.610E+00. Samples per second 67582.51\n",
      "Epoch 8 complete in in 119.30 seconds, loss is 2.534E+00. Samples per second 67056.61\n",
      "Epoch 9 complete in in 118.39 seconds, loss is 2.473E+00. Samples per second 67574.40\n",
      "Epoch 10 complete in in 118.81 seconds, loss is 2.422E+00. Samples per second 67335.62\n"
     ]
    }
   ],
   "source": [
    "item_item_S, loss, samples_per_second = train_multiple_epochs(URM_train, learning_rate, 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",
    "def 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 wose 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 best 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}